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
79  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
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
114  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
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,
134  const Vector<double> &interpolated_x,
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)
147  {AxisymmetricFluidInterfaceElement::output(outfile,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)
154  {AxisymmetricFluidInterfaceElement::output(file_pt,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.
251  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
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
332  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
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)
347  {AxisymmetricFluidInterfaceElement::output(outfile,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)
354  {AxisymmetricFluidInterfaceElement::output(file_pt,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,
365  const Vector<double> &interpolated_x,
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.
Specialise the Elastic update case to axisymmetric equations.
void output(std::ostream &outfile, const unsigned &n_plot)
Output the element.
void output(std::ostream &outfile)
Overload the output function.
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)
Fill in contribution to residuals and Jacobian.
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 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...
void hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the Lagrange...
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...
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...
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.
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.
Spine version of the PointFluidInterfaceBoundingElement.