surface_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_SURFACE_INTERFACE_ELEMENTS_HEADER
34 #define OOMPH_SURFACE_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 /// Two-dimensional interface elements that are used with a spine mesh,
52 /// i.e. the mesh deformation is handled by Kistler & Scriven's "method
53 /// of spines". These elements are FaceElements are attached to 3D bulk
54 /// Fluid elements and the particular type of fluid element is passed
55 /// as a template parameter to the element. It
56 /// shouldn't matter whether the passed
57 /// element is the underlying (fixed) element or the templated
58 /// SpineElement<Element>.
59 /// Optionally, an external pressure may be specified, which must be
60 /// passed to the element as external data. If there is no such object,
61 /// the external pressure is assumed to be zero.
62 //======================================================================
63  template<class ELEMENT>
64  class SpineSurfaceFluidInterfaceElement :
65  public virtual Hijacked<SpineElement<FaceGeometry<ELEMENT> > >,
66  public virtual SurfaceFluidInterfaceElement
67  {
68  private:
69 
70  /// \short In spine elements, the kinematic condition is the equation
71  /// used to determine the unknown spine heights. Overload the
72  /// function accordingly
73  int kinematic_local_eqn(const unsigned &n)
74  {return this->spine_local_eqn(n);}
75 
76  /// \short Hijacking the kinematic condition corresponds to hijacking the
77  /// variables associated with thespine heights.
78  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
79  {
80  //Loop over all the node numbers that are passed in
81  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
82  it!=bulk_node_number.end();++it)
83  {
84  //Hijack the spine heights. (and delete the returned data object)
85  delete this->hijack_nodal_spine_value(*it,0);
86  }
87  }
88 
89  public:
90 
91 
92  /// \short Constructor, the arguments are a pointer to the "bulk" element
93  /// and the index of the face to be created
94  SpineSurfaceFluidInterfaceElement(FiniteElement* const &element_pt,
95  const int &face_index) :
96  Hijacked<SpineElement<FaceGeometry<ELEMENT> > >(),
97  SurfaceFluidInterfaceElement()
98  {
99  //Attach the geometrical information to the element, by
100  //making the face element from the bulk element
101  element_pt->build_face_element(face_index,this);
102 
103 #ifdef PARANOID
104  //Is it refineable
105  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(element_pt);
106  if(ref_el_pt!=0)
107  {
108  if (this->has_hanging_nodes())
109  {
110  throw OomphLibError(
111  "This interface element will not work correctly if nodes are hanging\n",
112  OOMPH_CURRENT_FUNCTION,
113  OOMPH_EXCEPTION_LOCATION);
114  }
115  }
116 #endif
117 
118  //Find the index at which the velocity unknowns are stored
119  //from the bulk element
120  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
121  this->U_index_interface.resize(3);
122  for(unsigned i=0;i<3;i++)
123  {
124  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
125  }
126  }
127 
128  /// Calculate the contribution to the residuals and the jacobian
129  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
130  DenseMatrix<double> &jacobian)
131  {
132  //Call the generic routine with the flag set to 1
134 
135  //Call the generic routine to handle the shape derivatives
136  this->fill_in_jacobian_from_geometric_data(jacobian);
137  }
138 
139  /// \short
140  /// Helper function to calculate the additional contributions
141  /// to be added at each integration point. Empty for this
142  /// implemenetation
144  Vector<double> &residuals,
145  DenseMatrix<double> &jacobian,
146  const unsigned &flag,
147  const Shape &psif,
148  const DShape &dpsifds,
149  const Vector<double> &interpolated_x,
150  const Vector<double> &interpolated_n,
151  const double &W,
152  const double &J)
153  {
154  }
155 
156  /// Overload the output function
157  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
158 
159  /// Output the element
160  void output(std::ostream &outfile, const unsigned &n_plot)
161  {SurfaceFluidInterfaceElement::output(outfile,n_plot);}
162 
163  ///Overload the C-style output function
164  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
165 
166  ///C-style Output function
167  void output(FILE* file_pt, const unsigned &n_plot)
168  {SurfaceFluidInterfaceElement::output(file_pt,n_plot);}
169 
170 
171 
172  /// \short Create an "bounding" element (here actually a 2D line element
173  /// of type SpineLineFluidInterfaceBoundingElement<ELEMENT> that allows
174  /// the application of a contact angle boundary condition on the
175  /// the specified face.
177  const int &face_index)
178  {
179  //Create a temporary pointer to the appropriate FaceElement
182 
183  //Attach the geometrical information to the new element
184  this->build_face_element(face_index,face_el_pt);
185 
186  //Set the index at which the velocity nodes are stored
187  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
188 
189  //Set the value of the nbulk_value, the node is not resized
190  //in this bounding element,
191  //so it will just be the actual nvalue here
192  // There is some ambiguity about what this means (however)
193  const unsigned n_node_bounding = face_el_pt->nnode();
194  for(unsigned n=0;n<n_node_bounding;n++)
195  {
196  face_el_pt->nbulk_value(n) =
197  face_el_pt->node_pt(n)->nvalue();
198  }
199 
200  //Set of unique geometric data that is used to update the bulk,
201  //but is not used to update the face
202  std::set<Data*> unique_additional_geom_data;
203 
204  //Get all the geometric data for this (bulk) element
205  this->assemble_set_of_all_geometric_data(unique_additional_geom_data);
206 
207  //Now assemble the set of geometric data for the face element
208  std::set<Data*> unique_face_geom_data_pt;
209  face_el_pt->assemble_set_of_all_geometric_data(unique_face_geom_data_pt);
210 
211  //Erase the face geometric data from the additional data
212  for(std::set<Data*>::iterator it=unique_face_geom_data_pt.begin();
213  it!=unique_face_geom_data_pt.end();++it)
214  {unique_additional_geom_data.erase(*it);}
215 
216  //Finally add all unique additional data as external data
217  for(std::set<Data*>::iterator it = unique_additional_geom_data.begin();
218  it!= unique_additional_geom_data.end();++it)
219  {
220  face_el_pt->add_external_data(*it);
221  }
222 
223  //Return the value of the pointer
224  return face_el_pt;
225  }
226 
227  };
228 
229 
230 
231 ///////////////////////////////////////////////////////////////////////
232 ///////////////////////////////////////////////////////////////////////
233 ///////////////////////////////////////////////////////////////////////
234 
235 
236 
237 //=======================================================================
238 /// Two-dimensional interface elements that are used when the mesh
239 /// deformation is handled by a set of equations that modify the nodal
240 /// positions. These elements are FaceElements attached to 3D bulk fluid
241 /// elements and the fluid element is passed as a template parameter to
242 /// the element.
243 /// Optionally an external pressure may be specified, which must be
244 /// passed to the element as external data. The default value of the external
245 /// pressure is zero.
246 //=======================================================================
247  template<class ELEMENT>
248  class ElasticSurfaceFluidInterfaceElement :
249  public virtual Hijacked<FaceGeometry<ELEMENT> >,
250  public SurfaceFluidInterfaceElement
251  {
252  private:
253 
254  /// \short ID of the Lagrange Lagrange multiplier (in the collection of
255  /// nodal values accomodated by resizing)
256  unsigned Id;
257 
258  /// \short Equation number of the kinematic BC associated with node j.
259  /// (This is the equation for the Lagrange multiplier)
260  int kinematic_local_eqn(const unsigned &j)
261  {
262  // Get the index of the nodal value associated with Lagrange multiplier
263  const unsigned lagr_index=
264  dynamic_cast<BoundaryNodeBase*>(this->node_pt(j))->
265  index_of_first_value_assigned_by_face_element(Id);
266 
267  return this->nodal_local_eqn(j,lagr_index);
268  }
269 
270 
271  /// \short Hijacking the kinematic condition corresponds to hijacking the
272  /// variables associated with the Lagrange multipliers that are assigned
273  /// on construction of this element.
274  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
275  {
276  //Loop over all the nodes that are passed in
277  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
278  it!=bulk_node_number.end();++it)
279  {
280  //Get the index associated with the Id for each node
281  //(the Lagrange multiplier)
282  unsigned n_lagr = dynamic_cast<BoundaryNodeBase*>(node_pt(*it))->
283  index_of_first_value_assigned_by_face_element(Id);
284 
285  //Hijack the appropriate value and delete the returned Node
286  delete this->hijack_nodal_value(*it,n_lagr);
287  }
288  }
289 
290  public:
291 
292 
293 
294  /// \short Constructor, pass a pointer to the bulk element and the face
295  /// index of the bulk element to which the element is to be attached to.
296  /// The optional identifier can be used
297  /// to distinguish the additional nodal value (Lagr mult) created by
298  /// this element from those created by other FaceElements.
299  ElasticSurfaceFluidInterfaceElement(FiniteElement* const &element_pt,
300  const int &face_index,
301  const unsigned &id=0) :
302  FaceGeometry<ELEMENT>(), SurfaceFluidInterfaceElement(), Id(id)
303  {
304  //Attach the geometrical information to the element
305  //This function also assigned nbulk_value from required_nvalue of the
306  //bulk element
307  element_pt->build_face_element(face_index,this);
308 
309 #ifdef PARANOID
310  //Is it refineable
311  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(element_pt);
312  if(ref_el_pt!=0)
313  {
314  if (this->has_hanging_nodes())
315  {
316  throw OomphLibError(
317  "This flux element will not work correctly if nodes are hanging\n",
318  OOMPH_CURRENT_FUNCTION,
319  OOMPH_EXCEPTION_LOCATION);
320  }
321  }
322 #endif
323 
324  //Find the index at which the velocity unknowns are stored
325  //from the bulk element
326  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
327  this->U_index_interface.resize(3);
328  for(unsigned i=0;i<3;i++)
329  {
330  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
331  }
332 
333  //Read out the number of nodes on the face
334  unsigned n_node_face = this->nnode();
335 
336  //Set the additional data values in the face
337  //There is one additional values at each node --- the Lagrange multiplier
338  Vector<unsigned> additional_data_values(n_node_face);
339  for(unsigned i=0;i<n_node_face;i++) additional_data_values[i] = 1;
340 
341  // Now add storage for Lagrange multipliers and set the map containing
342  // the position of the first entry of this face element's
343  // additional values.
344  add_additional_values(additional_data_values,id);
345  }
346 
347  /// \short The "global" intrinsic coordinate of the element when
348  /// viewed as part of a geometric object should be given by
349  /// the FaceElement representation, by default
350  double zeta_nodal(const unsigned &n, const unsigned &k,
351  const unsigned &i) const
352  {return FaceElement::zeta_nodal(n,k,i);}
353 
354  /// Return the lagrange multiplier at local node j
355  double &lagrange(const unsigned &j)
356  {
357  // Get the index of the nodal value associated with Lagrange multiplier
358  unsigned lagr_index=dynamic_cast<BoundaryNodeBase*>(node_pt(j))->
359  index_of_first_value_assigned_by_face_element(Id);
360  return *node_pt(j)->value_pt(lagr_index);
361  }
362 
363 
364  /// Fill in contribution to residuals and Jacobian
365  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
366  DenseMatrix<double> &jacobian)
367  {
368  //Call the generic routine with the flag set to 1
370 
371  //Call the generic finite difference routine for the solid variables
372  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
373  }
374 
375  /// Overload the output function
376  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
377 
378  /// Output the element
379  void output(std::ostream &outfile, const unsigned &n_plot)
380  {SurfaceFluidInterfaceElement::output(outfile,n_plot);}
381 
382  ///Overload the C-style output function
383  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
384 
385  ///C-style Output function
386  void output(FILE* file_pt, const unsigned &n_plot)
387  {SurfaceFluidInterfaceElement::output(file_pt,n_plot);}
388 
389 
390  /// \short Helper function to calculate the additional contributions
391  /// to be added at each integration point. This deals with
392  /// Lagrange multiplier contribution
394  Vector<double> &residuals,
395  DenseMatrix<double> &jacobian,
396  const unsigned &flag,
397  const Shape &psif,
398  const DShape &dpsifds,
399  const Vector<double> &interpolated_x,
400  const Vector<double> &interpolated_n,
401  const double &W,
402  const double &J)
403  {
404  //Assemble Lagrange multiplier by loop over the shape functions
405  unsigned n_node = this->nnode();
406  double interpolated_lagrange = 0.0;
407  for(unsigned l=0;l<n_node;l++)
408  {
409  //Note same shape functions used for lagrange multiplier field
410  interpolated_lagrange += lagrange(l)*psif(l);
411  }
412 
413  int local_eqn=0, local_unknown = 0;
414 
415  //Loop over the shape functions
416  for(unsigned l=0;l<n_node;l++)
417  {
418  //Loop over the directions
419  for(unsigned i=0;i<3;i++)
420  {
421  //Now using the same shape functions for the elastic equations,
422  //so we can stay in the loop
423  local_eqn = this->position_local_eqn(l,0,i);
424  if(local_eqn >= 0)
425  {
426  //Add in the Lagrange multiplier contribution
427  residuals[local_eqn] -=
428  interpolated_lagrange*interpolated_n[i]*psif(l)*J*W;
429 
430  //Do the Jacobian calculation
431  if(flag)
432  {
433  //Loop over the nodes
434  for(unsigned l2=0;l2<n_node;l2++)
435  {
436  // Dependence on solid positions will be handled by FDs
437  //That leaves the Lagrange multiplier contribution
438  local_unknown = kinematic_local_eqn(l2);
439  if(local_unknown >= 0)
440  {
441  jacobian(local_eqn,local_unknown) -=
442  psif(l2)*interpolated_n[i]*psif(l)*J*W;
443  }
444  }
445  } //End of Jacobian calculation
446  }
447  }
448 
449  }
450  }
451 
452 
453 
454  /// \short Create an "bounding" element (here actually a 2D line element
455  /// of type ElasticLineFluidInterfaceBoundingElement<ELEMENT> that allows
456  /// the application of a contact angle boundary condition on the
457  /// the specified face.
459  const int &face_index)
460  {
461  //Create a temporary pointer to the appropriate FaceElement
464 
465  //Attach the geometrical information to the new element
466  this->build_face_element(face_index,face_el_pt);
467 
468  //Set the index at which the velocity nodes are stored
469  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
470 
471  //Set the value of the nbulk_value, the node is not resized
472  //in this bounding element,
473  //so it will just be the actual nvalue here
474  // There is some ambiguity about what this means (however)
475  const unsigned n_node_bounding = face_el_pt->nnode();
476  for(unsigned n=0;n<n_node_bounding;n++)
477  {
478  face_el_pt->nbulk_value(n) =
479  face_el_pt->node_pt(n)->nvalue();
480  }
481 
482  //Pass the ID down
483  face_el_pt->set_id(Id);
484 
485  //Find the nodes
486  std::set<SolidNode*> set_of_solid_nodes;
487  const unsigned n_node = this->nnode();
488  for(unsigned n=0;n<n_node;n++)
489  {
490  set_of_solid_nodes.insert(static_cast<SolidNode*>(this->node_pt(n)));
491  }
492 
493  //Delete the nodes from the face
494  //n_node = face_el_pt->nnode();
495  for(unsigned n=0;n<n_node_bounding;n++)
496  {
497  //Set the value of the nbulk_value, from the present element
498  face_el_pt->nbulk_value(n) =
499  this->nbulk_value(face_el_pt->bulk_node_number(n));
500 
501  //Now delete the nodes from the set
502  set_of_solid_nodes.erase(static_cast<SolidNode*>(
503  face_el_pt->node_pt(n)));
504  }
505 
506  //Now add these as external data
507  for(std::set<SolidNode*>::iterator it=set_of_solid_nodes.begin();
508  it!=set_of_solid_nodes.end();++it)
509  {
510  face_el_pt->add_external_data((*it)->variable_position_pt());
511  }
512 
513 
514  //Return the value of the pointer
515  return face_el_pt;
516  }
517 
518  };
519 
520 }
521 
522 #endif
523 
524 
525 
526 
527 
528 
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the Lagrange...
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, 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...
void output(FILE *file_pt)
Overload the C-style output function.
Spine version of the LineFluidInterfaceBoundingElement.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with thespine hei...
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, 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...
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 ...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian.
int kinematic_local_eqn(const unsigned &j)
Equation number of the kinematic BC associated with node j. (This is the equation for the Lagrange mu...
double & lagrange(const unsigned &j)
Return the lagrange multiplier at local node j.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt)
Overload the C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 2D line element of type SpineLineFluidInterfaceBounding...
Pseudo-elasticity version of the LineFluidInterfaceBoundingElement.
void output(std::ostream &outfile)
Overload the output function.
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 2D line element of type ElasticLineFluidInterfaceBoundi...
SpineSurfaceFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, the arguments are a pointer to the "bulk" element and the index of the face to be create...
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
ElasticSurfaceFluidInterfaceElement(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...
Vector< unsigned > & u_index_interface_boundary()
Access for nodal index at which the velocity components are stored.
int kinematic_local_eqn(const unsigned &n)
In spine elements, the kinematic condition is the equation used to determine the unknown spine height...
Vector< unsigned > U_index_interface
Nodal index at which the i-th velocity component is stored.
unsigned Id
ID of the Lagrange Lagrange multiplier (in the collection of nodal values accomodated by resizing) ...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contribution to the residuals and the jacobian.
virtual void fill_in_generic_residual_contribution_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Helper function to calculate the residuals and (if flag==1) the Jacobian of the equations. This is implemented generically using the surface divergence information that is overloaded in each element i.e. axisymmetric, two- or three-dimensional.