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