axisym_poroelasticity_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 apply surface loads to
27// the Darcy equations
28
29#ifndef OOMPH_AXISYM_POROELASITICTY_FACE_ELEMENTS_HEADER
30#define OOMPH_AXISYM_POROELASITICTY_FACE_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "generic/Qelements.h"
41
42namespace oomph
43{
44 //=======================================================================
45 /// Namespace containing the zero pressure function for Darcy pressure
46 /// elements
47 //=======================================================================
48 namespace AxisymmetricPoroelasticityTractionElementHelper
49 {
50 //=======================================================================
51 /// Default load function (zero traction)
52 //=======================================================================
53 void Zero_traction_fct(const double& time,
54 const Vector<double>& x,
55 const Vector<double>& N,
56 Vector<double>& load)
57 {
58 unsigned n_dim = load.size();
59 for (unsigned i = 0; i < n_dim; i++)
60 {
61 load[i] = 0.0;
62 }
63 }
64
65 //=======================================================================
66 /// Default load function (zero pressure)
67 //=======================================================================
68 void Zero_pressure_fct(const double& time,
69 const Vector<double>& x,
70 const Vector<double>& N,
71 double& load)
72 {
73 load = 0.0;
74 }
75
76 /// Public boolean to allow gap between poro-elastic and Navier
77 /// Stokes element in FSI computations. Useful in hybrid linear/nonlinear
78 /// geometry runs where this will happen
79 bool Allow_gap_in_FSI = false;
80
81 } // namespace AxisymmetricPoroelasticityTractionElementHelper
82
83
84 //======================================================================
85 /// A class for elements that allow the imposition of an applied combined
86 /// traction and pore fluid pressure in the axisym poroelasticity equations.
87 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
88 /// class and thus, we can be generic enough without the need to have
89 /// a separate equations class.
90 //======================================================================
91 template<class ELEMENT>
93 : public virtual FaceGeometry<ELEMENT>,
94 public virtual FaceElement
95 {
96 protected:
97 /// Pointer to an imposed traction function. Arguments:
98 /// Eulerian coordinate; outer unit normal; applied traction.
99 /// (Not all of the input arguments will be required for all specific load
100 /// functions but the list should cover all cases)
101 void (*Traction_fct_pt)(const double& time,
102 const Vector<double>& x,
103 const Vector<double>& n,
104 Vector<double>& result);
105
106 /// Pointer to an imposed pressure function. Arguments:
107 /// Eulerian coordinate; outer unit normal; applied pressure.
108 /// (Not all of the input arguments will be required for all specific load
109 /// functions but the list should cover all cases)
110 void (*Pressure_fct_pt)(const double& time,
111 const Vector<double>& x,
112 const Vector<double>& n,
113 double& result);
114
115
116 /// Get the traction vector: Pass number of integration point
117 /// (dummy), Eulerrian coordinate and normal vector and return the pressure
118 /// (not all of the input arguments will be required for all specific load
119 /// functions but the list should cover all cases). This function is virtual
120 /// so it can be overloaded for FSI.
121 virtual void get_traction(const double& time,
122 const unsigned& intpt,
123 const Vector<double>& x,
124 const Vector<double>& n,
126 {
127 Traction_fct_pt(time, x, n, traction);
128 }
129
130 /// Get the pressure value: Pass number of integration point (dummy),
131 /// Eulerrian coordinate and normal vector and return the pressure
132 /// (not all of the input arguments will be required for all specific load
133 /// functions but the list should cover all cases). This function is virtual
134 /// so it can be overloaded for FSI.
135 virtual void get_pressure(const double& time,
136 const unsigned& intpt,
137 const Vector<double>& x,
138 const Vector<double>& n,
139 double& pressure)
140 {
141 Pressure_fct_pt(time, x, n, pressure);
142 }
143
144
145 /// Helper function that actually calculates the residuals
146 // This small level of indirection is required to avoid calling
147 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
148 // which causes all kinds of pain if overloading later on
150 Vector<double>& residuals);
151
152
153 public:
154 /// Constructor, which takes a "bulk" element and the value of the
155 /// index and its limit
157 const int& face_index)
158 : FaceGeometry<ELEMENT>(), FaceElement()
159 {
160#ifdef PARANOID
161 {
162 // Check that the element is not a refineable 3d element
163 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
164 // If it's three-d
165 if (elem_pt->dim() == 3)
166 {
167 // Is it refineable
168 RefineableElement* ref_el_pt =
169 dynamic_cast<RefineableElement*>(elem_pt);
170 if (ref_el_pt != 0)
171 {
172 if (this->has_hanging_nodes())
173 {
174 throw OomphLibError("This flux element will not work correctly "
175 "if nodes are hanging\n",
176 OOMPH_CURRENT_FUNCTION,
177 OOMPH_EXCEPTION_LOCATION);
178 }
179 }
180 }
181 }
182#endif
183
184 // Attach the geometrical information to the element. N.B. This function
185 // also assigns nbulk_value from the required_nvalue of the bulk element
186 element_pt->build_face_element(face_index, this);
187
188 // Zero traction
191
192 // Zero pressure
195 }
196
197
198 /// Default constructor
200
201 /// Reference to the traction function pointer
202 void (*&traction_fct_pt())(const double& time,
203 const Vector<double>& x,
204 const Vector<double>& n,
206 {
207 return Traction_fct_pt;
208 }
209
210
211 /// Reference to the pressure function pointer
212 void (*&pressure_fct_pt())(const double& time,
213 const Vector<double>& x,
214 const Vector<double>& n,
215 double& pressure)
216 {
217 return Pressure_fct_pt;
218 }
219
220
221 /// Return the residuals
223 {
225 }
226
227
228 /// Fill in contribution from Jacobian
230 DenseMatrix<double>& jacobian)
231 {
232 // Call the residuals (element makes no contribution to Jacobian)
234 }
235
236 /// Specify the value of nodal zeta from the face geometry
237 /// The "global" intrinsic coordinate of the element when
238 /// viewed as part of a geometric object should be given by
239 /// the FaceElement representation, by default (needed to break
240 /// indeterminacy if bulk element is SolidElement)
241 double zeta_nodal(const unsigned& n,
242 const unsigned& k,
243 const unsigned& i) const
244 {
245 return FaceElement::zeta_nodal(n, k, i);
246 }
247
248 /// Output function
249 void output(std::ostream& outfile)
250 {
251 unsigned n_plot = 5;
252 output(outfile, n_plot);
253 }
254
255
256 /// Output function
257 void output(std::ostream& outfile, const unsigned& n_plot)
258 {
259 // Get continuous time from timestepper of first node
260 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
261 unsigned n_dim = this->nodal_dimension();
262
263 // Find out how many nodes there are
264 unsigned n_node = nnode();
265
266 Vector<double> x(n_dim);
267 Vector<double> s(n_dim - 1);
268 Vector<double> s_bulk(n_dim);
269 Vector<double> disp(n_dim);
270 Shape psi(n_node);
271 DShape dpsids(n_node, n_dim - 1);
272
273 // Tecplot header info
274 outfile << this->tecplot_zone_string(n_plot);
275
276 // Loop over plot points
277 unsigned num_plot_points = this->nplot_points(n_plot);
278 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
279 {
280 // Get local coordinates of plot point
281 this->get_s_plot(iplot, n_plot, s);
282
283 // Call the derivatives of the shape function at the knot point
284 this->dshape_local(s, psi, dpsids);
285
286 // Get pointer to bulk element
287 ELEMENT* bulk_pt = dynamic_cast<ELEMENT*>(bulk_element_pt());
288 s_bulk = local_coordinate_in_bulk(s);
289
290
291 // Get Eulerian coordinates
292 this->interpolated_x(s, x);
293
294 // Outer unit normal
295 Vector<double> unit_normal(n_dim);
296 outer_unit_normal(s, unit_normal);
297
298 /// Calculate the FE representation of u -- the skeleton displacement
299 bulk_pt->interpolated_u(s_bulk, disp);
300
301 /// Skeleton velocity
302 Vector<double> du_dt(2);
303 bulk_pt->interpolated_du_dt(s_bulk, du_dt);
304
305 // Porous seepage flux
306 Vector<double> q(2);
307 bulk_pt->interpolated_q(s_bulk, q);
308
309 // Get permeability from the bulk poroelasticity element
310 const double permeability = bulk_pt->permeability();
311
312
313 // Surface area: S = 2 \pi \int r \sqrt((dr/ds)^2+(dz/ds)^2) ds
314 // = 2 \pi \int r \sqrt( 1 + ( (dr/ds)/(dz/ds) )^2 )
315 // dz/ds ds = 2 \pi \int J dz
316 // where J is an objective measure of the length of the line element
317 // (indep of local coordinates) so should be the same from fluid
318 // and solid.
319
320 // Get deformed and undeformed tangent vectors
321 Vector<double> interpolated_t1(2, 0.0);
322 Vector<double> interpolated_T1(2, 0.0);
323 for (unsigned l = 0; l < n_node; l++)
324 {
325 // Loop over directional components
326 for (unsigned i = 0; i < 2; i++)
327 {
328 // Index at which the nodal value is stored
329 unsigned u_nodal_index = bulk_pt->u_index_axisym_poroelasticity(i);
330
331 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
332 interpolated_T1[i] +=
333 (this->nodal_position(l, i) + nodal_value(l, u_nodal_index)) *
334 dpsids(l, 0);
335 }
336 }
337
338 // Set the Jacobian of the undeformed line element
339 double J_undef =
340 sqrt(1.0 + (interpolated_t1[0] * interpolated_t1[0]) /
341 (interpolated_t1[1] * interpolated_t1[1])) *
342 x[0];
343
344
345 // Set the Jacobian of the deformed line element
346 double J_def = sqrt(1.0 + (interpolated_T1[0] * interpolated_T1[0]) /
347 (interpolated_T1[1] * interpolated_T1[1])) *
348 (x[0] + disp[0]);
349
350 // Dummy
351 unsigned ipt = 0;
352
353 // Now calculate the traction load
355 get_traction(time, ipt, x, unit_normal, traction);
356
357 // Now calculate the load
358 double pressure;
359 get_pressure(time, ipt, x, unit_normal, pressure);
360
361 // Get correction factor for geometry
362 double lagr_euler_translation_factor =
364
365 // Output the x,y,..
366 for (unsigned i = 0; i < n_dim; i++)
367 {
368 outfile << x[i] << " ";
369 } // column 1,2
370
371 // Output displacement
372 for (unsigned i = 0; i < n_dim; i++)
373 {
374 outfile << disp[i] << " "; // column 3,4
375 }
376
377 // Output imposed traction
378 for (unsigned i = 0; i < n_dim; i++)
379 {
380 outfile << traction[i] << " "; // column 5,6
381 }
382
383 // Output imposed pressure
384 outfile << pressure << " "; // column 7
385
386 // Output seepage flux
387 outfile << permeability * q[0] << " " // column 8
388 << permeability * q[1] << " "; // column 9
389
390 // Output skeleton velocity
391 outfile << du_dt[0] << " " // column 10
392 << du_dt[1] << " "; // column 11
393
394 // Total veloc
395 outfile << du_dt[0] + permeability * q[0] << " " // column 12
396 << du_dt[1] + permeability * q[1] << " "; // column 13
397
398 // Outer unit normal
399 outfile << unit_normal[0] << " " // column 14
400 << unit_normal[1] << " "; // column 15
401
402 // Output FE representation of div u at s_bulk
403 outfile << bulk_pt->interpolated_div_q(s_bulk) << " "; // column 16
404
405 // Output FE representation of p at s_bulk
406 outfile << bulk_pt->interpolated_p(s_bulk) << " "; // column 17
407
408 // Output undeformed jacobian
409 outfile << J_undef << " "; // column 18
410
411 // Output deformed jacobian
412 outfile << J_def << " "; // column 19
413
414 // Lagranian/Eulerian translation factor
415 outfile << lagr_euler_translation_factor << " "; // column 20
416
417 outfile << std::endl;
418 }
419
420 // Write tecplot footer (e.g. FE connectivity lists)
421 this->write_tecplot_zone_footer(outfile, n_plot);
422 }
423
424 /// C_style output function
425 void output(FILE* file_pt)
426 {
428 }
429
430 /// C-style output function
431 void output(FILE* file_pt, const unsigned& n_plot)
432 {
433 FaceGeometry<ELEMENT>::output(file_pt, n_plot);
434 }
435
436
437 /// Ratio of lengths of line elements (or annular surface areas)
438 /// in the undeformed and deformed configuration (needed to translate
439 /// normal flux correctly between small displacement formulation used
440 /// here and the possibly large displacement NSt formulation used in
441 /// FSI problems. Maths is as follows:
442 /// Surface area: A = 2 \pi \int r \sqrt((dr/ds)^2+(dz/ds)^2) ds
443 /// = 2 \pi \int J ds
444 /// This function returns the ratio of J in the undeformed
445 /// configuration (used here) to that in the deformed configuration
446 /// where (r,z) = (r,z)_undef + (u_r,u_z).
448 {
449 // Get continuous time from timestepper of first node
450 unsigned n_dim = this->nodal_dimension();
451
452 // Find out how many nodes there are
453 unsigned n_node = nnode();
454
455 Vector<double> x(n_dim);
456 Vector<double> s_bulk(n_dim);
457 Vector<double> disp(n_dim);
458 Shape psi(n_node);
459 DShape dpsids(n_node, n_dim - 1);
460
461 // Call the derivatives of the shape function at the knot point
462 this->dshape_local(s, psi, dpsids);
463
464 // Get pointer to bulk element
465 ELEMENT* bulk_pt = dynamic_cast<ELEMENT*>(bulk_element_pt());
466 s_bulk = local_coordinate_in_bulk(s);
467
468 // Get Eulerian coordinates
469 this->interpolated_x(s, x);
470
471 // Outer unit normal
472 Vector<double> unit_normal(n_dim);
473 outer_unit_normal(s, unit_normal);
474
475 /// Calculate the FE representation of u -- the skeleton displacement
476 bulk_pt->interpolated_u(s_bulk, disp);
477
478 // Get deformed and undeformed tangent vectors
479 Vector<double> interpolated_t1(2, 0.0);
480 Vector<double> interpolated_T1(2, 0.0);
481 for (unsigned l = 0; l < n_node; l++)
482 {
483 // Loop over directional components
484 for (unsigned i = 0; i < 2; i++)
485 {
486 // Index at which the nodal value is stored
487 unsigned u_nodal_index = bulk_pt->u_index_axisym_poroelasticity(i);
488
489 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
490 interpolated_T1[i] +=
491 (this->nodal_position(l, i) + nodal_value(l, u_nodal_index)) *
492 dpsids(l, 0);
493 }
494 }
495
496 // Set the Jacobian of the undeformed line element
497 double J_undef = sqrt(interpolated_t1[0] * interpolated_t1[0] +
498 interpolated_t1[1] * interpolated_t1[1]) *
499 x[0];
500
501
502 // Set the Jacobian of the deformed line element
503 double J_def = sqrt(interpolated_T1[0] * interpolated_T1[0] +
504 interpolated_T1[1] * interpolated_T1[1]) *
505 (x[0] + disp[0]);
506
507 double return_val = 1.0;
508 if (J_def != 0.0) return_val = J_undef / J_def;
509 return return_val;
510 }
511
512 /// Compute traction vector at specified local coordinate
513 /// Should only be used for post-processing; ignores dependence
514 /// on integration point!
515 void traction(const double& time,
516 const Vector<double>& s,
518
519 /// Compute pressure value at specified local coordinate
520 /// Should only be used for post-processing; ignores dependence
521 /// on integration point!
522 void pressure(const double& time,
523 const Vector<double>& s,
524 double& pressure);
525
526
527 /// Compute contributions to integrated porous flux over boundary:
528 /// q_skeleton = \int \partial u_displ / \partial t \cdot n ds
529 /// q_seepage = \int k q \cdot n ds
530 void contribution_to_total_porous_flux(double& skeleton_flux_contrib,
531 double& seepage_flux_contrib)
532 {
533 // Get pointer to bulk element
534 ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(bulk_element_pt());
535
536 // Get permeability from the bulk poroelasticity element
537 const double permeability = bulk_el_pt->permeability();
538
539 // Find out how many nodes there are
540 const unsigned n_node = this->nnode();
541
542 // Set up memeory for the shape functions
543 Shape psi(n_node);
544 DShape dpsids(n_node, 1);
545
546 // Get the value of Nintpt
547 const unsigned n_intpt = integral_pt()->nweight();
548
549 // Set the Vector to hold local coordinates
550 Vector<double> s(1);
551 Vector<double> s_bulk(2);
552
553 // Loop over the integration points
554 skeleton_flux_contrib = 0.0;
555 seepage_flux_contrib = 0.0;
556 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
557 {
558 // Assign values of s
559 s[0] = integral_pt()->knot(ipt, 0);
560
561 // Get the integral weight
562 double W = this->integral_pt()->weight(ipt);
563
564 // Call the derivatives of the shape function at the knot point
565 this->dshape_local_at_knot(ipt, psi, dpsids);
566
567 // Get position and tangent vector
568 Vector<double> interpolated_t1(2, 0.0);
570 for (unsigned l = 0; l < n_node; l++)
571 {
572 // Loop over directional components
573 for (unsigned i = 0; i < 2; i++)
574 {
575 interpolated_x[i] += this->nodal_position(l, i) * psi(l);
576 interpolated_t1[i] += this->nodal_position(l, i) * dpsids(l, 0);
577 }
578 }
579
580 // Calculate the length of the tangent Vector
581 double tlength = interpolated_t1[0] * interpolated_t1[0] +
582 interpolated_t1[1] * interpolated_t1[1];
583
584 // Set the Jacobian of the line element
585 double J = sqrt(tlength) * interpolated_x[0];
586
587 // Get the outer unit normal
588 Vector<double> interpolated_normal(2);
589 outer_unit_normal(s, interpolated_normal);
590
591 // Get coordinate in bulk element
592 s_bulk = this->local_coordinate_in_bulk(s);
593
594 /// Calculate the FE representation of q
595 Vector<double> q(2);
596 bulk_el_pt->interpolated_q(s_bulk, q);
597
598 // Skeleton velocity from bulk
599 Vector<double> du_dt(2);
600 bulk_el_pt->interpolated_du_dt(s_bulk, du_dt);
601
602#ifdef PARANOID
603 Vector<double> x_bulk(2);
604 x_bulk[0] = bulk_el_pt->interpolated_x(s_bulk, 0);
605 x_bulk[1] = bulk_el_pt->interpolated_x(s_bulk, 1);
606 double error = sqrt(
607 (interpolated_x[0] - x_bulk[0]) * (interpolated_x[0] - x_bulk[0]) +
608 (interpolated_x[1] - x_bulk[1]) * (interpolated_x[1] - x_bulk[1]));
609 double tol = 1.0e-10;
610 if (error > tol)
611 {
612 std::stringstream junk;
613 junk << "Gap between bulk and face element coordinate\n"
614 << "is suspiciously large: " << error
615 << "\nBulk at: " << x_bulk[0] << " " << x_bulk[1] << "\n"
616 << "Face at: " << interpolated_x[0] << " " << interpolated_x[1]
617 << "\n";
619 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
620 }
621#endif
622
623 // Get net flux through boundary
624 double q_flux = 0.0;
625 double dudt_flux = 0.0;
626 for (unsigned i = 0; i < 2; i++)
627 {
628 q_flux += permeability * q[i] * interpolated_normal[i];
629 dudt_flux += du_dt[i] * interpolated_normal[i];
630 }
631
632 // Add
633 seepage_flux_contrib +=
634 2.0 * MathematicalConstants::Pi * q_flux * W * J;
635 skeleton_flux_contrib +=
636 2.0 * MathematicalConstants::Pi * dudt_flux * W * J;
637 }
638 }
639 };
640
641 /// ////////////////////////////////////////////////////////////////////
642 /// ////////////////////////////////////////////////////////////////////
643 /// ////////////////////////////////////////////////////////////////////
644
645 //=====================================================================
646 /// Compute traction vector at specified local coordinate
647 /// Should only be used for post-processing; ignores dependence
648 /// on integration point!
649 //=====================================================================
650 template<class ELEMENT>
652 const double& time, const Vector<double>& s, Vector<double>& traction)
653 {
654 unsigned n_dim = this->nodal_dimension();
655
656 // Position vector
657 Vector<double> x(n_dim);
658 interpolated_x(s, x);
659
660 // Outer unit normal
661 Vector<double> unit_normal(n_dim);
662 outer_unit_normal(s, unit_normal);
663
664 // Dummy
665 unsigned ipt = 0;
666
667 // Pressure value
668 get_traction(time, ipt, x, unit_normal, traction);
669 }
670
671 //=====================================================================
672 /// Compute pressure value at specified local coordinate
673 /// Should only be used for post-processing; ignores dependence
674 /// on integration point!
675 //=====================================================================
676 template<class ELEMENT>
678 const double& time, const Vector<double>& s, double& pressure)
679 {
680 unsigned n_dim = this->nodal_dimension();
681
682 // Position vector
683 Vector<double> x(n_dim);
684 interpolated_x(s, x);
685
686 // Outer unit normal
687 Vector<double> unit_normal(n_dim);
688 outer_unit_normal(s, unit_normal);
689
690 // Dummy
691 unsigned ipt = 0;
692
693 // Pressure value
694 get_pressure(time, ipt, x, unit_normal, pressure);
695 }
696
697
698 //=====================================================================
699 /// Return the residuals for the AxisymmetricPoroelasticityTractionElement
700 /// equations
701 //=====================================================================
702 template<class ELEMENT>
705 Vector<double>& residuals)
706 {
707 // Find out how many nodes there are
708 unsigned n_node = nnode();
709
710 // Get continuous time from timestepper of first node
711 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
712
713#ifdef PARANOID
714 // Find out how many positional dofs there are
715 unsigned n_position_type = this->nnodal_position_type();
716 if (n_position_type != 1)
717 {
718 throw OomphLibError("Poroelasticity equations are not yet implemented "
719 "for more than one position type",
720 OOMPH_CURRENT_FUNCTION,
721 OOMPH_EXCEPTION_LOCATION);
722 }
723#endif
724
725 // Find out the dimension of the node
726 unsigned n_dim = this->nodal_dimension();
727
728
729 // Get bulk element
730 ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(bulk_element_pt());
731
732 unsigned n_q_basis = bulk_el_pt->nq_basis();
733 unsigned n_q_basis_edge = bulk_el_pt->nq_basis_edge();
734
735 // Integer to hold the local equation number
736 int local_eqn = 0;
737
738 // Set up memory for the shape functions
739 // Note that in this case, the number of lagrangian coordinates is always
740 // equal to the dimension of the nodes
741 Shape psi(n_node);
742 DShape dpsids(n_node, n_dim - 1);
743 Shape q_basis(n_q_basis, n_dim);
744
745 // Set the value of n_intpt
746 unsigned n_intpt = integral_pt()->nweight();
747
748 // Storage for the local coordinates
749 Vector<double> s_face(n_dim - 1, 0.0), s_bulk(n_dim, 0.0);
750
751 // Loop over the integration points
752 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
753 {
754 // Get the integral weight
755 double w = integral_pt()->weight(ipt);
756
757 // Only need to call the local derivatives
758 dshape_local_at_knot(ipt, psi, dpsids);
759
760 // Assign values of s in FaceElement and local coordinates in bulk element
761 for (unsigned i = 0; i < n_dim - 1; i++)
762 {
763 s_face[i] = integral_pt()->knot(ipt, i);
764 }
765 s_bulk = local_coordinate_in_bulk(s_face);
766
767 // Get bulk element
768 ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(bulk_element_pt());
769
770 // Get the q basis at bulk local coordinate s_bulk,
771 // corresponding to face local
772 // coordinate s_face
773 bulk_el_pt->get_q_basis(s_bulk, q_basis);
774
775 // Calculate the Eulerian and Lagrangian coordinates
776 Vector<double> interpolated_x(n_dim, 0.0);
777
778 // Also calculate the surface Vectors (derivatives wrt local coordinates)
779 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
780
781 // Calculate displacements and derivatives
782 for (unsigned l = 0; l < n_node; l++)
783 {
784 // Loop over directions
785 for (unsigned i = 0; i < n_dim; i++)
786 {
787 // Calculate the Eulerian coords
788 const double x_local = nodal_position(l, i);
789 interpolated_x[i] += x_local * psi(l);
790
791 // Loop over LOCAL derivative directions, to calculate the tangent(s)
792 for (unsigned j = 0; j < n_dim - 1; j++)
793 {
794 interpolated_A(j, i) += x_local * dpsids(l, j);
795 }
796 }
797 }
798
799 // Now find the local metric tensor from the tangent vectors
800 DenseMatrix<double> A(n_dim - 1);
801 for (unsigned i = 0; i < n_dim - 1; i++)
802 {
803 for (unsigned j = 0; j < n_dim - 1; j++)
804 {
805 // Initialise surface metric tensor to zero
806 A(i, j) = 0.0;
807
808 // Take the dot product
809 for (unsigned k = 0; k < n_dim; k++)
810 {
811 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
812 }
813 }
814 }
815
816 // Get the outer unit normal
817 Vector<double> interpolated_normal(n_dim);
818 outer_unit_normal(ipt, interpolated_normal);
819
820 // Find the determinant of the metric tensor
821 double Adet = 0.0;
822 switch (n_dim)
823 {
824 case 2:
825 Adet = A(0, 0);
826 break;
827 case 3:
828 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
829 break;
830 default:
831 throw OomphLibError(
832 "Wrong dimension in AxisymmetricPoroelasticityTractionElement",
833 OOMPH_CURRENT_FUNCTION,
834 OOMPH_EXCEPTION_LOCATION);
835 }
836
837 // Premultiply the weights and the square-root of the determinant of
838 // the metric tensor
839 double W = w * sqrt(Adet);
840
841 // Now calculate the traction load
842 Vector<double> traction(n_dim);
843 get_traction(time, ipt, interpolated_x, interpolated_normal, traction);
844
845 // Now calculate the load
846 double pressure;
847 get_pressure(time, ipt, interpolated_x, interpolated_normal, pressure);
848
849 // Loop over the test functions, nodes of the element
850 for (unsigned l = 0; l < n_node; l++)
851 {
852 // Loop over the displacement components
853 for (unsigned i = 0; i < n_dim; i++)
854 {
855 local_eqn = this->nodal_local_eqn(
856 l, bulk_el_pt->u_index_axisym_poroelasticity(i));
857 /*IF it's not a boundary condition*/
858 if (local_eqn >= 0)
859 {
860 // Add the traction loading terms to the residuals
861 residuals[local_eqn] -=
862 traction[i] * psi(l) * interpolated_x[0] * W;
863 } // End of if not boundary condition
864 }
865 } // End of loop over shape functions
866
867 // Loop over the q edge test functions only (internal basis functions
868 // have zero normal component on the boundary)
869 for (unsigned l = 0; l < n_q_basis_edge; l++)
870 {
871 local_eqn = this->nodal_local_eqn(1, bulk_el_pt->q_edge_index(l));
872
873 /*IF it's not a boundary condition*/
874 if (local_eqn >= 0)
875 {
876 // Loop over the displacement components
877 for (unsigned i = 0; i < n_dim; i++)
878 {
879 // Add the loading terms to the residuals
880 residuals[local_eqn] += pressure * q_basis(l, i) *
881 interpolated_normal[i] * interpolated_x[0] *
882 W;
883 }
884 } // End of if not boundary condition
885 } // End of loop over shape functions
886 } // End of loop over integration points
887 }
888
889 /// ////////////////////////////////////////////////////////////////////
890 /// ////////////////////////////////////////////////////////////////////
891 /// ////////////////////////////////////////////////////////////////////
892
893 //======================================================================
894 /// A class for elements that allow the imposition of an applied combined
895 /// traction and pore fluid pressure in the poroelasticity equations.
896 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
897 /// class and thus, we can be generic enough without the need to have
898 /// a separate equations class.
899 //======================================================================
900 template<class POROELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
903 POROELASTICITY_BULK_ELEMENT>,
904 public virtual ElementWithExternalElement
905 {
906 protected:
907 /// Pointer to the ratio, \f$ Q \f$, of the stress used to
908 /// non-dimensionalise the fluid stresses to the stress used to
909 /// non-dimensionalise the poroelastic stresses.
910 double* Q_pt;
911
912 /// Static default value for the ratio of stress scales
913 /// used in the fluid and poroelasticity equations (default is 1.0)
914 static double Default_Q_Value;
915
916
917 public:
918 /// Return the ratio of the stress scales used to non-dimensionalise
919 /// the fluid and poroelasticity equations.
920 const double& q() const
921 {
922 return *Q_pt;
923 }
924
925 /// Return a pointer the ratio of stress scales used to
926 /// non-dimensionalise the fluid and poroelastic equations.
927 double*& q_pt()
928 {
929 return Q_pt;
930 }
931
932 /// Get the (combined) traction from the neighbouring Navier-Stokes
933 /// bulk element's stress
934 void get_traction(const double& time,
935 const unsigned& intpt,
936 const Vector<double>& x,
937 const Vector<double>& n,
939 {
940 // Get traction from Navier-Stokes
941 NAVIER_STOKES_BULK_ELEMENT* ext_el_pt =
942 dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(
943 external_element_pt(0, intpt));
945
946 // Find the dimension of the problem
947 unsigned n_dim = this->nodal_dimension();
948
949 // Get the outer unit normal
950 Vector<double> interpolated_normal(n_dim);
951 this->outer_unit_normal(intpt, interpolated_normal);
952
953 // Get all 3 (r,z,theta) components of the nst traction
954 Vector<double> traction_nst(3);
955 ext_el_pt->traction(s_ext, interpolated_normal, traction_nst);
956
957
958#ifdef PARANOID
960 {
961 // Get own coordinates:
962 Vector<double> s(n_dim - 1);
963 for (unsigned i = 0; i < (n_dim - 1); i++)
964 {
965 s[i] = integral_pt()->knot(intpt, i);
966 }
967 Vector<double> x_local(n_dim);
968 this->interpolated_x(s, x_local);
969
970 // Get bulk coordinates in external element
971 Vector<double> x_bulk(n_dim);
972 x_bulk[0] = ext_el_pt->interpolated_x(s_ext, 0);
973 x_bulk[1] = ext_el_pt->interpolated_x(s_ext, 1);
974 double error =
975 sqrt((x_local[0] - x_bulk[0]) * (x_local[0] - x_bulk[0]) +
976 (x_local[1] - x_bulk[1]) * (x_local[1] - x_bulk[1]));
977 double tol = 1.0e-10;
978 if (error > tol)
979 {
980 std::stringstream junk;
981 junk << "Gap between external and face element coordinate\n"
982 << "is suspiciously large:" << error << " ( tol = " << tol
983 << " ) "
984 << "\nExternal/bulk at: " << x_bulk[0] << " " << x_bulk[1]
985 << "\n"
986 << "Face at: " << x_local[0] << " " << x_local[1] << "\n";
987 throw OomphLibError(
988 junk.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
989 }
990 }
991#endif
992
993
994 // Get FSI parameter
995 const double q_value = q();
996
997 // Just assign the r and z components - we require the NST "swirl"
998 // velocity to be 0 since axisym_poroelasticity doesn't support a theta
999 // component. Premulitply by the FSI parameter Q
1000 traction[0] = q_value * traction_nst[0];
1001 traction[1] = q_value * traction_nst[1];
1002 }
1003
1004 /// Get the pore fluid pressure from the neighbouring Navier-Stokes
1005 /// bulk element's stress
1006 void get_pressure(const double& time,
1007 const unsigned& intpt,
1008 const Vector<double>& x,
1009 const Vector<double>& n,
1010 double& pressure)
1011 {
1012 // Get pressure from Navier-Stokes
1013 NAVIER_STOKES_BULK_ELEMENT* ext_el_pt =
1014 dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(
1015 external_element_pt(0, intpt));
1017
1018 // Find the dimension of the problem
1019 unsigned n_dim = this->nodal_dimension();
1020
1021 // Get the outer unit normal
1022 Vector<double> interpolated_normal(n_dim);
1023 this->outer_unit_normal(intpt, interpolated_normal);
1024
1025 // Get all 3 components of the nst traction
1026 Vector<double> traction_nst(3);
1027 ext_el_pt->traction(s_ext, interpolated_normal, traction_nst);
1028
1029 // Get FSI parameter
1030 const double q_value = q();
1031
1032 // Just use the r and z components the calculate the pore pressure - we
1033 // require the "swirl" nst velocity to be 0 since axisym_poroelasticity
1034 // doesn't support a theta component. Premultiply by the FSI parameter Q
1035 pressure = -q_value * (interpolated_normal[0] * traction_nst[0] +
1036 interpolated_normal[1] * traction_nst[1]);
1037 }
1038
1039 public:
1040 /// Constructor, which takes a "bulk" element and the
1041 /// value of the index and its limit
1043 FiniteElement* const& element_pt, const int& face_index)
1044 : AxisymmetricPoroelasticityTractionElement<POROELASTICITY_BULK_ELEMENT>(
1045 element_pt, face_index),
1047 {
1048 // Set source element storage: one interaction with an external
1049 // element -- the Navier Stokes bulk element that provides the traction
1050 this->set_ninteraction(1);
1051 }
1052
1053 /// Default constructor
1055
1056 /// Output function -- overloaded version -- ignores
1057 /// n_plot since fsi elements can only evaluate traction at
1058 /// Gauss points.
1059 void output(std::ostream& outfile, const unsigned& n_plot)
1060 {
1061 // Get continuous time from timestepper of first node
1062 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1063 unsigned n_dim = this->nodal_dimension();
1064
1065 Vector<double> x(n_dim);
1066 Vector<double> s(n_dim - 1);
1067 Vector<double> s_bulk(n_dim);
1068 Vector<double> disp(n_dim);
1069
1070 // Set the value of n_intpt
1071 unsigned n_intpt = integral_pt()->nweight();
1072
1073 // Tecplot header info
1074 outfile << this->tecplot_zone_string(n_intpt);
1075
1076 // Loop over the integration points
1077 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1078 {
1079 // Assign values of s in FaceElement and local coordinates in bulk
1080 // element
1081 for (unsigned i = 0; i < n_dim - 1; i++)
1082 {
1083 s[i] = integral_pt()->knot(ipt, i);
1084 }
1085
1086 // Get Eulerian coordinates
1087 this->interpolated_x(s, x);
1088
1089 // Outer unit normal
1090 Vector<double> unit_normal(n_dim);
1091 this->outer_unit_normal(s, unit_normal);
1092
1093 // Get pointer to bulk element
1094 POROELASTICITY_BULK_ELEMENT* bulk_pt =
1095 dynamic_cast<POROELASTICITY_BULK_ELEMENT*>(this->bulk_element_pt());
1096 s_bulk = this->local_coordinate_in_bulk(s);
1097
1098 // Get permeability from the bulk poroelasticity element
1099 const double local_permeability = bulk_pt->permeability();
1100
1101 // Porous seepage flux
1102 Vector<double> q(2);
1103 bulk_pt->interpolated_q(s_bulk, q);
1104
1105 /// Calculate the FE representation of u
1106 bulk_pt->interpolated_u(s_bulk, disp);
1107
1108 // Now calculate the traction load
1109 Vector<double> traction(n_dim);
1110 this->get_traction(time, ipt, x, unit_normal, traction);
1111
1112 // Now calculate the load
1113 double pressure;
1114 this->get_pressure(time, ipt, x, unit_normal, pressure);
1115
1116
1117 // Output the x,y,..
1118 for (unsigned i = 0; i < n_dim; i++)
1119 {
1120 outfile << x[i] << " ";
1121 } // column 1,2
1122
1123 // Output displacement
1124 for (unsigned i = 0; i < n_dim; i++)
1125 {
1126 outfile << disp[i] << " "; // column 3,4
1127 }
1128
1129 // Output traction
1130 for (unsigned i = 0; i < n_dim; i++)
1131 {
1132 outfile << traction[i] << " "; // column 5,6
1133 }
1134
1135 // Output pressure
1136 outfile << pressure << " "; // column 7
1137
1138 // Output seepage flux
1139 outfile << local_permeability * q[0] << " " // column 8
1140 << local_permeability * q[1] << " "; // column 9
1141
1142 // Output skeleton velocity
1143 Vector<double> du_dt(2);
1144 bulk_pt->interpolated_du_dt(s_bulk, du_dt);
1145 outfile << du_dt[0] << " " // column 10
1146 << du_dt[1] << " "; // column 11
1147
1148 // Total veloc
1149 outfile << du_dt[0] + local_permeability * q[0] << " " // column 12
1150 << du_dt[1] + local_permeability * q[1] << " "; // column 13
1151
1152 // Outer unit normal
1153 outfile << unit_normal[0] << " " // column 14
1154 << unit_normal[1] << " "; // column 15
1155
1156 // Output FE representation of div u at s_bulk
1157 outfile << bulk_pt->interpolated_div_q(s_bulk) << " "; // column 16
1158
1159 // Output FE representation of p at s_bulk
1160 outfile << bulk_pt->interpolated_p(s_bulk) << " "; // column 17
1161 outfile << std::endl;
1162 }
1163
1164 // Write tecplot footer (e.g. FE connectivity lists)
1165 this->write_tecplot_zone_footer(outfile, n_plot);
1166 }
1167
1168
1169 /// Fill in contribution from Jacobian
1171 DenseMatrix<double>& jacobian)
1172 {
1173 // Call the residuals
1175 residuals);
1176
1177 // Derivatives w.r.t. external data
1179 }
1180 };
1181
1182 //=================================================================
1183 /// Static default value for the ratio of stress scales
1184 /// used in the fluid and solid equations (default is 1.0)
1185 //=================================================================
1186 template<class POROELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
1187 double FSILinearisedAxisymPoroelasticTractionElement<
1188 POROELASTICITY_BULK_ELEMENT,
1189 NAVIER_STOKES_BULK_ELEMENT>::Default_Q_Value = 1.0;
1190
1191} // namespace oomph
1192#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class for elements that allow the imposition of an applied combined traction and pore fluid pressur...
virtual void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the traction vector: Pass number of integration point (dummy), Eulerrian coordinate and normal ve...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
virtual void get_pressure(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, double &pressure)
Get the pressure value: Pass number of integration point (dummy), Eulerrian coordinate and normal vec...
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
void contribution_to_total_porous_flux(double &skeleton_flux_contrib, double &seepage_flux_contrib)
Compute contributions to integrated porous flux over boundary: q_skeleton = \int \partial u_displ / \...
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
double lagrangian_eulerian_translation_factor(const Vector< double > &s)
Ratio of lengths of line elements (or annular surface areas) in the undeformed and deformed configura...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
void(* Pressure_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, double &result)
Pointer to an imposed pressure function. Arguments: Eulerian coordinate; outer unit normal; applied p...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void fill_in_contribution_to_residuals_axisym_poroelasticity_face(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void pressure(const double &time, const Vector< double > &s, double &pressure)
Compute pressure value at specified local coordinate Should only be used for post-processing; ignores...
AxisymmetricPoroelasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void traction(const double &time, const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, double &pressure) pressure_fct_pt()
Reference to the pressure function pointer.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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...
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void get_pressure(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, double &pressure)
Get the pore fluid pressure from the neighbouring Navier-Stokes bulk element's stress.
FSILinearisedAxisymPoroelasticTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress use...
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and poroelastic equa...
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and poroelasticity equatio...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and poroelasticity equations (d...
void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the (combined) traction from the neighbouring Navier-Stokes bulk element's stress.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – overloaded version – ignores n_plot since fsi elements can only evaluate traction a...
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
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Definition: elements.h:4497
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
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 double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
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 dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
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
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
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.
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....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
bool Allow_gap_in_FSI
Public boolean to allow gap between poro-elastic and Navier Stokes element in FSI computations....
void Zero_pressure_fct(const double &time, const Vector< double > &x, const Vector< double > &N, double &load)
Default load function (zero pressure)
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void output()
Doc the command line arguments.
const double Pi
50 digits from maple
//////////////////////////////////////////////////////////////////// ////////////////////////////////...