specific_node_update_interface_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 specific (two-dimensional) fluid free surface elements
31 
32 //Include guards, to prevent multiple includes
33 #ifndef OOMPH_SPECIFIC_NODE_UPDATE_INTERFACE_ELEMENTS_HEADER
34 #define OOMPH_SPECIFIC_NODE_UPDATE_INTERFACE_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <oomph-lib-config.h>
39 #endif
40 
41 //OOMPH-LIB headers
42 #include "../generic/Qelements.h"
43 #include "../generic/spines.h"
44 #include "../generic/hijacked_elements.h"
45 #include "interface_elements.h"
46 
47 namespace oomph
48 {
49 
50 //=======================================================================
51 /// \short This policy class is used to associate specific bounding
52 /// elements with specific FluidInterface elements. It must be filled
53 /// in for every class that uses the SpineUpdateFluidInterface<...>
54 /// or ElasticUpdateFluidInterface<....> generic template classes.
55 /// Examples for our default Line, Axisymmetric and Surface types
56 /// are included below
57 //=======================================================================
58  template<class ELEMENT>
60  {};
61 
62 //======================================================================
63 /// \short This policy class is used to allow additional values to be
64 /// added to the nodes from new surface equations, for examples of
65 /// usage see the SurfactantTransportFluidInterfaceElements.
66 /// The use of this class avoids issues with calling virtual
67 /// functions in constructors and avoids having a global look-up
68 /// able, although it functions in much the same way.
69 /// Typically, this will only be filled in by "expert users" and
70 /// is only required if you want to write generic surface-element
71 /// classes. Specific classes can always be overloaded on a
72 /// case-by-case basis.
73 //=====================================================================
74  template<class ELEMENT>
76  {
77  private:
78  //Issue an error message if this base class is called
80  {
81  std::ostringstream error_message;
82  error_message
83  <<
84  "The generic FluidInterfaceAdditionalValues<ELEMENT> class has been\n"
85  <<
86  "called. This should never happen. If you are creating your own\n"
87  <<
88  "surface equations you must overload the policy class to specify\n"
89  <<
90  "how many additional values are required at the surface nodes.\n"
91  <<
92  "For an example, see src/fluid_interface/surfactant_transport_elements.h\n"
93  << std::endl;
94  throw OomphLibError(error_message.str(),
95  OOMPH_CURRENT_FUNCTION,
96  OOMPH_EXCEPTION_LOCATION);
97  }
98 
99  public:
100 
101  ///Empty constructor
103 
104  ///\short Specific interface that states how many additional values are
105  ///required for the n-th node. Default is zero, but issue error_message.
106  inline unsigned nadditional_values(const unsigned &n)
107  {error_message(); return 0;}
108 
109  ///\short Specify any additional index setup information that is required;
110  ///i.e. the look-up schemes for the additional values.
111  ///Default is empty with error message
112  inline void setup_equation_indices(ELEMENT* const &element_pt, const unsigned &id)
113  {error_message();}
114  };
115 
116 
117 //======================================================================
118 /// \short Specific policy class for the FluidInterfaceElemetnts,
119 /// which do not require any additional values at the nodes.
120  //=====================================================================
121  template<>
123  {
124  public:
125 
126  ///Empty constructor
128 
129  ///\short Specific interface that states how many additional values are
130  ///required for the n-th node. No additional values
131  inline unsigned nadditional_values(const unsigned &n)
132  {return 0;}
133 
134  ///\short Specify any additional index setup information that is required;
135  ///i.e. the look-up schemes for the additional values.
136  ///Empty
137  inline void setup_equation_indices(FluidInterfaceElement* const &element_pt, const unsigned &id)
138  {}
139  };
140 
141 
142 
143 //-------------SPINE NODE UPDATE CLASSES-------------------------------
144 //---------------------------------------------------------------------
145 
146 //======================================================================
147 /// \short Generic Spine node update interface template class that can
148 /// be combined with a given surface equations class and surface
149 /// derivative class to provide
150 /// a concrete implementation of any surface element that uses spines.
151 //======================================================================
152  template<class EQUATION_CLASS, class DERIVATIVE_CLASS, class ELEMENT>
154  public virtual Hijacked<SpineElement<FaceGeometry<ELEMENT> > >,
155  public virtual EQUATION_CLASS, public DERIVATIVE_CLASS
156  {
157  private:
158 
159  /// \short In spine elements, the kinematic condition is the equation
160  /// used to determine the unknown spine heights. Overload the
161  /// function accordingly
162  int kinematic_local_eqn(const unsigned &n)
163  {return this->spine_local_eqn(n);}
164 
165  /// \short Hijacking the kinematic condition corresponds to hijacking the
166  /// variables associated with the spine heights.
167  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
168  {
169  //Loop over all the node numbers that are passed in
170  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
171  it!=bulk_node_number.end();++it)
172  {
173  //Hijack the spine heights. (and delete the returned data object)
174  delete this->hijack_nodal_spine_value(*it,0);
175  }
176  }
177 
178  protected:
179 
180  ///\short Fill in the specific surface derivative calculations
181  /// by calling the appropriate class function
182  double compute_surface_derivatives(const Shape &psi, const DShape &dpsids,
183  const DenseMatrix<double> &interpolated_t,
184  const Vector<double> &interpolated_x,
185  DShape &surface_gradient,
186  DShape &surface_divergence)
187  {return DERIVATIVE_CLASS::compute_surface_derivatives(
188  psi,dpsids,interpolated_t,interpolated_x,surface_gradient,
189  surface_divergence);}
190 
191 
192  public:
193 
194 
195  /// \short Constructor, the arguments are a pointer to the "bulk" element
196  /// and the index of the face to be created
197  SpineUpdateFluidInterfaceElement(FiniteElement* const &element_pt,
198  const int &face_index,
199  const unsigned &id=0) :
200  Hijacked<SpineElement<FaceGeometry<ELEMENT> > >(),
201  EQUATION_CLASS(), DERIVATIVE_CLASS()
202  {
203  //Attach the geometrical information to the element, by
204  //making the face element from the bulk element
205  element_pt->build_face_element(face_index,this);
206 
207 #ifdef PARANOID
208  //Is it refineable
209  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(element_pt);
210  if(ref_el_pt!=0)
211  {
212  if (this->has_hanging_nodes())
213  {
214  throw OomphLibError(
215  "This interface element will not work correctly if nodes are hanging\n",
216  OOMPH_CURRENT_FUNCTION,
217  OOMPH_EXCEPTION_LOCATION);
218  }
219  }
220 #endif
221 
222  //Find the index at which the velocity unknowns are stored
223  //from the bulk element and allocate the local storage
224  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
225  //Find number of momentum equation required
226  const unsigned n_u_index = cast_element_pt->n_u_nst();
227  this->U_index_interface.resize(n_u_index);
228  for(unsigned i=0;i<n_u_index;i++)
229  {
230  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
231  }
232 
233  //Add any additional values required by the equations class
234  unsigned n_node_face = this->nnode();
235  //Create an instance of the policy class that determines
236  //how many additional values are required
238  *interface_additional_values_pt =
240  //Do we need to add any values
241  //By default the spines don't need to so we can save
242  //some effort here
243  bool add_values = false;
244  //Storage for the number of additional values required
245  //for each node
246  Vector<unsigned> additional_data_values(n_node_face);
247  for(unsigned i=0;i<n_node_face;++i)
248  {
249  //Read out additional values for each node
250  additional_data_values[i] = interface_additional_values_pt->nadditional_values(i);
251  //If any of these are non-zero it will set the boolean to true
252  add_values |= additional_data_values[i];
253  }
254 
255  // If storage is needed allocate it and call
256  // the associated local index setup function
257  if(add_values)
258  {
259  this->add_additional_values(additional_data_values,id);
260  //Call any local setup from the interface policy class
261  interface_additional_values_pt->setup_equation_indices(this,id);
262  }
263 
264  //Delete the policy class, it's work is done.
265  delete interface_additional_values_pt; interface_additional_values_pt=0;
266  }
267 
268  /// Calculate the contribution to the residuals and the jacobian
269  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
270  DenseMatrix<double> &jacobian)
271  {
272  //Call the generic routine with the flag set to 1
273  EQUATION_CLASS::fill_in_generic_residual_contribution_interface(residuals,jacobian,1);
274 
275  //Call the generic routine to handle the shape derivatives
276  this->fill_in_jacobian_from_geometric_data(jacobian);
277  }
278 
279  /// \short
280  /// Helper function to calculate the additional contributions
281  /// These are those filled in by the particular equations
283  Vector<double> &residuals,
284  DenseMatrix<double> &jacobian,
285  const unsigned &flag,
286  const Shape &psif,
287  const DShape &dpsifds,
288  const DShape &dpsifdS,
289  const DShape &dpsifdS_div,
290  const Vector<double> &s,
291  const Vector<double> &interpolated_x,
292  const Vector<double> &interpolated_n,
293  const double &W,
294  const double &J)
295  {
296  EQUATION_CLASS::add_additional_residual_contributions_interface(
297  residuals,jacobian,flag,psif,dpsifds,dpsifdS,dpsifdS_div,s,interpolated_x,interpolated_n,W,J);
298  }
299 
300  /// Overload the output function
301  void output(std::ostream &outfile) {EQUATION_CLASS::output(outfile);}
302 
303  /// Output the element
304  void output(std::ostream &outfile, const unsigned &n_plot)
305  {EQUATION_CLASS::output(outfile,n_plot);}
306 
307  ///Overload the C-style output function
308  void output(FILE* file_pt) {EQUATION_CLASS::output(file_pt);}
309 
310  ///C-style Output function
311  void output(FILE* file_pt, const unsigned &n_plot)
312  {EQUATION_CLASS::output(file_pt,n_plot);}
313 
314 
315  /// \short Create an "bounding" element of the type
316  /// specified by the BoundingElementType policy class
317  /// Here, this allows
318  /// the application of a contact angle boundary condition on the
319  /// the specified face.
321  const int &face_index)
322  {
323  //Create a temporary pointer to the appropriate FaceElement read our from
324  //our policy class
326  DERIVATIVE_CLASS,ELEMENT> > *face_el_pt =
328  DERIVATIVE_CLASS,ELEMENT> >;
329 
330  //Attach the geometrical information to the new element
331  this->build_face_element(face_index,face_el_pt);
332 
333  //Set the index at which the velocity nodes are stored
334  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
335 
336  //Set the value of the nbulk_value, the node is not resized
337  //in this bounding element,
338  //so it will just be the actual nvalue here
339  // There is some ambiguity about what this means (however)
340  // We are interpreting it to mean the number of
341  // values in this FaceElement before creating the new
342  // bounding element.
343  const unsigned n_node_bounding = face_el_pt->nnode();
344  for(unsigned n=0;n<n_node_bounding;n++)
345  {
346  face_el_pt->nbulk_value(n) =
347  face_el_pt->node_pt(n)->nvalue();
348  }
349 
350  //Set of unique geometric data that is used to update the bulk,
351  //but is not used to update the face
352  std::set<Data*> unique_additional_geom_data;
353 
354  //Get all the geometric data for this (bulk) element
355  this->assemble_set_of_all_geometric_data(unique_additional_geom_data);
356 
357  //Now assemble the set of geometric data for the face element
358  std::set<Data*> unique_face_geom_data_pt;
359  face_el_pt->assemble_set_of_all_geometric_data(unique_face_geom_data_pt);
360 
361  //Erase the face geometric data from the additional data
362  for(std::set<Data*>::iterator it=unique_face_geom_data_pt.begin();
363  it!=unique_face_geom_data_pt.end();++it)
364  {unique_additional_geom_data.erase(*it);}
365 
366  //Finally add all unique additional data as external data
367  for(std::set<Data*>::iterator it = unique_additional_geom_data.begin();
368  it!= unique_additional_geom_data.end();++it)
369  {
370  face_el_pt->add_external_data(*it);
371  }
372 
373  //Return the value of the pointer
374  return face_el_pt;
375  }
376 
377  };
378 
379 
380 //=====================================================================
381 ///Spine version of the PointFluidInterfaceBoundingElement
382 //=====================================================================
383  template<class ELEMENT>
385  public Hijacked< SpineElement<FaceGeometry<FaceGeometry<ELEMENT> > > >,
387 
388  {
389  public:
390 
391 
392  /// Constructor
394  Hijacked<SpineElement<FaceGeometry<FaceGeometry<ELEMENT> > > >(),
396 
397  /// Overload the output function
398  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
399 
400  /// Output the element
401  void output(std::ostream &outfile, const unsigned &n_plot)
402  {FluidInterfaceBoundingElement::output(outfile,n_plot);}
403 
404  ///Overload the C-style output function
405  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
406 
407  ///C-style Output function
408  void output(FILE* file_pt, const unsigned &n_plot)
409  {FluidInterfaceBoundingElement::output(file_pt,n_plot);}
410 
411  /// Calculate the elemental residual vector and the Jacobian
412  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
413  DenseMatrix<double> &jacobian)
414  {
415  //Call the generic routine with the flag set to 1
417  jacobian,1);
418  //Call generic FD routine for the external data
419  this->fill_in_jacobian_from_external_by_fd(jacobian);
420 
421  //Call the generic routine to handle the spine variables
422  this->fill_in_jacobian_from_geometric_data(jacobian);
423  }
424 
425  /// \short Return local equation number associated with the kinematic
426  /// constraint for local node n
427  int kinematic_local_eqn(const unsigned &n)
428  {return this->spine_local_eqn(n);}
429 
430  };
431 
432 
433 
434 //=========================================================================
435 /// Spine version of the LineFluidInterfaceBoundingElement
436 //========================================================================
437  template<class ELEMENT>
439  Hijacked<SpineElement<FaceGeometry<FaceGeometry<ELEMENT> > > >,
441 
442  {
443  public:
444 
445  /// Constructor
447  Hijacked<SpineElement<FaceGeometry<FaceGeometry<ELEMENT> > > >(),
449 
450  /// Overload the output function
451  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
452 
453  /// Output the element
454  void output(std::ostream &outfile, const unsigned &n_plot)
455  {FluidInterfaceBoundingElement::output(outfile,n_plot);}
456 
457  ///Overload the C-style output function
458  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
459 
460  ///C-style Output function
461  void output(FILE* file_pt, const unsigned &n_plot)
462  {FluidInterfaceBoundingElement::output(file_pt,n_plot);}
463 
464 
465  /// Calculate the jacobian
466  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
467  DenseMatrix<double> &jacobian)
468  {
469  //Call the generic routine with the flag set to 1
471  jacobian,1);
472  //Call generic FD routine for the external data
473  this->fill_in_jacobian_from_external_by_fd(jacobian);
474 
475  //Call the generic routine to handle the spine variables
476  this->fill_in_jacobian_from_geometric_data(jacobian);
477  }
478 
479 
480  /// Local eqn number of the kinematic bc associated with local node n
481  int kinematic_local_eqn(const unsigned &n)
482  {
483  //Kinematic bc is always associated with the n-th spine height
484  return this->spine_local_eqn(n);
485  }
486 
487  };
488 
489 
490 //============GEOMETRIC SPECIALISATIONS============================
491 
492 
493 //Specialise the spine update template class to concrete 1D case
494  template<class ELEMENT>
495  class SpineLineFluidInterfaceElement :
496  public SpineUpdateFluidInterfaceElement<FluidInterfaceElement,
497  LineDerivatives,ELEMENT>
498  {
499  public:
500 
502  FiniteElement* const &element_pt,
503  const int &face_index) :
505  LineDerivatives,ELEMENT>(element_pt,face_index) {}
506  };
507 
508 
509 //Define the BoundingElement type associated with the 1D surface element
510  template<class ELEMENT>
514  {
515  public:
516 
520  };
521 
522 
523 //Specialise Spine update case to concrete axisymmetric case
524  template<class ELEMENT>
525  class SpineAxisymmetricFluidInterfaceElement :
526  public SpineUpdateFluidInterfaceElement<FluidInterfaceElement,
527  AxisymmetricDerivatives,ELEMENT>
528  {
529  public:
530 
532  FiniteElement* const &element_pt,
533  const int &face_index) :
535  AxisymmetricDerivatives,ELEMENT>(element_pt,face_index) {}
536  };
537 
538 //Define the bounding element associated with the axisymmetric fluid interface element
539  template<class ELEMENT>
543  {
544  public:
545 
549  };
550 
551 //Specialise Spine update case to concrete 2D case
552  template<class ELEMENT>
554  public SpineUpdateFluidInterfaceElement<FluidInterfaceElement,
555  SurfaceDerivatives,ELEMENT>
556  {
557  public:
558 
560  FiniteElement* const &element_pt,
561  const int &face_index) :
563  SurfaceDerivatives,ELEMENT>(element_pt,face_index) {}
564  };
565 
566 //Define the bounding element type for the 2D surface
567  template<class ELEMENT>
571  {
572  public:
573 
577  };
578 
579 
580 ///////////////////////////////////////////////////////////////////////
581 ///////////////////////////////////////////////////////////////////////
582 ///////////////////////////////////////////////////////////////////////
583 
584 //-------------ELASTIC NODE UPDATE CLASSES-------------------------------
585 //---------------------------------------------------------------------
586 
587 
588 //=======================================================================
589 /// \short Generic Elastic node update interface template class that can
590 /// be combined with a given surface equations class and surface derivative
591 /// class to provide a concrete implementation of any surface element
592 /// that uses elastic node updates.
593 //======================================================================
594  template<class EQUATION_CLASS, class DERIVATIVE_CLASS, class ELEMENT>
596  public virtual Hijacked<FaceGeometry<ELEMENT> >,
597  public EQUATION_CLASS, public DERIVATIVE_CLASS
598  {
599  private:
600 
601  /// \short Storage for the location of the Lagrange multiplier
602  /// (If other additional values have been added we need
603  /// to add the Lagrange multiplier at the end)
604  Vector<unsigned> Lagrange_index;
605 
606  /// \short Return the index at which the lagrange multiplier is
607  ///stored at the n-th node
608  inline unsigned lagrange_index(const unsigned &n)
609  {return this->Lagrange_index[n];}
610 
611  /// \short Equation number of the kinematic BC associated with node j.
612  /// (This is the equation for the Lagrange multiplier)
613  inline int kinematic_local_eqn(const unsigned &n)
614  {
615  // Get the index of the nodal value associated with Lagrange multiplier
616  return this->nodal_local_eqn(n,this->lagrange_index(n));
617  }
618 
619  /// \short Hijacking the kinematic condition corresponds to hijacking the
620  /// variables associated with the Lagrange multipliers that are assigned
621  /// on construction of this element.
622  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
623  {
624  //Loop over all the nodes that are passed in
625  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
626  it!=bulk_node_number.end();++it)
627  {
628  //Hijack the appropriate value and delete the returned Node
629  delete this->hijack_nodal_value(*it,this->lagrange_index(*it));
630  }
631  }
632 
633  protected:
634 
635  ///\short Fill in the specific surface derivative calculations
636  /// by calling the appropriate function from the derivative class
637  double compute_surface_derivatives(const Shape &psi, const DShape &dpsids,
638  const DenseMatrix<double> &interpolated_t,
639  const Vector<double> &interpolated_x,
640  DShape &surface_gradient,
641  DShape &surface_divergence)
642  {return DERIVATIVE_CLASS::compute_surface_derivatives(
643  psi,dpsids,interpolated_t,interpolated_x,surface_gradient,
644  surface_divergence);}
645 
646 
647  public:
648 
649  /// \short Constructor, pass a pointer to the bulk element and the face
650  /// index of the bulk element to which the element is to be attached to.
651  /// The optional identifier can be used
652  /// to distinguish the additional nodal value (Lagr mult) created by
653  /// this element from those created by other FaceElements.
654  ElasticUpdateFluidInterfaceElement(FiniteElement* const &element_pt,
655  const int &face_index,
656  const unsigned &id=0) :
657  FaceGeometry<ELEMENT>(), EQUATION_CLASS(), DERIVATIVE_CLASS()
658  {
659  //Attach the geometrical information to the element
660  //This function also assigned nbulk_value from required_nvalue of the
661  //bulk element
662  element_pt->build_face_element(face_index,this);
663 
664 #ifdef PARANOID
665  //Is it refineable
666  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(element_pt);
667  if(ref_el_pt!=0)
668  {
669  if (this->has_hanging_nodes())
670  {
671  throw OomphLibError(
672  "This flux element will not work correctly if nodes are hanging\n",
673  OOMPH_CURRENT_FUNCTION,
674  OOMPH_EXCEPTION_LOCATION);
675  }
676  }
677 #endif
678 
679  //Find the index at which the velocity unknowns are stored
680  //from the bulk element and resize the local storage scheme
681  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
682  const unsigned n_u_index = cast_element_pt->n_u_nst();
683  this->U_index_interface.resize(n_u_index);
684  for(unsigned i=0;i<n_u_index;i++)
685  {
686  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
687  }
688 
689  //Read out the number of nodes on the face
690  unsigned n_node_face = this->nnode();
691 
692  //Create an instance of the policy class that determines
693  //how many additional values are required
695  *interface_additional_values_pt =
697 
698  //Set the additional data values in the face
699  //There is always also one additional values at each node --- the Lagrange multiplier
700  Vector<unsigned> additional_data_values(n_node_face);
701  for(unsigned n=0;n<n_node_face;n++)
702  {
703  //Now add one to the addtional values at every single node
704  additional_data_values[n] = interface_additional_values_pt->nadditional_values(n) + 1;
705  }
706 
707  // Now add storage for Lagrange multipliers and set the map containing
708  // the position of the first entry of this face element's
709  // additional values.
710  this->add_additional_values(additional_data_values,id);
711 
712  //Now I can just store the lagrange index offset to give the storage
713  //location of the nodes
714  Lagrange_index.resize(n_node_face);
715  for(unsigned n=0;n<n_node_face;++n)
716  {
717  Lagrange_index[n] = additional_data_values[n] -1 +
718  dynamic_cast<BoundaryNodeBase*>(this->node_pt(n))->index_of_first_value_assigned_by_face_element(id);
719  }
720 
721  //Call any local setup from the interface policy class
722  interface_additional_values_pt->setup_equation_indices(this,id);
723 
724  //Can now delete the policy class
725  delete interface_additional_values_pt; interface_additional_values_pt=0;
726  }
727 
728  /// \short The "global" intrinsic coordinate of the element when
729  /// viewed as part of a geometric object should be given by
730  /// the FaceElement representation, by default
731  double zeta_nodal(const unsigned &n, const unsigned &k,
732  const unsigned &i) const
733  {return FaceElement::zeta_nodal(n,k,i);}
734 
735  /// Return the lagrange multiplier at local node n
736  double &lagrange(const unsigned &n)
737  {
738  return *this->node_pt(n)->value_pt(this->lagrange_index(n));
739  }
740 
741 
742  /// Fill in contribution to residuals and Jacobian
743  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
744  DenseMatrix<double> &jacobian)
745  {
746  //Call the generic routine with the flag set to 1
747  EQUATION_CLASS::fill_in_generic_residual_contribution_interface(residuals,jacobian,1);
748 
749  //Call the generic finite difference routine for the solid variables
750  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
751  }
752 
753  /// Overload the output function
754  void output(std::ostream &outfile) {EQUATION_CLASS::output(outfile);}
755 
756  /// Output the element
757  void output(std::ostream &outfile, const unsigned &n_plot)
758  {EQUATION_CLASS::output(outfile,n_plot);}
759 
760  ///Overload the C-style output function
761  void output(FILE* file_pt) {EQUATION_CLASS::output(file_pt);}
762 
763  ///C-style Output function
764  void output(FILE* file_pt, const unsigned &n_plot)
765  {EQUATION_CLASS::output(file_pt,n_plot);}
766 
767 
768  /// \short Helper function to calculate the additional contributions
769  /// to be added at each integration point. This deals with
770  /// Lagrange multiplier contribution, as well as any
771  /// additional contributions by the other equations.
773  Vector<double> &residuals,
774  DenseMatrix<double> &jacobian,
775  const unsigned &flag,
776  const Shape &psif,
777  const DShape &dpsifds,
778  const DShape &dpsifdS,
779  const DShape &dpsifdS_div,
780  const Vector<double> &s,
781  const Vector<double> &interpolated_x,
782  const Vector<double> &interpolated_n,
783  const double &W,
784  const double &J)
785  {
786  EQUATION_CLASS::add_additional_residual_contributions_interface(
787  residuals,jacobian,flag,psif,dpsifds,dpsifdS,dpsifdS_div,s,interpolated_x,interpolated_n,W,J);
788 
789  //Assemble Lagrange multiplier by loop over the shape functions
790  const unsigned n_node = this->nnode();
791  //Read out the dimension of the node
792  const unsigned nodal_dimension = this->nodal_dimension();
793 
794  double interpolated_lagrange = 0.0;
795  for(unsigned l=0;l<n_node;l++)
796  {
797  //Note same shape functions used for lagrange multiplier field
798  interpolated_lagrange += lagrange(l)*psif(l);
799  }
800 
801  int local_eqn=0, local_unknown = 0;
802 
803  //Loop over the shape functions
804  for(unsigned l=0;l<n_node;l++)
805  {
806  //Loop over the directions
807  for(unsigned i=0;i<nodal_dimension;i++)
808  {
809  //Now using the same shape functions for the elastic equations,
810  //so we can stay in the loop
811  local_eqn = this->position_local_eqn(l,0,i);
812  if(local_eqn >= 0)
813  {
814  //Add in the Lagrange multiplier contribution
815  residuals[local_eqn] -=
816  interpolated_lagrange*interpolated_n[i]*psif(l)*J*W;
817 
818  //Do the Jacobian calculation
819  if(flag)
820  {
821  //Loop over the nodes
822  for(unsigned l2=0;l2<n_node;l2++)
823  {
824  // Dependence on solid positions will be handled by FDs
825  //That leaves the Lagrange multiplier contribution
826  local_unknown = this->kinematic_local_eqn(l2);
827  if(local_unknown >= 0)
828  {
829  jacobian(local_eqn,local_unknown) -=
830  psif(l2)*interpolated_n[i]*psif(l)*J*W;
831  }
832  }
833  } //End of Jacobian calculation
834  }
835  }
836 
837  }
838  }
839 
840 
841 
842  /// \short Create an "bounding" element (here actually a 2D line element
843  /// of type ElasticLineFluidInterfaceBoundingElement<ELEMENT> that allows
844  /// the application of a contact angle boundary condition on the
845  /// the specified face.
847  const int &face_index)
848  {
849  //Create a temporary pointer to the appropriate FaceElement
851  DERIVATIVE_CLASS,ELEMENT> > *face_el_pt =
853 
854  //Attach the geometrical information to the new element
855  this->build_face_element(face_index,face_el_pt);
856 
857  //Set the index at which the velocity nodes are stored
858  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
859 
860  //Set the value of the nbulk_value, the node is not resized
861  //in this bounding element,
862  //so it will just be the actual nvalue here
863  // There is some ambiguity about what this means (however)
864  const unsigned n_node_bounding = face_el_pt->nnode();
865  //Storage for lagrange multiplier index at the face nodes
866  Vector<unsigned> local_lagrange_index(n_node_bounding);
867  for(unsigned n=0;n<n_node_bounding;n++)
868  {
869  face_el_pt->nbulk_value(n) =
870  face_el_pt->node_pt(n)->nvalue();
871  //At the moment the assumption is that it is stored at all nodes, but that is consistent with
872  //the assumption in this element
873  local_lagrange_index[n] = this->Lagrange_index[face_el_pt->bulk_node_number(n)];
874  }
875 
876  //Pass the ID and offset down
877  face_el_pt->set_lagrange_index(local_lagrange_index);
878 
879  //Find the nodes
880  std::set<SolidNode*> set_of_solid_nodes;
881  const unsigned n_node = this->nnode();
882  for(unsigned n=0;n<n_node;n++)
883  {
884  set_of_solid_nodes.insert(static_cast<SolidNode*>(this->node_pt(n)));
885  }
886 
887  //Delete the nodes from the face
888  //n_node = face_el_pt->nnode();
889  for(unsigned n=0;n<n_node_bounding;n++)
890  {
891  //Set the value of the nbulk_value, from the present element
892  face_el_pt->nbulk_value(n) =
893  this->nbulk_value(face_el_pt->bulk_node_number(n));
894 
895  //Now delete the nodes from the set
896  set_of_solid_nodes.erase(static_cast<SolidNode*>(
897  face_el_pt->node_pt(n)));
898  }
899 
900  //Now add these as external data
901  for(std::set<SolidNode*>::iterator it=set_of_solid_nodes.begin();
902  it!=set_of_solid_nodes.end();++it)
903  {
904  face_el_pt->add_external_data((*it)->variable_position_pt());
905  }
906 
907 
908  //Return the value of the pointer
909  return face_el_pt;
910  }
911 
912  };
913 
914 
915 //=========================================================================
916 /// Pseudo-elasticity version of the PointFluidInterfaceBoundingElement
917 //========================================================================
918  template<class ELEMENT>
920  public FaceGeometry<FaceGeometry<ELEMENT> > ,
921  public PointFluidInterfaceBoundingElement, public virtual SolidFiniteElement
922 
923  {
924  private:
925 
926  /// Short Storage for the index of the Lagrange multiplier at the chosen nodes
927  Vector<unsigned> Lagrange_index;
928 
929  public:
930 
931  ///\short Set the Id and offset
932  void set_lagrange_index(const Vector<unsigned> &lagrange_index)
933  {Lagrange_index = lagrange_index;}
934 
935 
936  /// \short Specify the value of nodal zeta from the face geometry
937  /// The "global" intrinsic coordinate of the element when
938  /// viewed as part of a geometric object should be given by
939  /// the FaceElement representation, by default
940  double zeta_nodal(const unsigned &n, const unsigned &k,
941  const unsigned &i) const
942  {return FaceElement::zeta_nodal(n,k,i);}
943 
944 
945  /// Constructor
947  FaceGeometry<FaceGeometry<ELEMENT> >(),
949 
950  /// Overload the output function
951  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
952 
953  /// Output the element
954  void output(std::ostream &outfile, const unsigned &n_plot)
955  {FluidInterfaceBoundingElement::output(outfile,n_plot);}
956 
957  ///Overload the C-style output function
958  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
959 
960  ///C-style Output function
961  void output(FILE* file_pt, const unsigned &n_plot)
962  {FluidInterfaceBoundingElement::output(file_pt,n_plot);}
963 
964  /// Calculate the element's residual vector and Jacobian
965  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
966  DenseMatrix<double> &jacobian)
967  {
968  //Call the generic routine with the flag set to 1
970  jacobian,1);
971  //Call the generic FD routine to get external data
972  this->fill_in_jacobian_from_external_by_fd(jacobian);
973 
974  //Call the generic finite difference routine to handle the solid variables
975  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
976  }
977 
978  ///Set the kinematic local equation
979  inline int kinematic_local_eqn(const unsigned &n)
980  {
981  return this->nodal_local_eqn(n,this->Lagrange_index[n]);
982  }
983 
984  };
985 
986 
987 //=========================================================================
988 /// Pseudo-elasticity version of the LineFluidInterfaceBoundingElement
989 //========================================================================
990  template<class ELEMENT>
992  public FaceGeometry<FaceGeometry<ELEMENT> > ,
993  public LineFluidInterfaceBoundingElement, public virtual SolidFiniteElement
994 
995  {
996 
997  /// Short Storage for the index of Lagrange multiplier
998  Vector<unsigned> Lagrange_index;
999 
1000  public:
1001 
1002  ///\short Set the Id
1003  void set_lagrange_index(const Vector<unsigned> &lagrange_index)
1004  {Lagrange_index = lagrange_index;}
1005 
1006  /// Constructor
1008  FaceGeometry<FaceGeometry<ELEMENT> >(),
1010 
1011  /// \short Specify the value of nodal zeta from the face geometry
1012  /// The "global" intrinsic coordinate of the element when
1013  /// viewed as part of a geometric object should be given by
1014  /// the FaceElement representation, by default
1015  double zeta_nodal(const unsigned &n, const unsigned &k,
1016  const unsigned &i) const
1017  {return FaceElement::zeta_nodal(n,k,i);}
1018 
1019 
1020  /// Overload the output function
1021  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
1022 
1023  /// Output the element
1024  void output(std::ostream &outfile, const unsigned &n_plot)
1025  {
1026  FluidInterfaceBoundingElement::output(outfile,n_plot);
1027  }
1028 
1029  ///Overload the C-style output function
1030  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
1031 
1032  ///C-style Output function
1033  void output(FILE* file_pt, const unsigned &n_plot)
1034  {FluidInterfaceBoundingElement::output(file_pt,n_plot);}
1035 
1036  /// Calculate the elemental residual vector and Jacobian
1037  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
1038  DenseMatrix<double> &jacobian)
1039  {
1040  //Call the generic routine with the flag set to 1
1042 
1043  //Call the generic FD routine to get externals
1044  this->fill_in_jacobian_from_external_by_fd(jacobian);
1045 
1046  //Call the generic finite difference routine to handle the solid variables
1047  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
1048  }
1049 
1050  /// Local eqn number of kinematic bc associated with local node n
1051  int kinematic_local_eqn(const unsigned &n)
1052  {
1053  //Read out the kinematic constraint from the Id which is passed down
1054  //from the constructing element
1055  return this->nodal_local_eqn(n,this->Lagrange_index[n]);
1056  }
1057 
1058  };
1059 
1060 
1061 //==================GEOMETRIC SPECIALISATIONS==========================
1062 
1063 
1064 ///Specialise the elastic update template class to concrete 1D case
1065  template<class ELEMENT>
1066  class ElasticLineFluidInterfaceElement :
1067  public ElasticUpdateFluidInterfaceElement<FluidInterfaceElement,
1068  LineDerivatives,ELEMENT>
1069  {
1070  public:
1071 
1073  FiniteElement* const &element_pt,
1074  const int &face_index, const unsigned &id=0) :
1076  LineDerivatives,ELEMENT>(element_pt,face_index,id) {}
1077  };
1078 
1079 ///Define the BoundingElement type associated with the 1D surface element
1080  template<class ELEMENT>
1084  {
1085  public:
1086 
1090  };
1091 
1092 
1093 ///Specialise the Elastic update case to axisymmetric equations
1094  template<class ELEMENT>
1095  class ElasticAxisymmetricFluidInterfaceElement :
1096  public ElasticUpdateFluidInterfaceElement<FluidInterfaceElement,
1097  AxisymmetricDerivatives,ELEMENT>
1098  {
1099  public:
1100 
1102  FiniteElement* const &element_pt,
1103  const int &face_index, const unsigned &id=0) :
1106  element_pt,face_index,id) {}
1107  };
1108 
1109 //Define the bounding element associated with the axsymmetric elastic fluid interface element
1110  template<class ELEMENT>
1114  {
1115  public:
1116 
1120  };
1121 
1122 
1123 ///Specialise Elastic update case to the concrete 2D case
1124  template<class ELEMENT>
1126  public ElasticUpdateFluidInterfaceElement<FluidInterfaceElement,
1127  SurfaceDerivatives,ELEMENT>
1128  {
1129  public:
1130 
1132  FiniteElement* const &element_pt,
1133  const int &face_index, const unsigned &id=0) :
1135  SurfaceDerivatives,ELEMENT>(element_pt,face_index,id) {}
1136  };
1137 
1138 
1139 //Define the bounding element associated with the 2D surface elements
1140  template<class ELEMENT>
1144  {
1145  public:
1146 
1150  };
1151 
1152 
1153 
1154 }
1155 
1156 #endif
1157 
1158 
1159 
1160 
1161 
1162 
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the element's residual vector and Jacobian.
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const DShape &dpsifdS, const DShape &dpsifdS_div, const Vector< double > &s, const Vector< double > &interpolated_x, const Vector< double > &interpolated_n, const double &W, const double &J)
Helper function to calculate the additional contributions These are those filled in by the particular...
Spine version of the LineFluidInterfaceBoundingElement.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian.
Specialisation of the interface boundary constraint to a line.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the spine he...
void output(std::ostream &outfile)
Overload the output function.
Specific policy class for the FluidInterfaceElemetnts, which do not require any additional values at ...
SpineUpdateFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor, the arguments are a pointer to the "bulk" element and the index of the face to be create...
unsigned nadditional_values(const unsigned &n)
Specific interface that states how many additional values are required for the n-th node...
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element of the type specified by the BoundingElementType policy class Here...
Specialise Elastic update case to the concrete 2D case.
void fill_in_generic_residual_contribution_interface_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overload the helper function to calculate the residuals and (if flag==1) the Jacobian – this function...
double & lagrange(const unsigned &n)
Return the lagrange multiplier at local node n.
Vector< unsigned > Lagrange_index
Short Storage for the index of the Lagrange multiplier at the chosen nodes.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
Specialisation of the interface boundary constraint to a point.
int kinematic_local_eqn(const unsigned &n)
Equation number of the kinematic BC associated with node j. (This is the equation for the Lagrange mu...
int kinematic_local_eqn(const unsigned &n)
In spine elements, the kinematic condition is the equation used to determine the unknown spine height...
void output(std::ostream &outfile)
Overload the output function.
int kinematic_local_eqn(const unsigned &n)
Local eqn number of kinematic bc associated with local node n.
void set_lagrange_index(const Vector< unsigned > &lagrange_index)
Set the Id.
This policy class is used to associate specific bounding elements with specific FluidInterface elemen...
Generic Spine node update interface template class that can be combined with a given surface equation...
void setup_equation_indices(ELEMENT *const &element_pt, const unsigned &id)
Specify any additional index setup information that is required; i.e. the look-up schemes for the add...
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const DShape &dpsifdS, const DShape &dpsifdS_div, const Vector< double > &s, const Vector< double > &interpolated_x, const Vector< double > &interpolated_n, const double &W, const double &J)
Helper function to calculate the additional contributions to be added at each integration point...
ElasticUpdateFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor, pass a pointer to the bulk element and the face index of the bulk element to which the e...
void setup_equation_indices(FluidInterfaceElement *const &element_pt, const unsigned &id)
Specify any additional index setup information that is required; i.e. the look-up schemes for the add...
void output(std::ostream &outfile)
Overload the output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
void output(FILE *file_pt)
Overload the C-style output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 2D line element of type ElasticLineFluidInterfaceBoundi...
SpineLineFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
void output(FILE *file_pt)
Overload the C-style output function.
void output(std::ostream &outfile)
Overload the output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contribution to the residuals and the jacobian.
ElasticLineFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
This policy class is used to allow additional values to be added to the nodes from new surface equati...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
Pseudo-elasticity version of the LineFluidInterfaceBoundingElement.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
void fill_in_generic_residual_contribution_interface_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overload the helper function to calculate the residuals and (if flag==true) the Jacobian – this funct...
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)
Calculate the elemental residual vector and Jacobian.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the jacobian.
void output(FILE *file_pt)
Overload the C-style output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
Pseudo-elasticity version of the PointFluidInterfaceBoundingElement.
void output(std::ostream &outfile)
Overload the output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
int kinematic_local_eqn(const unsigned &n)
Set the kinematic local equation.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
int kinematic_local_eqn(const unsigned &n)
Local eqn number of the kinematic bc associated with local node n.
Generic Elastic node update interface template class that can be combined with a given surface equati...
SpineSurfaceFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
ElasticSurfaceFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the Lagrange...
void set_lagrange_index(const Vector< unsigned > &lagrange_index)
Set the Id and offset.
void output(std::ostream &outfile)
Overload the output function.
int kinematic_local_eqn(const unsigned &n)
Return local equation number associated with the kinematic constraint for local node n...
unsigned nadditional_values(const unsigned &n)
Specific interface that states how many additional values are required for the n-th node...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental residual vector and the Jacobian.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
unsigned lagrange_index(const unsigned &n)
Return the index at which the lagrange multiplier is stored at the n-th node.
Vector< unsigned > Lagrange_index
Short Storage for the index of Lagrange multiplier.
void output(FILE *file_pt)
Overload the C-style output function.
ElasticAxisymmetricFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
void output(FILE *file_pt)
Overload the C-style output function.
void output(FILE *file_pt)
Overload the C-style output function.
double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &surface_gradient, DShape &surface_divergence)
Fill in the specific surface derivative calculations by calling the appropriate class function...
double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &surface_gradient, DShape &surface_divergence)
Fill in the specific surface derivative calculations by calling the appropriate function from the der...
Vector< unsigned > Lagrange_index
Storage for the location of the Lagrange multiplier (If other additional values have been added we ne...
SpineAxisymmetricFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
Spine version of the PointFluidInterfaceBoundingElement.
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 ...