impose_impenetrability_element.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 #ifndef OOMPH_IMPOSE_IMPENETRABLE_ELEMENTS_HEADER
31 #define OOMPH_IMPOSE_IMPENETRABLE_ELEMENTS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35 #include <oomph-lib-config.h>
36 #endif
37 
38 namespace oomph
39 {
40 //========================================================================
41 /// \short ImposeImpenetrabilityElement
42 /// are elements that coincide with the faces of
43 /// higher-dimensional "bulk" elements. They are used on
44 /// boundaries where we would like to impose impenetrability.
45 //========================================================================
46  template <class ELEMENT>
48  public virtual FaceGeometry<ELEMENT>,
49  public virtual FaceElement
50  {
51  private :
52 
53  /// Lagrange Id
54  unsigned Id;
55 
56  public:
57 
58  /// \short Constructor takes a "bulk" element, the
59  /// index that identifies which face the
60  /// ImposeImpenetrabilityElement is supposed
61  /// to be attached to, and the face element ID
63  (FiniteElement* const &element_pt,
64  const int &face_index,const unsigned &id=0) :
66  {
67  // set the Id
68  Id=id;
69 
70  //Build the face element
71  element_pt->build_face_element(face_index,this);
72 
73  // we need 1 additional values for each FaceElement node
74  Vector<unsigned> n_additional_values(nnode(), 1);
75 
76  // add storage for lagrange multipliers and set the map containing
77  // the position of the first entry of this face element's
78  // additional values.
79  add_additional_values(n_additional_values,id);
80 
81  }
82 
83 
84  /// Fill in the residuals
86  {
87  //Call the generic routine with the flag set to 0
90  }
91 
92  /// Fill in contribution from Jacobian
94  DenseMatrix<double> &jacobian)
95  {
96  //Call the generic routine with the flag set to 1
98  residuals,jacobian,1);
99  }
100 
101  ///Overload the output function
102  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
103 
104  ///Output function: x,y,[z],u,v,[w],p in tecplot format
105  void output(std::ostream &outfile, const unsigned &nplot)
106  {FiniteElement::output(outfile,nplot);}
107 
108  /// \short The "global" intrinsic coordinate of the element when
109  /// viewed as part of a geometric object should be given by
110  /// the FaceElement representation, by default
111  /// This final over-ride is required because both SolidFiniteElements
112  /// and FaceElements overload zeta_nodal
113  double zeta_nodal(const unsigned &n, const unsigned &k,
114  const unsigned &i) const
115  {return FaceElement::zeta_nodal(n,k,i);}
116 
117  protected:
118 
119  /// \short Helper function to compute the residuals and, if flag==1, the
120  /// Jacobian
122  Vector<double> &residuals,
123  DenseMatrix<double> &jacobian,
124  const unsigned& flag)
125  {
126  //Find out how many nodes there are
127  unsigned n_node = nnode();
128 
129  // Dimension of element
130  unsigned dim_el=dim();
131 
132  //Set up memory for the shape functions
133  Shape psi(n_node);
134 
135  //Set the value of n_intpt
136  unsigned n_intpt = integral_pt()->nweight();
137 
138  //to store local equation number
139  int local_eqn=0;
140  int local_unknown=0;
141 
142  //to store normal vector
143  Vector<double> norm_vec(dim_el+1);
144 
145  //get the value at which the velocities are stored
146  Vector<unsigned> u_index(dim_el+1);
147  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
148  for(unsigned i=0;i<dim_el+1;i++) {u_index[i] = el_pt->u_index_nst(i);}
149 
150  //Loop over the integration points
151  for(unsigned ipt=0;ipt<n_intpt;ipt++)
152  {
153 
154  //Get the integral weight
155  double w = integral_pt()->weight(ipt);
156 
157  //Jacobian of mapping
158  double J=J_eulerian_at_knot(ipt);
159 
160  //Premultiply the weights and the Jacobian
161  double W = w*J;
162 
163  //Calculate the shape functions
164  shape_at_knot(ipt,psi);
165 
166  //compute the velocity and the Lagrange multiplier
167  Vector<double> interpolated_u(dim_el+1,0.0);
168  double lambda=0.0;
169 
170  // Loop over nodes
171  for(unsigned j=0;j<n_node;j++)
172  {
173  //Assemble the velocity component
174  for(unsigned i=0;i<dim_el+1;i++)
175  {
176  interpolated_u[i] += nodal_value(j,u_index[i])*psi(j);
177  }
178 
179  // Cast to a boundary node
180  BoundaryNodeBase *bnod_pt =
181  dynamic_cast<BoundaryNodeBase*>(node_pt(j));
182 
183  // get the node
184  Node* nod_pt = node_pt(j);
185 
186  // Get the index of the first nodal value associated with
187  // this FaceElement
188  unsigned first_index=
190 
191  //Assemble the Lagrange multiplier
192  lambda+=nod_pt->value(first_index) * psi(j);
193  }
194 
195  // compute the normal vector
196  outer_unit_normal(ipt,norm_vec);
197 
198  // Assemble residuals and jacobian
199 
200  //Loop over the nodes
201  for(unsigned j=0;j<n_node;j++)
202  {
203  // Cast to a boundary node
204  BoundaryNodeBase *bnod_pt =
205  dynamic_cast<BoundaryNodeBase*>(node_pt(j));
206 
207  // Get the index of the first nodal value associated with
208  // this FaceElement
209  unsigned first_index=
211 
212  // Local eqn number for the l-th component of lamdba
213  //in the j-th element
214  local_eqn=nodal_local_eqn(j,first_index);
215 
216  if (local_eqn>=0)
217  {
218  for(unsigned i=0; i<dim_el+1; i++)
219  {
220  // Assemble residual for lagrange multiplier
221  residuals[local_eqn]+=
222  interpolated_u[i] * norm_vec[i] * psi(j)* W;
223 
224  // Assemble Jacobian for Lagrange multiplier:
225  if (flag==1)
226  {
227  // Loop over the nodes again for unknowns
228  for(unsigned jj=0;jj<n_node;jj++)
229  {
230  // Local eqn number for the i-th component
231  //of the velocity in the jj-th element
232  local_unknown=nodal_local_eqn(jj,u_index[i]);
233  if (local_unknown>=0)
234  {
235  jacobian(local_eqn,local_unknown)+=
236  norm_vec[i]*psi(jj)*psi(j)*W;
237  }
238  }
239  }
240  }
241  }
242 
243  //Loop over the directions
244  for(unsigned i=0;i<dim_el+1;i++)
245  {
246  // Local eqn number for the i-th component of the
247  //velocity in the j-th element
248  local_eqn = nodal_local_eqn(j,u_index[i]);
249 
250  if (local_eqn>=0)
251  {
252  // Add lagrange multiplier contribution to the bulk equation
253  // Add to residual
254  residuals[local_eqn]+=norm_vec[i]*psi(j)*lambda*W;
255 
256  // Do Jacobian too?
257  if (flag==1)
258  {
259  // Loop over the nodes again for unknowns
260  for(unsigned jj=0;jj<n_node;jj++)
261  {
262  // Cast to a boundary node
263  BoundaryNodeBase *bnode_pt =
264  dynamic_cast<BoundaryNodeBase*>(node_pt(jj));
265 
266  // Local eqn number for the l-th component of lamdba
267  // in the jj-th element
268  local_unknown=nodal_local_eqn
269  (jj,bnode_pt->
270  index_of_first_value_assigned_by_face_element(Id));
271  if (local_unknown>=0)
272  {
273  jacobian(local_eqn,local_unknown)+=
274  norm_vec[i]*psi(jj)*psi(j)*W;
275  }
276  }
277  }
278  }
279  }
280  }
281  }
282  }
283 
284 
285  /// \short The number of "DOF types" that degrees of freedom in this element
286  /// are sub-divided into: Just the solid degrees of freedom themselves.
287  unsigned ndof_types() const
288  {
289  // There is only ever one normal. Plus the constrained velocities.
290  //unsigned ndofndof = 1 + additional_ndof_types();
291  //std::cout << "ndof: " << ndofndof << std::endl;
292 
293  return (1 + additional_ndof_types());
294  }
295 
296 
297  unsigned additional_ndof_types() const
298  {
299  // Additional dof types for the constained bulk velocities
300  // two velocities for a 2D problem, 3 for 3D.
301  return (this->dim() + 1);
302  }
303 
304  /// \short Create a list of pairs for all unknowns in this element,
305  /// so that the first entry in each pair contains the global equation
306  /// number of the unknown, while the second one contains the number
307  /// of the "DOF type" that this unknown is associated with.
308  /// (Function can obviously only be called if the equation numbering
309  /// scheme has been set up.)
311  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
312  {
313  // temporary pair (used to store dof lookup prior to
314  // being added to list)
315  std::pair<unsigned,unsigned> dof_lookup;
316 
317  // number of nodes
318  const unsigned n_node = this->nnode();
319 
320  //Loop over directions in this Face(!)Element
321  unsigned dim_el = this->dim();
322  //for(unsigned i=0;i<dim_el;i++)
323  {
324  unsigned i=0;
325  //Loop over the nodes
326  for(unsigned j=0;j<n_node;j++)
327  {
328  // Cast to a boundary node
329  BoundaryNodeBase *bnod_pt =
330  dynamic_cast<BoundaryNodeBase*>(node_pt(j));
331 
332  // Local eqn number:
333  int local_eqn=nodal_local_eqn
335  if (local_eqn>=0)
336  {
337  // store dof lookup in temporary pair: First entry in pair
338  // is global equation number; second entry is dof type
339  dof_lookup.first = this->eqn_number(local_eqn);
340  dof_lookup.second = i+additional_ndof_types();
341 
342  // add to list
343  dof_lookup_list.push_front(dof_lookup);
344  }
345  }
346  }
347 
348  //*
349  // Now we do the bulk elements. Each velocity component of a constrained dof
350  // of a different type of FaceElement has a different dof_type. E.g. Consider
351  // the Navier Stokes equations in three spatial dimensions with parallel
352  // outflow (using ImposeParallelOutflowElement with Boundary_id = 1) and
353  // tangential flow (using ImposeTangentialFlowElement with Boundary_id = 2)
354  // imposed along two different boundaries.
355  // There will be 10 dof types:
356  // 0 1 2 3 4 5 6 7 8 9
357  // u v w p u1 v1 w1 u2 v2 w2
358 
359  // Loop over only the nodes of the "bulk" element that are associated
360  // with this "face" element.
361  //cout << "n_node: " << n_node << endl;
362  unsigned const bulk_dim = dim_el + 1;
363  //cout << "bulk_dim: " << bulk_dim << endl;
364  for(unsigned node_i = 0; node_i < n_node; node_i++)
365  {
366  // Loop over the velocity components
367  for(unsigned velocity_i = 0; velocity_i < bulk_dim; velocity_i++)
368  {
369  // Calculating the offset for this Boundary_id.
370  // 0 1 2 3 4 5 6 7 8 9
371  // u v w p u1 v1 w1 u2 v2 w2
372  //
373  // for the first surface mesh, offset = 4
374  // for the second surface mesh, offset = 7
375  //unsigned offset = bulk_dim * Boundary_id + 1;
376 
377  // The local equation number is required to check if the value is pinned,
378  // if it is not pinned, the local equation number is required to get the
379  // global equation number.
380  int local_eqn = Bulk_element_pt
382  velocity_i);
383 
384  // ignore pinned values
385  if(local_eqn >= 0)
386  {
387  // store the dof lookup in temporary pair: First entry in pair
388  // is the global equation number; second entry is the dof type
389  dof_lookup.first = Bulk_element_pt->eqn_number(local_eqn);
390  dof_lookup.second = velocity_i;
391  dof_lookup_list.push_front(dof_lookup);
392 
393  //RRRcout << "Face v: " << dof_lookup.first
394  //RRR << ", doftype: " << dof_lookup.second << endl;
395 
396  } // ignore pinned nodes "if(local-eqn>=0)"
397  } // for loop over the velocity components
398  } // for loop over bulk nodes only
399 
400  } // get unknowns...
401  };
402 }
403 #endif
ImposeImpenetrabilityElement are elements that coincide with the faces of higher-dimensional "bulk" e...
void output(std::ostream &outfile)
Overload the output function.
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
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
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
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Return the index of the first value associated with the i-th face element value. If no argument is sp...
Definition: nodes.h:1924
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Just the soli...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:5695
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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 fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
Vector< unsigned > Bulk_node_number
List of indices of the local node numbers in the "bulk" element that correspond to the local node num...
Definition: elements.h:4156
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
FaceElement()
Constructor: Initialise all appropriate member data.
Definition: elements.h:4201
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
FiniteElement * Bulk_element_pt
Pointer to the associated higher-dimensional "bulk" element.
Definition: elements.h:4152
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
ImposeImpenetrabilityElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor takes a "bulk" element, the index that identifies which face the ImposeImpenetrabilityEle...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
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
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_generic_contribution_to_residuals_parall_lagr_multiplier(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to compute the residuals and, if flag==1, the Jacobian.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228