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