axisym_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 free surface/interface elements in an axisymmetric geometry
31 
32 //Include guards to prevent multiple inclusions of the file
33 #ifndef OOMPH_AXISYMMETRIC_INTERFACE_ELEMENTS_HEADER
34 #define OOMPH_AXISYMMETRIC_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 header files
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 /// Axisymmetric 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 attached to 2D 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>
65  public virtual Hijacked<SpineElement<FaceGeometry<ELEMENT> > >,
66  public virtual AxisymmetricFluidInterfaceElement
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 
77  /// \short Hijacking the kinematic condition corresponds to hijacking the
78  /// spine heights -- used for strong imposition of contact angle condition
80  {
81  //Loop over all the node numbers that are passed in
82  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
83  it!=bulk_node_number.end();++it)
84  {
85  //Hijack the spine heights and delete the returned data objects
86  delete this->hijack_nodal_spine_value(*it,0);
87  }
88  }
89 
90  public:
91 
92  /// \short Constructor, the arguments are a pointer to the "bulk" element
93  /// and the face index
95  FiniteElement* const &element_pt, const int &face_index) :
96  Hijacked<SpineElement<FaceGeometry<ELEMENT> > >(),
97  AxisymmetricFluidInterfaceElement()
98  {
99  //Attach the geometrical information to the element
100  element_pt->build_face_element(face_index,this);
101 
102  //Find the index at which the velocity unknowns are stored
103  //from the bulk element
104  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
105  this->U_index_interface.resize(3);
106  for(unsigned i=0;i<3;i++)
107  {
108  this->U_index_interface[i] =
109  cast_element_pt->u_index_nst(i);
110  }
111  }
112 
113  /// Calculate the contribution to the residuals and the jacobian
115  DenseMatrix<double> &jacobian)
116  {
117  //Call the generic residuals routine with the flag set to 1
119 
120  //Call the generic routine to evaluate shape derivatives
121  this->fill_in_jacobian_from_geometric_data(jacobian);
122  }
123 
124  /// \short
125  /// Helper function to calculate the additional contributions
126  /// to be added at each integration point. Empty for this
127  /// implemenetation
129  Vector<double> &residuals,
130  DenseMatrix<double> &jacobian,
131  const unsigned &flag,
132  const Shape &psif,
133  const DShape &dpsifds,
135  const Vector<double> &interpolated_n,
136  const double &W,
137  const double &J)
138  {
139  }
140 
141 
142  /// Overload the output function
143  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
144 
145  /// Output the element
146  void output(std::ostream &outfile, const unsigned &n_plot)
148 
149  ///Overload the C-style output function
150  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
151 
152  ///C-style Output function
153  void output(FILE* file_pt, const unsigned &n_plot)
155 
156 
157  /// \short Create an "bounding" element (here actually a 1D point element
158  /// of type SpinePointFluidInterfaceBoundingElement<ELEMENT>) that allows
159  /// the application of a contact angle boundary condition on the
160  /// the specified face.
162  const int& face_index)
163  {
164  //Create a temporary pointer to the appropriate FaceElement
167 
168  //Attach geometric information to the new element
169  this->build_face_element(face_index,face_el_pt);
170 
171  //Set the index at which the unknowns are stored from the element
172  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
173 
174  //Set the value of the nbulk_value, the node is not resized
175  //in this problem, so it will just be the actual nvalue
176  face_el_pt->nbulk_value(0) = face_el_pt->node_pt(0)->nvalue();
177 
178  //Set of unique geometric data that is used to update the bulk,
179  //but is not used to update the face
180  std::set<Data*> unique_additional_geom_data;
181 
182  //Get all the geometric data for this (bulk) element
183  this->assemble_set_of_all_geometric_data(unique_additional_geom_data);
184 
185  //Now assemble the set of geometric data for the face element
186  std::set<Data*> unique_face_geom_data_pt;
187  face_el_pt->assemble_set_of_all_geometric_data(unique_face_geom_data_pt);
188 
189  //Erase the face geometric data from the additional data
190  for(std::set<Data*>::iterator it=unique_face_geom_data_pt.begin();
191  it!=unique_face_geom_data_pt.end();++it)
192  {unique_additional_geom_data.erase(*it);}
193 
194  //Finally add all unique additional data as external data
195  for(std::set<Data*>::iterator it = unique_additional_geom_data.begin();
196  it!= unique_additional_geom_data.end();++it)
197  {
198  face_el_pt->add_external_data(*it);
199  }
200 
201  //Return the value of the pointer
202  return face_el_pt;
203  }
204 
205  };
206 
207 
208 
209 ////////////////////////////////////////////////////////////////////////
210 ////////////////////////////////////////////////////////////////////////
211 ////////////////////////////////////////////////////////////////////////
212 
213 
214 
215 //========================================================================
216 /// One-dimensional interface elements that are used when the mesh
217 /// deformation is handled by a set of equations that modify the nodal
218 /// positions. These elements are FaceElements attached to axisymmetric
219 /// bulk fluid elements and the fluid element is passed as a
220 /// template parameter to the element.
221 /// Optionally an external pressure may be specified, which must be
222 /// passed to the element as external data. The default value of the external
223 /// pressure is zero.
224 //======================================================================
225  template<class ELEMENT>
227  public virtual Hijacked<FaceGeometry<ELEMENT> >,
228  public AxisymmetricFluidInterfaceElement
229  {
230  protected:
231 
232  /// \short ID of the Lagrange Lagrange multiplier (in the collection of
233  /// nodal values accomodated by resizing)
234  unsigned Id;
235 
236  /// \short Equation number of the kinematic BC associated with node j.
237  /// (This is the equation for the Lagrange multiplier)
238  int kinematic_local_eqn(const unsigned &j)
239  {
240  // Get the index of the nodal value associated with Lagrange multiplier
241  unsigned lagr_index=dynamic_cast<BoundaryNodeBase*>(node_pt(j))->
242  index_of_first_value_assigned_by_face_element(Id);
243 
244  // Return nodal value
245  return this->nodal_local_eqn(j,lagr_index);
246  }
247 
248  /// \short Hijacking the kinematic condition corresponds to hijacking the
249  /// variables associated with the Lagrange multipliers that are assigned
250  /// on construction of this element.
252  {
253 
254  //Loop over all the passed nodes and delete the returned values
255  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
256  it!=bulk_node_number.end();++it)
257  {
258  //Get the index associated with the Id for each node
259  //(the Lagrange multiplier)
260  unsigned n_lagr = dynamic_cast<BoundaryNodeBase*>(node_pt(*it))->
261  index_of_first_value_assigned_by_face_element(Id);
262 
263  //Hijack the appropriate value and delete the returned Node
264  delete this->hijack_nodal_value(*it,n_lagr);
265  }
266  }
267 
268  public:
269 
270  /// \short The "global" intrinsic coordinate of the element when
271  /// viewed as part of a geometric object should be given by
272  /// the FaceElement representation, by default
273  double zeta_nodal(const unsigned &n, const unsigned &k,
274  const unsigned &i) const
275  {return FaceElement::zeta_nodal(n,k,i);}
276 
277 
278  /// \short Constructor, pass a pointer to the bulk element
279  /// and the face index
280  /// of the bulk element to which the element is to be attached to.
281  /// The optional identifier can be used
282  /// to distinguish the additional nodal value (Lagr mult) created by
283  /// this element from those created by other FaceElements.
285  FiniteElement* const &element_pt,
286  const int &face_index,
287  const unsigned &id=0) : FaceGeometry<ELEMENT> (),
288  AxisymmetricFluidInterfaceElement(), Id(id)
289  {
290  //Attach the geometrical information to the element
291  //This function also assigned nbulk_value from required_nvalue of the
292  //bulk element
293  element_pt->build_face_element(face_index,this);
294 
295  //Find the index at which the velocity unknowns are stored
296  //from the bulk element
297  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
298  this->U_index_interface.resize(3);
299  for(unsigned i=0;i<3;i++)
300  {
301  this->U_index_interface[i] = cast_element_pt->
302  u_index_nst(i);
303  }
304 
305  //Read out the number of nodes on the face
306  //For some reason I need to specify the this pointer here(!)
307  unsigned n_node_face = this->nnode();
308 
309  //Set the additional data values in the face
310  //There is one additional values at each node, in this case
311  Vector<unsigned> additional_data_values(n_node_face);
312  for(unsigned i=0;i<n_node_face;i++) additional_data_values[i] = 1;
313 
314  // Now add storage for Lagrange multipliers and set the map containing
315  // the position of the first entry of this face element's
316  // additional values.
317  add_additional_values(additional_data_values,id);
318  }
319 
320 
321  /// Return the Lagrange multiplier at local node j
322  double &lagrange(const unsigned &j)
323  {
324  // Get the index of the nodal value associated with Lagrange multiplier
325  unsigned lagr_index=dynamic_cast<BoundaryNodeBase*>(node_pt(j))->
326  index_of_first_value_assigned_by_face_element(Id);
327  return *node_pt(j)->value_pt(lagr_index);
328  }
329 
330 
331  /// Fill in contribution to residuals and Jacobian
333  DenseMatrix<double> &jacobian)
334  {
335  //Call the generic routine with the flag set to 1
337 
338  //Call the generic finite difference routine for the solid variables
339  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
340  }
341 
342  /// Overload the output function
343  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
344 
345  /// Output the element
346  void output(std::ostream &outfile, const unsigned &n_plot)
348 
349  ///Overload the C-style output function
350  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
351 
352  ///C-style Output function
353  void output(FILE* file_pt, const unsigned &n_plot)
355 
356  /// \short Helper function to calculate the additional contributions
357  /// to be added at each integration point. This deals with
358  /// Lagrange multiplier contribution.
360  Vector<double> &residuals,
361  DenseMatrix<double> &jacobian,
362  const unsigned &flag,
363  const Shape &psif,
364  const DShape &dpsifds,
366  const Vector<double> &interpolated_n,
367  const double &W, const double &J)
368  {
369  //Loop over the shape functions
370  unsigned n_node = this->nnode();
371 
372  //Assemble Lagrange multiplier by loop over the shape functions
373  double interpolated_lagrange = 0.0;
374  for(unsigned l=0;l<n_node;l++)
375  {
376  //Note same shape functions used for lagrange multiplier field
377  interpolated_lagrange += lagrange(l)*psif[l];
378  }
379 
380  //Cache the radial position
381  const double r = interpolated_x[0];
382 
383  int local_eqn=0, local_unknown = 0;
384  //Loop over the shape functions
385  for(unsigned l=0;l<n_node;l++)
386  {
387  //Loop over the directions
388  for(unsigned i=0;i<2;i++)
389  {
390  //Now using the same shape functions for the elastic equations,
391  //so we can stay in the loop
392  local_eqn = this->position_local_eqn(l,0,i);
393  if(local_eqn >= 0)
394  {
395  //Add in Lagrange multiplier contribution
396  residuals[local_eqn] -=
397  r*interpolated_lagrange*
398  interpolated_n[i]*psif[l]*W*J;
399 
400  //Do the Jacobian calculation
401  if(flag)
402  {
403  //Loop over the nodes
404  for(unsigned l2=0;l2<n_node;l2++)
405  {
406  // Dependence upon solid coordinates will be handled by FDs
407  //That leaves the Lagrange multiplier contribution
408  local_unknown = kinematic_local_eqn(l2);
409  if(local_unknown >= 0)
410  {
411  jacobian(local_eqn,local_unknown) -=
412  r*psif[l2]*interpolated_n[i]*psif[l]*W*J;
413  }
414  }
415  } //End of Jacobian calculation
416  }
417  }
418 
419  }
420  }
421 
422 
423  /// \short Create an "bounding" element (here actually a 1D point element
424  /// of type ElasticPointFluidInterfaceBoundingElement<ELEMENT>) that allows
425  /// the application of a contact angle boundary condition on the
426  /// the specified face.
428  const int& face_index)
429  {
430  //Create a temporary pointer to the appropriate FaceElement
433 
434  //Attach the geometrical information to the new element
435  this->build_face_element(face_index,face_el_pt);
436 
437  //Set the index at which the unknowns are stored from the element
438  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
439 
440  //Set the value of the nbulk_value, the node is not resized
441  //in this problem, so it will just be the actual nvalue - 1
442  face_el_pt->nbulk_value(0) = face_el_pt->node_pt(0)->nvalue();
443 
444  //Pass down the ID of the Lagrange multiplier
445  face_el_pt->set_id(Id);
446 
447  //Find the nodes
448  std::set<SolidNode*> set_of_solid_nodes;
449  unsigned n_node = this->nnode();
450  for(unsigned n=0;n<n_node;n++)
451  {
452  set_of_solid_nodes.insert(static_cast<SolidNode*>(this->node_pt(n)));
453  }
454 
455  //Delete the nodes from the face
456  n_node = face_el_pt->nnode();
457  for(unsigned n=0;n<n_node;n++)
458  {
459  set_of_solid_nodes.erase(static_cast<SolidNode*>(
460  face_el_pt->node_pt(n)));
461  }
462 
463  //Now add these as external data
464  for(std::set<SolidNode*>::iterator it=set_of_solid_nodes.begin();
465  it!=set_of_solid_nodes.end();++it)
466  {
467  face_el_pt->add_external_data((*it)->variable_position_pt());
468  }
469 
470 
471  //Return the value of the pointer
472  return face_el_pt;
473  }
474  };
475 
476 }
477 
478 #endif
479 
480 
481 
482 
483 
484 
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 1D point element of type SpinePointFluidInterfaceBoundi...
unsigned Id
ID of the Lagrange Lagrange multiplier (in the collection of nodal values accomodated by resizing) ...
void output(FILE *file_pt)
Overload the C-style output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
Data * hijack_nodal_value(const unsigned &n, const unsigned &i)
Hijack the i-th value stored at node n. Return a custom-made (copied) data object that contains only ...
Specialise the Elastic update case to axisymmetric equations.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void output(std::ostream &outfile)
Overload the output function.
int spine_local_eqn(const unsigned &n)
Return the local equation number corresponding to the height of the spine at the n-th node...
Definition: spines.h:467
A general Finite Element class.
Definition: elements.h:1271
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
void output(FILE *file_pt, const unsigned &n_plot)
C-style Output function.
unsigned & bulk_node_number(const unsigned &n)
Return the bulk node number that corresponds to the n-th local node number.
Definition: elements.h:4539
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian.
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
void add_additional_values(const Vector< unsigned > &nadditional_values, const unsigned &id)
Helper function adding additional values for the unknowns associated with the FaceElement. This function also sets the map containing the position of the first entry of this face element's additional values.The inputs are the number of additional values and the face element's ID. Note the same number of additonal values are allocated at ALL nodes.
Definition: elements.h:4181
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the spine heights – used for strong imposi...
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 1D point element of type ElasticPointFluidInterfaceBoun...
void output(FILE *file_pt)
Overload the C-style output function.
void output(std::ostream &outfile)
Overload the output function.
Pseudo-elasticity version of the PointFluidInterfaceBoundingElement.
int kinematic_local_eqn(const unsigned &n)
In spine elements, the kinematic condition is the equation used to determine the unknown spine height...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the Lagrange...
void output(std::ostream &outfile)
Output with default number of plot points.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
Vector< unsigned > & u_index_interface_boundary()
Access for nodal index at which the velocity components are stored.
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)
Calculate the contribution to the residuals and the 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...
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...
The SpineElement<ELEMENT> class takes an existing element as a template parameter and adds the necess...
Definition: spines.h:423
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void assemble_set_of_all_geometric_data(std::set< Data * > &unique_geom_data_pt)
Return a set of all geometric data associated with the element.
Hijacked elements are elements in which one or more Data values that affect the element's residuals...
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void 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...
ElasticAxisymmetricFluidInterfaceElement(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...
unsigned & nbulk_value(const unsigned &n)
Return the number of values originally stored at local node n (before the FaceElement added additiona...
Definition: elements.h:4552
Vector< unsigned > U_index_interface
Nodal index at which the i-th velocity component is stored.
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.
void fill_in_jacobian_from_geometric_data(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the Jacobian matrix from the geometric data. This version assumes that...
SpineAxisymmetricFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, the arguments are a pointer to the "bulk" element and the face index.
double & lagrange(const unsigned &j)
Return the Lagrange multiplier at local node j.
Data * hijack_nodal_spine_value(const unsigned &n, const unsigned &i)
Hijack the i-th value stored at the spine that affects local node n. Return a custom-made (copied) da...
Spine version of the PointFluidInterfaceBoundingElement.