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