line_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 (one-dimensional) free surface elements
31 //Include guards, to prevent multiple includes
32 #ifndef OOMPH_LINE_INTERFACE_ELEMENTS_HEADER
33 #define OOMPH_LINE_INTERFACE_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 //OOMPH-LIB headers
41 #include "../generic/Qelements.h"
42 #include "../generic/spines.h"
43 #include "../generic/hijacked_elements.h"
44 #include "interface_elements.h"
45 
46 namespace oomph
47 {
48 
49 //======================================================================
50 /// One-dimensional fluid interface elements that are used with a spine mesh,
51 /// i.e. the mesh deformation is handled by Kistler & Scriven's "method
52 /// of spines". These elements are FaceElements attached to 2D bulk
53 /// fluid elements and the particular type of fluid element is passed
54 /// as a template parameter to the element. It
55 /// shouldn't matter whether the passed
56 /// element is the underlying (fixed) element or the templated
57 /// SpineElement<Element>.
58 /// Optionally, an external pressure may be specified, which must be
59 /// passed to the element as external data. If there is no such object,
60 /// the external pressure is assumed to be zero.
61 //======================================================================
62 template<class ELEMENT>
64 public virtual Hijacked<SpineElement<FaceGeometry<ELEMENT> > >,
65 public virtual LineFluidInterfaceElement
66 {
67  private:
68 
69  /// \short In spine elements, the kinematic condition is the equation
70  /// used to determine the unknown spine heights. Overload the
71  /// function accordingly
72  int kinematic_local_eqn(const unsigned &n)
73  {return this->spine_local_eqn(n);}
74 
75 
76  /// \short Hijacking the kinematic condition corresponds to hijacking the
77  /// variables associated with the spine heights --- only used when
78  /// strong imposition of the contact angle condition is required.
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 object)
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 index of the face to be created.
94  SpineLineFluidInterfaceElement(FiniteElement* const &element_pt,
95  const int &face_index) :
96  Hijacked<SpineElement<FaceGeometry<ELEMENT> > >(),
97  LineFluidInterfaceElement()
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  //Find the index at which the velocity unknowns are stored
104  //from the bulk element
105  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
106  this->U_index_interface.resize(2);
107  for(unsigned i=0;i<2;i++)
108  {
109  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
110  }
111  } //End of constructor
112 
113 
114 
115  /// Calculate the contribution to the residuals and the jacobian
116  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
117  DenseMatrix<double> &jacobian)
118  {
119  //Call the generic routine with the flag set to 1
121 
122  //Call the generic routine to evaluate shape derivatives
123  this->fill_in_jacobian_from_geometric_data(jacobian);
124  } //End of jacobian contribution
125 
126 
127  /// \short
128  /// Helper function to calculate the additional contributions
129  /// to be added at each integration point. Empty for this
130  /// implemenetation
132  Vector<double> &residuals,
133  DenseMatrix<double> &jacobian,
134  const unsigned &flag,
135  const Shape &psif,
136  const DShape &dpsifds,
137  const Vector<double> &interpolated_x,
138  const Vector<double> &interpolated_n,
139  const double &W,
140  const double &J)
141  {
142  }
143 
144  /// Overload the output function
145  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
146 
147  /// Output the element
148  void output(std::ostream &outfile, const unsigned &n_plot)
149  {LineFluidInterfaceElement::output(outfile,n_plot);}
150 
151  ///Overload the C-style output function
152  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
153 
154  ///C-style Output function
155  void output(FILE* file_pt, const unsigned &n_plot)
156  {LineFluidInterfaceElement::output(file_pt,n_plot);}
157 
158 
159  /// \short Create an "bounding" element (here actually a 1D point element
160  /// of type SpinePointFluidInterfaceBoundingElement<ELEMENT> that allows
161  /// the application of a contact angle boundary condition on the
162  /// the specified face.
164  const int &face_index)
165  {
166  //Create a temporary pointer to the appropriate FaceElement
169 
170  //Attach the geometrical information to the new element
171  this->build_face_element(face_index,face_el_pt);
172 
173  //Set the index at which the unknowns are stored from the element
174  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
175 
176  //Set the value of the nbulk_value, the node is not resized
177  //in this problem, so it will just be the actual nvalue
178  face_el_pt->nbulk_value(0) = face_el_pt->node_pt(0)->nvalue();
179 
180  //Set of unique geometric data that is used to update the bulk,
181  //but is not used to update the face
182  std::set<Data*> unique_additional_geom_data;
183 
184  //Get all the geometric data for this (bulk) element
185  this->assemble_set_of_all_geometric_data(unique_additional_geom_data);
186 
187  //Now assemble the set of geometric data for the face element
188  std::set<Data*> unique_face_geom_data_pt;
189  face_el_pt->assemble_set_of_all_geometric_data(unique_face_geom_data_pt);
190 
191  //Erase the face geometric data from the additional data
192  for(std::set<Data*>::iterator it=unique_face_geom_data_pt.begin();
193  it!=unique_face_geom_data_pt.end();++it)
194  {unique_additional_geom_data.erase(*it);}
195 
196  //Finally add all unique additional data as external data
197  for(std::set<Data*>::iterator it = unique_additional_geom_data.begin();
198  it!= unique_additional_geom_data.end();++it)
199  {
200  face_el_pt->add_external_data(*it);
201  }
202 
203  //Return the value of the pointer
204  return face_el_pt;
205  }
206 
207 };
208 
209 
210 /////////////////////////////////////////////////////////////////////////
211 /////////////////////////////////////////////////////////////////////////
212 /////////////////////////////////////////////////////////////////////////
213 
214 
215 
216 //=======================================================================
217 /// One-dimensional interface elements that are used when the mesh
218 /// deformation is handled by a set of equations that modify the nodal
219 /// positions. These elements are FaceElements attached to 2D bulk fluid
220 /// elements and the fluid element is passed as a template parameter to
221 /// the element.
222 /// Optionally an external pressure may be specified, which must be
223 /// passed to the element as external data. The default value of the external
224 /// pressure is zero.
225 //=======================================================================
226 template<class ELEMENT>
228 public virtual Hijacked<FaceGeometry<ELEMENT> >,
229 public LineFluidInterfaceElement
230 {
231  private:
232 
233  /// \short ID of the Lagrange Lagrange multiplier (in the collection of nodal
234  /// values accomodated by resizing)
235  unsigned Id;
236 
237  /// \short Equation number of the kinematic BC associated with node j.
238  /// (This is the equation for the Lagrange multiplier)
239  int kinematic_local_eqn(const unsigned &j)
240  {
241  // Get the index of the nodal value associated with Lagrange multiplier
242  const unsigned lagr_index=
243  dynamic_cast<BoundaryNodeBase*>(node_pt(j))->
244  index_of_first_value_assigned_by_face_element(Id);
245 
246  // Return nodal value
247  return this->nodal_local_eqn(j,lagr_index);
248  }
249 
250 
251  /// \short Hijacking the kinematic condition corresponds to hijacking the
252  /// variables associated with the Lagrange multipliers that are assigned
253  /// on construction of this element.
254  void hijack_kinematic_conditions(const Vector<unsigned> &bulk_node_number)
255  {
256  //Loop over all the nodes that are passed in
257  for(Vector<unsigned>::const_iterator it=bulk_node_number.begin();
258  it!=bulk_node_number.end();++it)
259  {
260  //Get the index associated with the Id for each node
261  //(the Lagrange multiplier)
262  unsigned n_lagr = dynamic_cast<BoundaryNodeBase*>(node_pt(*it))->
263  index_of_first_value_assigned_by_face_element(Id);
264 
265  //Hijack the appropriate value and delete the returned Node
266  delete this->hijack_nodal_value(*it,n_lagr);
267  }
268  }
269 
270  public:
271 
272  /// \short The "global" intrinsic coordinate of the element when
273  /// viewed as part of a geometric object should be given by
274  /// the FaceElement representation, by default
275  double zeta_nodal(const unsigned &n, const unsigned &k,
276  const unsigned &i) const
277  {return FaceElement::zeta_nodal(n,k,i);}
278 
279 
280  /// \short Constructor, pass a pointer to the bulk element and the face index
281  /// of the bulk element to which the element is to be attached to.
282  /// The optional identifier can be used
283  /// to distinguish the additional nodal value (Lagr mult) created by
284  /// this element from those created by other FaceElements.
285  ElasticLineFluidInterfaceElement(FiniteElement* const &element_pt,
286  const int &face_index,
287  const unsigned &id=0) :
288  // Final optional argument to specify id for lagrange multiplier
289  FaceGeometry<ELEMENT>(), LineFluidInterfaceElement(), Id(id)
290  {
291  //Attach the geometrical information to the element
292  //This function also assigned nbulk_value from required_nvalue of the
293  //bulk element
294  element_pt->build_face_element(face_index,this);
295 
296  //Find the index at which the velocity unknowns are stored
297  //from the bulk element
298  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
299  this->U_index_interface.resize(2);
300  for(unsigned i=0;i<2;i++)
301  {
302  this->U_index_interface[i] = cast_element_pt->u_index_nst(i);
303  }
304 
305  //Read out the number of nodes on the face
306  unsigned n_node_face = this->nnode();
307 
308  //Set the additional data values in the face
309  //There is one additional values at each node --- the lagrange multiplier
310  Vector<unsigned> additional_data_values(n_node_face);
311  for(unsigned i=0;i<n_node_face;i++)
312  {
313  additional_data_values[i] = 1;
314  }
315 
316  // Now add storage for Lagrange multipliers and set the map containing
317  // the position of the first entry of this face element's
318  // additional values.
319  add_additional_values(additional_data_values,id);
320 
321  } //End of constructor
322 
323  /// Return the Lagrange multiplier at local node j
324  double &lagrange(const unsigned &j)
325  {
326  // Get the index of the nodal value associated with Lagrange multiplier
327  unsigned lagr_index=dynamic_cast<BoundaryNodeBase*>(node_pt(j))->
328  index_of_first_value_assigned_by_face_element(Id);
329  return *node_pt(j)->value_pt(lagr_index);
330  }
331 
332  /// Fill in contribution to residuals and Jacobian
333  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
334  DenseMatrix<double> &jacobian)
335  {
336  //Call the generic routine with the flag set to 1
338 
339  //Call the generic finite difference routine for the solid variables
340  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
341  }
342 
343  /// Overload the output function
344  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
345 
346  /// Output the element
347  void output(std::ostream &outfile, const unsigned &n_plot)
348  {LineFluidInterfaceElement::output(outfile,n_plot);}
349 
350  ///Overload the C-style output function
351  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
352 
353  ///C-style Output function
354  void output(FILE* file_pt, const unsigned &n_plot)
355  {LineFluidInterfaceElement::output(file_pt,n_plot);}
356 
357 
358  /// \short Helper function to calculate the additional contributions
359  /// to be added at each integration point. This deals with
360  /// Lagrange multiplier contribution
362  Vector<double> &residuals,
363  DenseMatrix<double> &jacobian,
364  const unsigned &flag,
365  const Shape &psif,
366  const DShape &dpsifds,
367  const Vector<double> &interpolated_x,
368  const Vector<double> &interpolated_n,
369  const double &W,
370  const double &J)
371  {
372  //Loop over the shape functions
373  unsigned n_node = this->nnode();
374 
375  // Assemble Lagrange multiplier in loop over the shape functions
376  double interpolated_lagrange = 0.0;
377  for(unsigned l=0;l<n_node;l++)
378  {
379  //Note same shape functions used for lagrange multiplier field
380  interpolated_lagrange += lagrange(l)*psif[l];
381  }
382 
383  int local_eqn=0, local_unknown = 0;
384 
385  //Loop over the shape functions to assemble contributions
386  for(unsigned l=0;l<n_node;l++)
387  {
388  //Loop over the directions
389  for(unsigned i=0;i<2;i++)
390  {
391  //Now using the same shape functions for the elastic equations,
392  //so we can stay in the loop
393  local_eqn = this->position_local_eqn(l,0,i);
394  if(local_eqn >= 0)
395  {
396  //Add in a "Lagrange multiplier"
397  residuals[local_eqn] -=
398  interpolated_lagrange*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  //Derivatives w.r.t. solid positions will be handled by FDing
407  //That leaves the "lagrange multipliers" only
408  local_unknown = kinematic_local_eqn(l2);
409  if(local_unknown >= 0)
410  {
411  jacobian(local_eqn,local_unknown) -=
412  psif[l2]*interpolated_n[i]*psif[l]*W*J;
413  }
414  }
415  } //End of Jacobian calculation
416  }
417  }
418  } //End of loop over shape functions
419  }
420 
421 
422  /// \short Create an "bounding" element (here actually a 1D point element
423  /// of type ElasticPointFluidInterfaceBoundingElement<ELEMENT> that allows
424  /// the application of a contact angle boundary condition on the
425  /// the specified face.
427  const int &face_index)
428  {
429  //Create a temporary pointer to the appropriate FaceElement
432 
433  //Attach the geometrical information to the new element
434  this->build_face_element(face_index,face_el_pt);
435 
436  //Set the index at which the unknowns are stored from the element
437  face_el_pt->u_index_interface_boundary() = this->U_index_interface;
438 
439  //Set the value of the nbulk_value, the node is not resized
440  //in this problem, so it will just be the actual nvalue
441  face_el_pt->nbulk_value(0) = this->node_pt(0)->nvalue();
442 
443  //Pass the ID down
444  face_el_pt->set_id(Id);
445 
446  //Find the nodes
447  std::set<SolidNode*> set_of_solid_nodes;
448  unsigned n_node = this->nnode();
449  for(unsigned n=0;n<n_node;n++)
450  {
451  set_of_solid_nodes.insert(static_cast<SolidNode*>(this->node_pt(n)));
452  }
453 
454  //Delete the nodes from the face
455  n_node = face_el_pt->nnode();
456  for(unsigned n=0;n<n_node;n++)
457  {
458  set_of_solid_nodes.erase(static_cast<SolidNode*>(face_el_pt->node_pt(n)));
459  }
460 
461  //Now add these as external data
462  for(std::set<SolidNode*>::iterator it=set_of_solid_nodes.begin();
463  it!=set_of_solid_nodes.end();++it)
464  {
465  face_el_pt->add_external_data((*it)->variable_position_pt());
466  }
467 
468  //Return the value of the pointer
469  return face_el_pt;
470  }
471 
472 };
473 
474 }
475 
476 #endif
477 
478 
479 
480 
481 
482 
void output(std::ostream &outfile)
Overload the output function.
Specialise the elastic update template class to concrete 1D case.
int kinematic_local_eqn(const unsigned &n)
In spine elements, the kinematic condition is the equation used to determine the unknown spine height...
unsigned Id
ID of the Lagrange Lagrange multiplier (in the collection of nodal values accomodated by resizing) ...
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(std::ostream &outfile, const unsigned &n_plot)
Output the element.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt)
Overload the C-style output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contribution to the residuals and the jacobian.
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
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
SpineLineFluidInterfaceElement(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...
ElasticLineFluidInterfaceElement(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 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 hijack_kinematic_conditions(const Vector< unsigned > &bulk_node_number)
Hijacking the kinematic condition corresponds to hijacking the variables associated with the spine he...
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, const unsigned &n_plot)
C-style Output function.
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...
Vector< unsigned > & u_index_interface_boundary()
Access for nodal index at which the velocity components are stored.
double & lagrange(const unsigned &j)
Return the Lagrange multiplier at local node j.
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
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.
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 1D point element of type SpinePointFluidInterfaceBoundi...
virtual FluidInterfaceBoundingElement * make_bounding_element(const int &face_index)
Create an "bounding" element (here actually a 1D point element of type ElasticPointFluidInterfaceBoun...
Spine version of the PointFluidInterfaceBoundingElement.