axisym_poroelastic_fsi_face_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for elements that are used to integrate fluid tractions
27#ifndef OOMPH_AXISYM_POROELASTIC_FSI_TRACTION_ELEMENTS_HEADER
28#define OOMPH_AXISYM_POROELASTIC_FSI_TRACTION_ELEMENTS_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35
36// OOMPH-LIB headers
37#include "../generic/shape.h"
38#include "../generic/elements.h"
39#include "../generic/element_with_external_element.h"
40
42
43
44namespace oomph
45{
46 //=======================================================================
47 /// Namespace containing the default Strouhal number of axisymmetric
48 /// linearised poroelastic FSI.
49 //=======================================================================
50 namespace LinearisedAxisymPoroelasticBJS_FSIHelper
51 {
52 /// Default for fluid Strouhal number
54
55 /// Default for inverse slip rate coefficient: no slip
57 } // namespace LinearisedAxisymPoroelasticBJS_FSIHelper
58
59
60 //======================================================================
61 /// A class for elements that allow the imposition of the linearised
62 /// poroelastic FSI
63 /// slip condition (according to the Beavers-Joseph-Saffman condition) from an
64 /// adjacent poroelastic axisymmetric medium. The element geometry is obtained
65 /// from the FaceGeometry<ELEMENT> policy class.
66 //======================================================================
67 template<class FLUID_BULK_ELEMENT, class POROELASTICITY_BULK_ELEMENT>
69 : public virtual FaceGeometry<FLUID_BULK_ELEMENT>,
70 public virtual FaceElement,
71 public virtual ElementWithExternalElement
72 {
73 public:
74 /// Constructor, takes the pointer to the "bulk" element and the
75 /// face index identifying the face to which the element is attached.
76 /// The optional identifier can be used
77 /// to distinguish the additional nodal values created by
78 /// this element from thos created by other FaceElements.
80 const int& face_index,
81 const unsigned& id = 0);
82
83 /// Default constructor
85
86 /// Broken copy constructor
88 const LinearisedAxisymPoroelasticBJS_FSIElement& dummy) = delete;
89
90 /// Broken assignment operator
92
93 /// Access function for the pointer to the fluid Strouhal number
94 /// (if not set, St defaults to 1)
95 double*& st_pt()
96 {
97 return St_pt;
98 }
99
100 /// Access function for the fluid Strouhal number
101 double st() const
102 {
103 return *St_pt;
104 }
105
106 /// Inverse slip rate coefficient
108 {
110 }
111
112
113 /// Pointer to inverse slip rate coefficient
115 {
117 }
118
119
120 /// Add the element's contribution to its residual vector
122 {
123 // Call the generic residuals function with flag set to 0
124 // using a dummy matrix argument
127 }
128
129
130 // hieher need to add derivs w.r.t external data (the
131 // bulk velocity dofs
132 /* /// Add the element's contribution to its residual vector and its
133 */
134 /* /// Jacobian matrix */
135 /* void fill_in_contribution_to_jacobian( */
136 /* Vector<double> &residuals, */
137 /* DenseMatrix<double> &jacobian) */
138 /* { */
139 /* //Call the generic routine with the flag set to 1 */
140 /* fill_in_generic_residual_contribution_fpsi_bjs_axisym */
141 /* (residuals,jacobian,1); */
142
143 /* //Derivatives w.r.t. external data */
144 /* fill_in_jacobian_from_external_interaction_by_fd(residuals,jacobian);
145 */
146 /* } */
147
148
149 /// Return this element's contribution to the total volume enclosed
150 /// by collection of these elements
152 {
153 // Initialise
154 double vol = 0.0;
155
156 // Find out how many nodes there are
157 const unsigned n_node = this->nnode();
158
159 // Set up memeory for the shape functions
160 Shape psi(n_node);
161 DShape dpsids(n_node, 1);
162
163 // Set the value of n_intpt
164 const unsigned n_intpt = this->integral_pt()->nweight();
165
166 // Storage for the local coordinate
167 Vector<double> s(1);
168
169 // Loop over the integration points
170 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
171 {
172 // Get the local coordinate at the integration point
173 s[0] = this->integral_pt()->knot(ipt, 0);
174
175 // Get the integral weight
176 double W = this->integral_pt()->weight(ipt);
177
178 // Call the derivatives of the shape function at the knot point
179 this->dshape_local_at_knot(ipt, psi, dpsids);
180
181 // Get position and tangent vector
182 Vector<double> interpolated_t1(2, 0.0);
184 for (unsigned l = 0; l < n_node; l++)
185 {
186 // Loop over directional components
187 for (unsigned i = 0; i < 2; i++)
188 {
189 interpolated_x[i] += this->nodal_position(l, i) * psi(l);
190 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
191 }
192 }
193
194 // Calculate the length of the tangent Vector
195 double tlength = interpolated_t1[0] * interpolated_t1[0] +
196 interpolated_t1[1] * interpolated_t1[1];
197
198 // Set the Jacobian of the line element
199 double J = sqrt(tlength) * interpolated_x[0];
200
201 // Now calculate the normal Vector
202 Vector<double> interpolated_n(2);
203 this->outer_unit_normal(ipt, interpolated_n);
204
205 // Assemble dot product
206 double dot = 0.0;
207 for (unsigned k = 0; k < 2; k++)
208 {
209 dot += interpolated_x[k] * interpolated_n[k];
210 }
211
212 // Add to volume with sign chosen so that...
213
214 // Factor of 1/3 comes from div trick
215 vol += 2.0 * MathematicalConstants::Pi * dot * W * J / 3.0;
216 }
217
218 return vol;
219 }
220
221
222 /// Output function
223 void output(std::ostream& outfile)
224 {
225 // Dummy
226 unsigned nplot = 0;
227 output(outfile, nplot);
228 }
229
230 /// Output function: Output at Gauss points; n_plot is ignored.
231 void output(std::ostream& outfile, const unsigned& n_plot)
232 {
233 // Find out how many nodes there are
234 unsigned n_node = nnode();
235
236 // Get the value of Nintpt
237 const unsigned n_intpt = integral_pt()->nweight();
238
239 // Tecplot header info
240 outfile << this->tecplot_zone_string(n_intpt);
241
242 // Set the Vector to hold local coordinates
243 Vector<double> s(Dim - 1);
244 Vector<double> x_bulk(Dim);
245 Shape psi(n_node);
246 DShape dpsids(n_node, Dim - 1);
247
248 // Cache the Strouhal number
249 const double local_st = st();
250
251 // Cache the slip rate coefficient
252 const double local_inverse_slip_rate_coeff =
254
255 // Loop over the integration points
256 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
257 {
258 // Assign values of s
259 for (unsigned i = 0; i < (Dim - 1); i++)
260 {
261 s[i] = integral_pt()->knot(ipt, i);
262 }
263
264 // Get the outer unit normal
265 Vector<double> interpolated_normal(Dim);
266 outer_unit_normal(ipt, interpolated_normal);
267
268 // Calculate the unit tangent vector
269 Vector<double> interpolated_tangent(Dim);
270 interpolated_tangent[0] = -interpolated_normal[1];
271 interpolated_tangent[1] = interpolated_normal[0];
272
273 // Get solid velocity and porous flux from adjacent solid
274 POROELASTICITY_BULK_ELEMENT* ext_el_pt =
275 dynamic_cast<POROELASTICITY_BULK_ELEMENT*>(
276 external_element_pt(0, ipt));
278 Vector<double> du_dt(3);
279 Vector<double> q(2);
280 ext_el_pt->interpolated_du_dt(s_ext, du_dt);
281 ext_el_pt->interpolated_q(s_ext, q);
282 x_bulk[0] = ext_el_pt->interpolated_x(s_ext, 0);
283 x_bulk[1] = ext_el_pt->interpolated_x(s_ext, 1);
284
285 // Get own coordinates:
287 this->interpolated_x(s, x);
288
289#ifdef PARANOID
291 {
292 double error = sqrt((x[0] - x_bulk[0]) * (x[0] - x_bulk[0]) +
293 (x[1] - x_bulk[1]) * (x[1] - x_bulk[1]));
294 double tol = 1.0e-10;
295 if (error > tol)
296 {
297 std::stringstream junk;
298 junk << "Gap between external and face element coordinate\n"
299 << "is suspiciously large: " << error
300 << "\nBulk/external at: " << x_bulk[0] << " " << x_bulk[1]
301 << "\n"
302 << "Face at: " << x[0] << " " << x[1] << "\n";
304 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
305 }
306 }
307#endif
308
309 // Get permeability from the bulk poroelasticity element
310 const double permeability = ext_el_pt->permeability();
311 const double local_permeability_ratio = ext_el_pt->permeability_ratio();
312
313 // Local coordinate in bulk element
314 Vector<double> s_bulk(Dim);
315 s_bulk = local_coordinate_in_bulk(s);
316
317 // Get the fluid traction from the NSt bulk element
318 Vector<double> traction_nst(3);
319 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
320 ->traction(s_bulk, interpolated_normal, traction_nst);
321
322 // Get fluid velocity from bulk element
323 Vector<double> fluid_veloc(Dim + 1, 0.0);
324 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
325 ->interpolated_u_axi_nst(s_bulk, fluid_veloc);
326
327 // Get fluid pressure from bulk element
328 double p_fluid = dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
329 ->interpolated_p_axi_nst(s_bulk);
330
331 // Calculate the normal components
332 double scaled_normal_wall_veloc = 0.0;
333 double scaled_normal_poro_veloc = 0.0;
334 double scaled_tangential_wall_veloc = 0.0;
335 double scaled_tangential_poro_veloc = 0.0;
336 double normal_nst_veloc = 0.0;
337 for (unsigned i = 0; i < Dim; i++)
338 {
339 scaled_normal_wall_veloc +=
340 local_st * du_dt[i] * interpolated_normal[i];
341
342 scaled_normal_poro_veloc +=
343 local_st * permeability * q[i] * interpolated_normal[i];
344
345 scaled_tangential_wall_veloc +=
346 local_st * du_dt[i] * interpolated_tangent[i];
347
348 scaled_tangential_poro_veloc +=
349 -traction_nst[i] * sqrt(local_permeability_ratio) *
350 local_inverse_slip_rate_coeff * interpolated_tangent[i];
351
352 normal_nst_veloc += fluid_veloc[i] * interpolated_normal[i];
353 }
354
355 // Calculate the combined poroelasticity "velocity" (RHS of BJS BC).
356 double total_poro_normal_component =
357 scaled_normal_wall_veloc + scaled_normal_poro_veloc;
358 double total_poro_tangential_component =
359 scaled_tangential_wall_veloc + scaled_tangential_poro_veloc;
360 Vector<double> poro_veloc(2, 0.0);
361 for (unsigned i = 0; i < Dim; i++)
362 {
363 poro_veloc[i] +=
364 total_poro_normal_component * interpolated_normal[i] +
365 total_poro_tangential_component * interpolated_tangent[i];
366 }
367
368
369 // Call the derivatives of the shape function at the knot point
370 this->dshape_local_at_knot(ipt, psi, dpsids);
371
372 // Get tangent vector
373 Vector<double> interpolated_t1(2, 0.0);
374 for (unsigned l = 0; l < n_node; l++)
375 {
376 // Loop over directional components
377 for (unsigned i = 0; i < 2; i++)
378 {
379 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
380 }
381 }
382
383 // Set the Jacobian of the line element
384 double J = sqrt(1.0 + (interpolated_t1[0] * interpolated_t1[0]) /
385 (interpolated_t1[1] * interpolated_t1[1])) *
386 x[0];
387
388
389 // Default geometry; evaluate everything in deformed (Nst) config.
390 double lagrangian_eulerian_translation_factor = 1.0;
391
392 // Get pointer to associated face element to get geometric information
393 // from (if set up)
395 POROELASTICITY_BULK_ELEMENT,
396 FLUID_BULK_ELEMENT>* ext_face_el_pt =
398 POROELASTICITY_BULK_ELEMENT,
399 FLUID_BULK_ELEMENT>*>(external_element_pt(1, ipt));
400
401 // Update geometry
402 if (ext_face_el_pt != 0)
403 {
405
406 // Get correction factor for geometry
407 lagrangian_eulerian_translation_factor =
408 ext_face_el_pt->lagrangian_eulerian_translation_factor(s_ext_face);
409 }
410
411
412 // Output
413 outfile << x_bulk[0] << " " // column 1
414 << x_bulk[1] << " " // column 2
415 << fluid_veloc[0] << " " // column 3
416 << fluid_veloc[1] << " " // column 4
417 << poro_veloc[0] << " " // column 5
418 << poro_veloc[1] << " " // column 6
419 << normal_nst_veloc * interpolated_normal[0] << " " // column 7
420 << normal_nst_veloc * interpolated_normal[1] << " " // column 8
421 << total_poro_normal_component * interpolated_normal[0]
422 << " " // column 9
423 << total_poro_normal_component * interpolated_normal[1]
424 << " " // column 10
425 << scaled_normal_wall_veloc * interpolated_normal[0]
426 << " " // column 11
427 << scaled_normal_wall_veloc * interpolated_normal[1]
428 << " " // column 12
429 << scaled_normal_poro_veloc * interpolated_normal[0]
430 << " " // column 13
431 << scaled_normal_poro_veloc * interpolated_normal[1]
432 << " " // column 14
433 << p_fluid << " " // column 15
434 << du_dt[0] << " " // column 16
435 << du_dt[1] << " " // column 17
436 << J << " " // column 18
437 << lagrangian_eulerian_translation_factor << " " // column 19
438 << std::endl;
439 }
440 }
441
442
443 /// Compute contributions to integrated porous flux over boundary:
444 /// q_skeleton = \int \partial u_displ / \partial t \cdot n ds
445 /// q_seepage = \int k q \cdot n ds
446 /// q_nst = \int u \cdot n ds
447 void contribution_to_total_porous_flux(double& skeleton_flux_contrib,
448 double& seepage_flux_contrib,
449 double& nst_flux_contrib)
450 {
451 // Get the value of Nintpt
452 const unsigned n_intpt = integral_pt()->nweight();
453
454 // Set the Vector to hold local coordinates
455 Vector<double> s(Dim - 1);
456 Vector<double> x_bulk(Dim);
457
458 // Find out how many nodes there are
459 const unsigned n_node = this->nnode();
460
461 // Set up memeory for the shape functions
462 Shape psi(n_node);
463 DShape dpsids(n_node, 1);
464
465 // Loop over the integration points
466 skeleton_flux_contrib = 0.0;
467 seepage_flux_contrib = 0.0;
468 nst_flux_contrib = 0.0;
469 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
470 {
471 // Assign values of s
472 for (unsigned i = 0; i < (Dim - 1); i++)
473 {
474 s[i] = integral_pt()->knot(ipt, i);
475 }
476
477 // Get the outer unit normal
478 Vector<double> interpolated_normal(Dim);
479 outer_unit_normal(ipt, interpolated_normal);
480
481 // Get the integral weight
482 double W = this->integral_pt()->weight(ipt);
483
484 // Call the derivatives of the shape function at the knot point
485 this->dshape_local_at_knot(ipt, psi, dpsids);
486
487 // Get position and tangent vector
488 Vector<double> interpolated_t1(2, 0.0);
490 for (unsigned l = 0; l < n_node; l++)
491 {
492 // Loop over directional components
493 for (unsigned i = 0; i < 2; i++)
494 {
495 interpolated_x[i] += this->nodal_position(l, i) * psi(l);
496 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
497 }
498 }
499
500 // Calculate the length of the tangent Vector
501 double tlength = interpolated_t1[0] * interpolated_t1[0] +
502 interpolated_t1[1] * interpolated_t1[1];
503
504 // Set the Jacobian of the line element
505 double J = sqrt(tlength) * interpolated_x[0];
506
507 // Get solid velocity and porous flux from adjacent solid
508 POROELASTICITY_BULK_ELEMENT* ext_el_pt =
509 dynamic_cast<POROELASTICITY_BULK_ELEMENT*>(
510 external_element_pt(0, ipt));
512 Vector<double> du_dt(3);
513 Vector<double> q(2);
514 ext_el_pt->interpolated_du_dt(s_ext, du_dt);
515 ext_el_pt->interpolated_q(s_ext, q);
516 x_bulk[0] = ext_el_pt->interpolated_x(s_ext, 0);
517 x_bulk[1] = ext_el_pt->interpolated_x(s_ext, 1);
518
519
520#ifdef PARANOID
522 {
523 // Get own coordinates:
525 this->interpolated_x(s, x);
526
527 double error = sqrt(
528 (interpolated_x[0] - x_bulk[0]) * (interpolated_x[0] - x_bulk[0]) +
529 (interpolated_x[1] - x_bulk[1]) * (interpolated_x[1] - x_bulk[1]));
530 double tol = 1.0e-10;
531 if (error > tol)
532 {
533 std::stringstream junk;
534 junk << "Gap between external and face element coordinate\n"
535 << "is suspiciously large: " << error
536 << "\nBulk/external at: " << x_bulk[0] << " " << x_bulk[1]
537 << "\n"
538 << "Face at: " << x[0] << " " << x[1] << "\n";
540 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
541 }
542 }
543#endif
544
545
546 // Default geometry; evaluate everything in deformed (Nst) config.
547 double lagrangian_eulerian_translation_factor = 1.0;
548
549 // Get the outer unit normal for poro
550 Vector<double> poro_normal(interpolated_normal);
551
552
553 // Get pointer to associated face element to get geometric information
554 // from (if set up)
556 POROELASTICITY_BULK_ELEMENT,
557 FLUID_BULK_ELEMENT>* ext_face_el_pt =
559 POROELASTICITY_BULK_ELEMENT,
560 FLUID_BULK_ELEMENT>*>(external_element_pt(1, ipt));
561
562 // Update geometry
563 if (ext_face_el_pt != 0)
564 {
566
567#ifdef PARANOID
568
569 Vector<double> x_face(2);
570 x_face[0] = ext_face_el_pt->interpolated_x(s_ext_face, 0);
571 x_face[1] = ext_face_el_pt->interpolated_x(s_ext_face, 1);
572
573 double tol = 1.0e-10;
574 double error =
575 std::fabs(x_bulk[0] - x_face[0]) + std::fabs(x_bulk[1] - x_face[1]);
576 if (error > tol)
577 {
578 std::stringstream junk;
579 junk << "Difference in Eulerian coordinates: " << error
580 << " is suspiciously large: "
581 << "Bulk: " << x_bulk[0] << " " << x_bulk[1] << " "
582 << "Face: " << x_face[0] << " " << x_face[1] << "\n";
584 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
585 }
586
587#endif
588
589 // Get correction factor for geometry
590 lagrangian_eulerian_translation_factor =
591 ext_face_el_pt->lagrangian_eulerian_translation_factor(s_ext_face);
592
593 // Get the outer unit normal
594 ext_face_el_pt->outer_unit_normal(s_ext_face, poro_normal);
595 poro_normal[0] = -poro_normal[0];
596 poro_normal[1] = -poro_normal[1];
597 }
598
599 // Get permeability from the bulk poroelasticity element
600 const double permeability = ext_el_pt->permeability();
601
602 // Local coordinate in bulk element
603 Vector<double> s_bulk(Dim);
604 s_bulk = local_coordinate_in_bulk(s);
605
606 // Get fluid velocity from bulk element
607 Vector<double> fluid_veloc(Dim + 1, 0.0);
608 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
609 ->interpolated_u_axi_nst(s_bulk, fluid_veloc);
610
611 // Get net flux through boundary
612 double q_flux = 0.0;
613 double dudt_flux = 0.0;
614 double nst_flux = 0.0;
615 for (unsigned i = 0; i < 2; i++)
616 {
617 q_flux += permeability * q[i] * poro_normal[i];
618 dudt_flux += du_dt[i] * interpolated_normal[i];
619 nst_flux += fluid_veloc[i] * interpolated_normal[i];
620 }
621
622 // Add
623 seepage_flux_contrib += 2.0 * MathematicalConstants::Pi * q_flux *
624 lagrangian_eulerian_translation_factor * W * J;
625 skeleton_flux_contrib +=
626 2.0 * MathematicalConstants::Pi * dudt_flux * W * J;
627 nst_flux_contrib += 2.0 * MathematicalConstants::Pi * nst_flux * W * J;
628 }
629 }
630
631
632 /// C-style output function
633 void output(FILE* file_pt)
634 {
636 }
637
638 /// C-style output function
639 void output(FILE* file_pt, const unsigned& n_plot)
640 {
642 }
643
644
645 protected:
646 /// Function to compute the shape and test functions and to return
647 /// the Jacobian of mapping between local and global (Eulerian)
648 /// coordinates
650 Shape& psi,
651 Shape& test) const
652 {
653 // Find number of nodes
654 unsigned n_node = nnode();
655
656 // Get the shape functions
657 shape(s, psi);
658
659 // Set the test functions to be the same as the shape functions
660 for (unsigned i = 0; i < n_node; i++)
661 {
662 test[i] = psi[i];
663 }
664
665 // Return the value of the jacobian
666 return J_eulerian(s);
667 }
668
669
670 /// Function to compute the shape and test functions and to return
671 /// the Jacobian of mapping between local and global (Eulerian)
672 /// coordinates
673 double shape_and_test_at_knot(const unsigned& ipt,
674 Shape& psi,
675 Shape& test) const
676 {
677 // Find number of nodes
678 unsigned n_node = nnode();
679
680 // Get the shape functions
681 shape_at_knot(ipt, psi);
682
683 // Set the test functions to be the same as the shape functions
684 for (unsigned i = 0; i < n_node; i++)
685 {
686 test[i] = psi[i];
687 }
688
689 // Return the value of the jacobian
690 return J_eulerian_at_knot(ipt);
691 }
692
693 private:
694 /// Add the element's contribution to its residual vector.
695 /// flag=1(or 0): do (or don't) compute the contribution to the
696 /// Jacobian as well.
698 Vector<double>& residuals,
699 DenseMatrix<double>& jacobian,
700 const unsigned& flag);
701
702 /// The spatial dimension of the problem
703 unsigned Dim;
704
705 /// The index at which the velocity unknowns are stored at the nodes
707
708 /// Lagrange Id
709 unsigned Id;
710
711 /// Pointer to fluid Strouhal number
712 double* St_pt;
713
714 /// Pointer to inverse slip rate coefficient
716 };
717
718 /// ///////////////////////////////////////////////////////////////////
719 /// ///////////////////////////////////////////////////////////////////
720 /// ///////////////////////////////////////////////////////////////////
721
722
723 //===========================================================================
724 /// Constructor, takes the pointer to the "bulk" element, and the
725 /// face index that identifies the face of the bulk element to which
726 /// this face element is to be attached.
727 /// The optional identifier can be used
728 /// to distinguish the additional nodal values created by
729 /// this element from thos created by other FaceElements.
730 //===========================================================================
731 template<class FLUID_BULK_ELEMENT, class POROELASTICITY_BULK_ELEMENT>
733 POROELASTICITY_BULK_ELEMENT>::
734 LinearisedAxisymPoroelasticBJS_FSIElement(FiniteElement* const& bulk_el_pt,
735 const int& face_index,
736 const unsigned& id)
737 : FaceGeometry<FLUID_BULK_ELEMENT>(), FaceElement()
738 {
739 // Set source element storage: one interaction with an external element
740 // that provides the velocity of the adjacent linear elasticity
741 // element; one with the associated face element that provides
742 // the geometric normalisation.
743 this->set_ninteraction(2);
744
745 // Store the ID of the FaceElement -- this is used to distinguish
746 // it from any others
747 Id = id;
748
749 // Initialise pointer to fluid Strouhal number. Defaults to 1
751
752 // Initialise pointer to inverse slip rate coefficient. Defaults to 0 (no
753 // slip)
756
757 // Let the bulk element build the FaceElement, i.e. setup the pointers
758 // to its nodes (by referring to the appropriate nodes in the bulk
759 // element), etc.
760 bulk_el_pt->build_face_element(face_index, this);
761
762 // Extract the dimension of the problem from the dimension of
763 // the first node
764 Dim = this->node_pt(0)->ndim();
765
766 // Upcast pointer to bulk element
767 FLUID_BULK_ELEMENT* cast_bulk_el_pt =
768 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_el_pt);
769
770 // Read the index from the (cast) bulk element.
772 for (unsigned i = 0; i < 3; i++)
773 {
774 U_index_axisym_poroelastic_fsi[i] = cast_bulk_el_pt->u_index_axi_nst(i);
775 }
776
777 // The velocities in the bulk affect the shear stress acting
778 // here so we must include them as external data
779 unsigned n = cast_bulk_el_pt->nnode();
780 for (unsigned j = 0; j < n; j++)
781 {
782 Node* nod_pt = cast_bulk_el_pt->node_pt(j);
783 bool do_it = true;
784 unsigned nn = nnode();
785 for (unsigned jj = 0; jj < nn; jj++)
786 {
787 if (nod_pt == node_pt(jj))
788 {
789 do_it = false;
790 break;
791 }
792 }
793 if (do_it) add_external_data(cast_bulk_el_pt->node_pt(j));
794 }
795
796 // We need Dim+1 additional values for each FaceElement node to store the
797 // Lagrange multipliers.
798 Vector<unsigned> n_additional_values(nnode(), Dim + 1);
799
800 // Now add storage for Lagrange multipliers and set the map containing the
801 // position of the first entry of this face element's additional values.
802 add_additional_values(n_additional_values, id);
803 }
804
805 //===========================================================================
806 /// Helper function to compute the element's residual vector and
807 /// the Jacobian matrix.
808 //===========================================================================
809 template<class FLUID_BULK_ELEMENT, class POROELASTICITY_BULK_ELEMENT>
810 void LinearisedAxisymPoroelasticBJS_FSIElement<FLUID_BULK_ELEMENT,
811 POROELASTICITY_BULK_ELEMENT>::
812 fill_in_generic_residual_contribution_axisym_poroelastic_fsi(
813 Vector<double>& residuals,
814 DenseMatrix<double>& jacobian,
815 const unsigned& flag)
816 {
817 // Find out how many nodes there are
818 const unsigned n_node = nnode();
819
820 // Set up memory for the shape and test functions
821 Shape psif(n_node), testf(n_node);
822
823 // Set the value of Nintpt
824 const unsigned n_intpt = integral_pt()->nweight();
825
826 // Set the Vector to hold local coordinates
827 Vector<double> s(Dim - 1);
828
829 // Cache the Strouhal number
830 const double local_st = st();
831
832 // Cache the slip rate coefficient
833 const double local_inverse_slip_rate_coeff =
834 inverse_slip_rate_coefficient();
835
836 // Integers to hold the local equation and unknown numbers
837 int local_eqn = 0;
838
839 // Loop over the integration points
840 // --------------------------------
841 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
842 {
843 // Assign values of s
844 for (unsigned i = 0; i < (Dim - 1); i++)
845 {
846 s[i] = integral_pt()->knot(ipt, i);
847 }
848
849 // Get the integral weight
850 double w = integral_pt()->weight(ipt);
851
852 // Find the shape and test functions and return the Jacobian
853 // of the mapping
854 double J = shape_and_test(s, psif, testf);
855
856 // Calculate the coordinates
857 double interpolated_r = 0;
858
859 // Premultiply the weights and the Jacobian
860 double W = w * J;
861
862 // Calculate the Lagrange multiplier and the fluid veloc
863 Vector<double> lambda(Dim + 1, 0.0);
864 Vector<double> fluid_veloc(Dim + 1, 0.0);
865
866 // Loop over nodes
867 for (unsigned j = 0; j < n_node; j++)
868 {
869 Node* nod_pt = node_pt(j);
870
871 // Cast to a boundary node
872 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(node_pt(j));
873
874 // Get the index of the first nodal value associated with
875 // this FaceElement
876 unsigned first_index =
878
879 // Work out radius
880 interpolated_r += nodal_position(j, 0) * psif(j);
881
882 // Assemble
883 for (unsigned i = 0; i < Dim + 1; i++)
884 {
885 lambda[i] += nod_pt->value(first_index + i) * psif(j);
886 fluid_veloc[i] +=
887 nod_pt->value(U_index_axisym_poroelastic_fsi[i]) * psif(j);
888 }
889 }
890
891 // Local coordinate in bulk element
892 Vector<double> s_bulk(Dim);
893 s_bulk = local_coordinate_in_bulk(s);
894
895#ifdef PARANOID
896 {
897 // Get fluid velocity from bulk element
898 Vector<double> fluid_veloc_from_bulk(Dim + 1, 0.0);
899 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
900 ->interpolated_u_axi_nst(s_bulk, fluid_veloc_from_bulk);
901
902 double error = 0.0;
903 for (unsigned i = 0; i < Dim + 1; i++)
904 {
905 error += (fluid_veloc[i] - fluid_veloc_from_bulk[i]) *
906 (fluid_veloc[i] - fluid_veloc_from_bulk[i]);
907 }
908 error = sqrt(error);
909 double tol = 1.0e-15;
910 if (error > tol)
911 {
912 std::stringstream junk;
913 junk << "Difference in Navier-Stokes velocities\n"
914 << "is suspiciously large: " << error
915 << "\nVeloc from bulk: " << fluid_veloc_from_bulk[0] << " "
916 << fluid_veloc_from_bulk[1] << "\n"
917 << "Veloc from face: " << fluid_veloc[0] << " " << fluid_veloc[1]
918 << "\n";
920 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
921 }
922 }
923#endif
924
925 // Get solid velocity from adjacent solid
926 POROELASTICITY_BULK_ELEMENT* ext_el_pt =
927 dynamic_cast<POROELASTICITY_BULK_ELEMENT*>(external_element_pt(0, ipt));
928 Vector<double> s_ext(external_element_local_coord(0, ipt));
929 Vector<double> du_dt(2), q(2);
930 ext_el_pt->interpolated_du_dt(s_ext, du_dt);
931 ext_el_pt->interpolated_q(s_ext, q);
932
933 // Get the outer unit normal
934 Vector<double> interpolated_normal(Dim);
935 outer_unit_normal(ipt, interpolated_normal);
936
937 // Calculate the unit tangent vector
938 Vector<double> interpolated_tangent(Dim);
939 interpolated_tangent[0] = -interpolated_normal[1];
940 interpolated_tangent[1] = interpolated_normal[0];
941
942 // Normal for poroelastic solid
943 Vector<double> poro_normal(interpolated_normal);
944
945 // Default geometry; evaluate everything in deformed (Nst) config.
946 double lagrangian_eulerian_translation_factor = 1.0;
947
948 // Get pointer to associated face element to get geometric information
949 // from (if set up)
950 FSILinearisedAxisymPoroelasticTractionElement<POROELASTICITY_BULK_ELEMENT,
951 FLUID_BULK_ELEMENT>*
952 ext_face_el_pt =
954 POROELASTICITY_BULK_ELEMENT,
955 FLUID_BULK_ELEMENT>*>(external_element_pt(1, ipt));
956
957 // Update geometry
958 if (ext_face_el_pt != 0)
959 {
960 Vector<double> s_ext_face(external_element_local_coord(1, ipt));
961
962 /* #ifdef PARANOID */
963
964 /* Vector<double> x_face(2); */
965 /* x_face[0]=ext_face_el_pt->interpolated_x(s_ext_face,0); */
966 /* x_face[1]=ext_face_el_pt->interpolated_x(s_ext_face,1); */
967
968 /* Vector<double> x_bulk(2); */
969 /* x_bulk[0]=this->interpolated_x(s,0); */
970 /* x_bulk[1]=this->interpolated_x(s,1); */
971
972 /* double tol=1.0e-10; */
973 /* double
974 * error=std::fabs(x_bulk[0]-x_face[0])+std::fabs(x_bulk[1]-x_face[1]);
975 */
976 /* if (error>tol) */
977 /* { */
978 /* std::stringstream junk; */
979 /* junk */
980 /* << "Difference in Eulerian coordinates: " << error */
981 /* << " is suspiciously large: " */
982 /* << "Bulk: " << x_bulk[0] << " " << x_bulk[1] << " " */
983 /* << "Face: " << x_face[0] << " " << x_face[1] << "\n"; */
984 /* OomphLibWarning(junk.str(), */
985 /* OOMPH_CURRENT_FUNCTION, */
986 /* OOMPH_EXCEPTION_LOCATION); */
987 /* } */
988
989 /* #endif */
990
991 // Get correction factor for geometry
992 lagrangian_eulerian_translation_factor =
993 ext_face_el_pt->lagrangian_eulerian_translation_factor(s_ext_face);
994
995 // Get the outer unit normal
996 ext_face_el_pt->outer_unit_normal(s_ext_face, poro_normal);
997 poro_normal[0] = -poro_normal[0];
998 poro_normal[1] = -poro_normal[1];
999 }
1000
1001
1002 // Get permeability from the bulk poroelasticity element
1003 const double permeability = ext_el_pt->permeability();
1004 const double local_permeability_ratio = ext_el_pt->permeability_ratio();
1005
1006 // We are given the normal and tangential components of the combined
1007 // poroelasticity "velocity" at the boundary from the BJS condition ---
1008 // calculate the vector in r-z coords from these
1009 double poro_normal_component = 0.0;
1010 double poro_tangential_component = 0.0;
1011
1012 // Get the fluid traction from the NSt bulk element
1013 Vector<double> traction_nst(3);
1014 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_element_pt())
1015 ->traction(s_bulk, interpolated_normal, traction_nst);
1016
1017 // Calculate the normal and tangential components
1018 for (unsigned i = 0; i < Dim; i++)
1019 {
1020 // Normal component computed with scaling factor for mass conservation
1021 poro_normal_component +=
1022 local_st *
1023 (du_dt[i] * interpolated_normal[i] +
1024 permeability * q[i] * lagrangian_eulerian_translation_factor *
1025 poro_normal[i]);
1026
1027 // Leave this one alone... There isn't really much point
1028 // in trying to correct this.
1029 poro_tangential_component +=
1030 (local_st * du_dt[i] - traction_nst[i] *
1031 sqrt(local_permeability_ratio) *
1032 local_inverse_slip_rate_coeff) *
1033 interpolated_tangent[i];
1034 }
1035
1036 // Get the normal and tangential nst components
1037 double nst_normal_component = 0.0;
1038 double nst_tangential_component = 0.0;
1039 for (unsigned i = 0; i < Dim; i++)
1040 {
1041 nst_normal_component += fluid_veloc[i] * interpolated_normal[i];
1042 nst_tangential_component += fluid_veloc[i] * interpolated_tangent[i];
1043 }
1044
1045 // Setup bjs terms
1046 Vector<double> bjs_term(2);
1047 bjs_term[0] = nst_normal_component - poro_normal_component;
1048 bjs_term[1] = nst_tangential_component - poro_tangential_component;
1049
1050 // Now add to the appropriate equations
1051
1052 // Loop over the test functions
1053 for (unsigned l = 0; l < n_node; l++)
1054 {
1055 // Loop over directions
1056 for (unsigned i = 0; i < Dim + 1; i++)
1057 {
1058 // Add contribution to bulk Navier Stokes equations where
1059 // ------------------------------------------------------
1060 // the Lagrange multiplier acts as a traction
1061 // ------------------------------------------
1062 local_eqn = nodal_local_eqn(l, U_index_axisym_poroelastic_fsi[i]);
1063
1064 /*IF it's not a boundary condition*/
1065 if (local_eqn >= 0)
1066 {
1067 // Add the Lagrange multiplier "traction" to the bulk
1068 residuals[local_eqn] -= lambda[i] * testf[l] * interpolated_r * W;
1069
1070 // Jacobian entries
1071 if (flag)
1072 {
1073 throw OomphLibError("Jacobian not done yet",
1074 OOMPH_CURRENT_FUNCTION,
1075 OOMPH_EXCEPTION_LOCATION);
1076
1077 // Loop over the lagrange multiplier unknowns
1078 for (unsigned l2 = 0; l2 < n_node; l2++)
1079 {
1080 // Cast to a boundary node
1081 BoundaryNodeBase* bnod_pt =
1082 dynamic_cast<BoundaryNodeBase*>(node_pt(l2));
1083
1084 // Local unknown
1085 int local_unknown = nodal_local_eqn(
1086 l2,
1088 i);
1089
1090 // If it's not pinned
1091 if (local_unknown >= 0)
1092 {
1093 jacobian(local_eqn, local_unknown) -=
1094 psif[l2] * testf[l] * interpolated_r * W;
1095 }
1096 }
1097 }
1098 }
1099
1100 // Now do the Lagrange multiplier equations
1101 //-----------------------------------------
1102 // Cast to a boundary node
1103 BoundaryNodeBase* bnod_pt =
1104 dynamic_cast<BoundaryNodeBase*>(node_pt(l));
1105
1106 // Local eqn number:
1107 int local_eqn = nodal_local_eqn(
1109
1110 // If it's not pinned
1111 if (local_eqn >= 0)
1112 {
1113#ifdef PARANOID
1114 if (i == Dim)
1115 {
1116 std::stringstream junk;
1117 junk << "Elements have not been validated for nonzero swirl!\n";
1119 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1120 }
1121#endif
1122
1123 residuals[local_eqn] += bjs_term[i] * testf(l) * interpolated_r * W;
1124
1125 // Jacobian entries
1126 if (flag)
1127 {
1128 throw OomphLibError("Jacobian not done yet",
1129 OOMPH_CURRENT_FUNCTION,
1130 OOMPH_EXCEPTION_LOCATION);
1131
1132 // Loop over the velocity unknowns [derivs w.r.t. to
1133 // wall velocity taken care of by fd-ing
1134 for (unsigned l2 = 0; l2 < n_node; l2++)
1135 {
1136 int local_unknown =
1137 nodal_local_eqn(l2, U_index_axisym_poroelastic_fsi[i]);
1138
1139 /*IF it's not a boundary condition*/
1140 if (local_unknown >= 0)
1141 {
1142 jacobian(local_eqn, local_unknown) -=
1143 psif[l2] * testf[l] * interpolated_r * W;
1144 }
1145 }
1146 }
1147 }
1148 }
1149 }
1150 }
1151 }
1152
1153} // namespace oomph
1154
1155
1156#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Return the index of the first value associated with the i-th face element value. If no argument is sp...
Definition: nodes.h:2061
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement.
Definition: elements.cc:6353
void add_additional_values(const Vector< unsigned > &nadditional_values, const unsigned &id)
Helper function adding additional values for the unknowns associated with the FaceElement....
Definition: elements.h:4428
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
Definition: elements.cc:5242
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5328
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3239
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
A class for elements that allow the imposition of the linearised poroelastic FSI slip condition (acco...
double inverse_slip_rate_coefficient() const
Inverse slip rate coefficient.
double *& st_pt()
Access function for the pointer to the fluid Strouhal number (if not set, St defaults to 1)
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Output at Gauss points; n_plot is ignored.
double st() const
Access function for the fluid Strouhal number.
LinearisedAxisymPoroelasticBJS_FSIElement(const LinearisedAxisymPoroelasticBJS_FSIElement &dummy)=delete
Broken copy constructor.
double * Inverse_slip_rate_coeff_pt
Pointer to inverse slip rate coefficient.
void fill_in_generic_residual_contribution_axisym_poroelastic_fsi(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed by collection of these elements.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void operator=(const LinearisedAxisymPoroelasticBJS_FSIElement &)=delete
Broken assignment operator.
Vector< unsigned > U_index_axisym_poroelastic_fsi
The index at which the velocity unknowns are stored at the nodes.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
double *& inverse_slip_rate_coefficient_pt()
Pointer to inverse slip rate coefficient.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void contribution_to_total_porous_flux(double &skeleton_flux_contrib, double &seepage_flux_contrib, double &nst_flux_contrib)
Compute contributions to integrated porous flux over boundary: q_skeleton = \int \partial u_displ / \...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
bool Allow_gap_in_FSI
Public boolean to allow gap between poro-elastic and Navier Stokes element in FSI computations....
void output()
Doc the command line arguments.
double Default_strouhal_number
Default for fluid Strouhal number.
double Default_inverse_slip_rate_coefficient
Default for inverse slip rate coefficient: no slip.
const double Pi
50 digits from maple
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:291
//////////////////////////////////////////////////////////////////// ////////////////////////////////...