impose_parallel_outflow_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: 1275 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-12-30 06:46:51 +0000 (Fri, 30 Dec 2016) $
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_PARALL_ELEMENTS_HEADER
31 #define OOMPH_IMPOSE_PARALL_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 ImposeParallelOutflowElement
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 parallel outflow and
45 /// impose the pressure.
46 //========================================================================
47  template <class ELEMENT>
49  public virtual FaceGeometry<ELEMENT>,
50  public virtual FaceElement
51  {
52  private :
53 
54  /// pointer to imposed pressure -- if null then no pressure imposed.
55  double* Pressure_pt;
56 
57  /// Lagrange Id
58  unsigned Id;
59 
60 
61  public:
62 
63  /// \short Constructor takes a "bulk" element, the
64  /// index that identifies which face the
65  /// ImposeParallelOutflowElement is supposed
66  /// to be attached to, and the face element ID
68  (FiniteElement* const &element_pt,
69  const int &face_index,const unsigned &id=0) :
71  {
72  // set the Id
73  Id=id;
74 
75  //Build the face element
76  element_pt->build_face_element(face_index,this);
77 
78  // dimension of the bulk element
79  unsigned dim=element_pt->dim();
80 
81  // we need dim-1 additional values for each FaceElement node
82  Vector<unsigned> n_additional_values(nnode(), dim-1);
83 
84  // add storage for lagrange multipliers and set the map containing
85  // the position of the first entry of this face element's
86  // additional values.
87  add_additional_values(n_additional_values,id);
88 
89  // set the pressure pointer to zero
90  Pressure_pt=0;
91  }
92 
93  /// Fill in the residuals
95  {
96  //Call the generic routine with the flag set to 0
99  }
100 
101  /// Fill in contribution from Jacobian
103  DenseMatrix<double> &jacobian)
104  {
105  //Call the generic routine with the flag set to 1
107  residuals,jacobian,1);
108  }
109 
110  ///Overload the output function
111  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
112 
113  ///Output function: x,y,[z],u,v,[w],p in tecplot format
114  void output(std::ostream &outfile, const unsigned &nplot)
115  {FiniteElement::output(outfile,nplot);}
116 
117  /// \short The "global" intrinsic coordinate of the element when
118  /// viewed as part of a geometric object should be given by
119  /// the FaceElement representation, by default
120  /// This final over-ride is required because both SolidFiniteElements
121  /// and FaceElements overload zeta_nodal
122  double zeta_nodal(const unsigned &n, const unsigned &k,
123  const unsigned &i) const
124  {return FaceElement::zeta_nodal(n,k,i);}
125 
126 
127  /// Access function for the pressure
128  double* &pressure_pt() {return Pressure_pt;}
129 
130  protected:
131 
132  /// \short Helper function to compute the residuals and, if flag==1, the
133  /// Jacobian
135  Vector<double> &residuals,
136  DenseMatrix<double> &jacobian,
137  const unsigned& flag)
138  {
139  //Find out how many nodes there are
140  unsigned n_node = nnode();
141 
142  // Dimension of element
143  unsigned dim_el=dim();
144 
145  //Set up memory for the shape functions
146  Shape psi(n_node);
147 
148  //Set the value of n_intpt
149  unsigned n_intpt = integral_pt()->nweight();
150 
151  //to store local equation number
152  int local_eqn=0;
153  int local_unknown=0;
154 
155  //to store normal vector
156  Vector<double> norm_vec(dim_el+1);
157 
158  //to store tangantial vectors
159  Vector<Vector<double> > tang_vec(dim_el,Vector<double>(dim_el+1));
160 
161  //get the value at which the velocities are stored
162  Vector<unsigned> u_index(dim_el+1);
163  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
164  for(unsigned i=0;i<dim_el+1;i++) {u_index[i] = el_pt->u_index_nst(i);}
165 
166  //Loop over the integration points
167  for(unsigned ipt=0;ipt<n_intpt;ipt++)
168  {
169  //Get the integral weight
170  double w = integral_pt()->weight(ipt);
171 
172  //Jacobian of mapping
173  double J=J_eulerian_at_knot(ipt);
174 
175  //Premultiply the weights and the Jacobian
176  double W = w*J;
177 
178  //Calculate the shape functions
179  shape_at_knot(ipt,psi);
180 
181  //compute the velocity and the Lagrange multiplier
182  Vector<double> interpolated_u(dim_el+1,0.0);
183  Vector<double> lambda(dim_el,0.0);
184  // Loop over nodes
185  for(unsigned j=0;j<n_node;j++)
186  {
187  //Assemble the velocity component
188  for(unsigned i=0;i<dim_el+1;i++)
189  {
190  interpolated_u[i] += nodal_value(j,u_index[i])*psi(j);
191  }
192 
193  // Cast to a boundary node
194  BoundaryNodeBase *bnod_pt =
195  dynamic_cast<BoundaryNodeBase*>(node_pt(j));
196 
197  // get the node
198  Node* nod_pt = node_pt(j);
199 
200  // Get the index of the first nodal value associated with
201  // this FaceElement
202  unsigned first_index=
204 
205  //Assemble the Lagrange multiplier
206  for(unsigned l=0;l<dim_el;l++)
207  {
208  lambda[l]+=nod_pt->value(first_index+l) * psi(j);
209  }
210  }
211 
212  this->continuous_tangent_and_outer_unit_normal(ipt,tang_vec,norm_vec);
213 
214  // Assemble residuals and jacobian
215 
216  //Loop over the nodes
217  for(unsigned j=0;j<n_node;j++)
218  {
219  // Cast to a boundary node
220  BoundaryNodeBase *bnod_pt =
221  dynamic_cast<BoundaryNodeBase*>(node_pt(j));
222 
223  // Get the index of the first nodal value associated with
224  // this FaceElement
225  unsigned first_index=
227 
228  //loop over the lagrange multiplier components
229  for(unsigned l=0;l<dim_el;l++)
230  {
231  // Local eqn number for the l-th component of lamdba
232  //in the j-th element
233  local_eqn=nodal_local_eqn(j,first_index+l);
234 
235  if (local_eqn>=0)
236  {
237  for(unsigned i=0; i<dim_el+1; i++)
238  {
239  // Assemble residual for lagrange multiplier
240  residuals[local_eqn]+=
241  interpolated_u[i] * tang_vec[l][i] * psi(j)* W;
242 
243  // Assemble Jacobian for Lagrange multiplier:
244  if (flag==1)
245  {
246  // Loop over the nodes again for unknowns
247  for(unsigned jj=0;jj<n_node;jj++)
248  {
249  // Local eqn number for the i-th component
250  //of the velocity in the jj-th element
251  local_unknown=nodal_local_eqn(jj,u_index[i]);
252  if (local_unknown>=0)
253  {
254  jacobian(local_eqn,local_unknown)+=
255  tang_vec[l][i]*psi(jj)*psi(j)*W;
256  }
257  }
258  }
259  }
260  }
261  }
262 
263  //Loop over the directions
264  for(unsigned i=0;i<dim_el+1;i++)
265  {
266  // Local eqn number for the i-th component of the
267  //velocity in the j-th element
268  local_eqn = nodal_local_eqn(j,u_index[i]);
269 
270  if (local_eqn>=0)
271  {
272  // Add the contribution of the imposed pressure
273  if (Pressure_pt!=0)
274  {
275  residuals[local_eqn]-= (*Pressure_pt)* norm_vec[i]*psi(j)*W;
276  }
277 
278  // Add lagrange multiplier contribution to the bulk equation
279  for(unsigned l=0;l<dim_el;l++)
280  {
281  // Add to residual
282  residuals[local_eqn]+=tang_vec[l][i]*psi(j)*lambda[l]*W;
283 
284  // Do Jacobian too?
285  if (flag==1)
286  {
287  // Loop over the nodes again for unknowns
288  for(unsigned jj=0;jj<n_node;jj++)
289  {
290  // Cast to a boundary node
291  BoundaryNodeBase *bnode_pt =
292  dynamic_cast<BoundaryNodeBase*>(node_pt(jj));
293 
294  // Local eqn number for the l-th component of lamdba
295  // in the jj-th element
296  local_unknown=nodal_local_eqn
297  (jj,bnode_pt->
298  index_of_first_value_assigned_by_face_element(Id)+l);
299  if (local_unknown>=0)
300  {
301  jacobian(local_eqn,local_unknown)+=
302  tang_vec[l][i]*psi(jj)*psi(j)*W;
303  }
304  }
305  }
306  }
307  }
308  }
309  }
310  }
311  }
312 
313  /// \short The number of degrees of freedom types that this element is
314  /// divided into: The ImposeParallelOutflowElement is responsible for 1(2) of
315  /// it's own degrees of freedom. In addition to this, it also classifies the
316  /// 2(3) constrained velocity bulk degrees of freedom which this face element
317  /// is attached to.
318  unsigned ndof_types() const
319  {
320  // The the dof types in this element is 1 in 2D, 2 in 3D, so we use the
321  // elemental dimension function this->dim().
322  // In addition, we have the number of constrained velocity dof types, this
323  // is this->dim()+1 dof types.
324  // In total we have 2*(this->dim()) + 1 dof types.
325  return 2 * (this->dim()) + 1;
326  }
327 
328  /// \short Create a list of pairs for all unknowns in this element,
329  /// so that the first entry in each pair contains the global equation
330  /// number of the unknown, while the second one contains the LOCAL number
331  /// of the "DOF types" that this unknown is associated with.
332  /// (Function can obviously only be called if the equation numbering
333  /// scheme has been set up.)
334  ///
335  /// For the ImposeParallelOutflowElement, we need to classify both the
336  /// constrained velocity from the bulk element and the dof types from this
337  /// element. In 3D, let up, vp, wp be the constrained velocity dof types
338  /// associated with the ImposeParallelOutflowElement. Let Lp1 and Lp2 be the
339  /// two Lagrange multiplier dof types associated with the
340  /// ImposeParallelOutflowElement. Then we have:
341  ///
342  /// 0 1 2 3 4
343  /// up vp wp L1 L2
345  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
346  {
347 
348  // Temporary pair (used to store dof lookup prior to being added to list)
349  std::pair<unsigned,unsigned> dof_lookup;
350 
351  // Number of nodes in this element.
352  const unsigned n_node = this->nnode();
353 
354  // The elemental dimension.
355  const unsigned dim_el = this->dim();
356 
357  // The spatial dimension, since this is a face element, as assume that the
358  // spatial dimension is one higher than the elemental dimension.
359  const unsigned dim_bulk = dim_el + 1;
360 
361  // The classification of the constrained bulk velocity dof types is
362  // i, i = 0, ..., dim_bulk.
363  //
364  // The classification of the lagrange multiplier dof types is
365  // i + dim_bulk, i = 0, ..., dim_el
366  // as it is offset by the spatial dimension.
367 
368  //Loop over the nodes
369  for(unsigned node_i =0; node_i < n_node; node_i++)
370  {
371  // First we classify the constrained velocity dof types from the bulk
372  // element. NOTE: This is a RE-CLASSIFICATION, since the dof types we are
373  // classifying should already be classified at least once by the bulk
374  // element.
375  // Furthermore, we assume that the velocity dof types comes first,
376  // then the pressure dof types. This is indeed true with OOMPH-LIB's
377  // implementation of Navier-Stokes elements.
378  for(unsigned velocity_i = 0; velocity_i < dim_bulk; velocity_i++)
379  {
380  // We only concern ourselves with the nodes of the "bulk" element
381  // that are associated with this "face" element. Bulk_node_number gives
382  // us the mapping from the face element's node to the bulk element's
383  // node.
384  // The local equation number is required to check if the value is pinned,
385  // if it is not pinned, the local equation number is required to get the
386  // global equation number.
387  int local_eqn = Bulk_element_pt
389  velocity_i);
390 
391  // Ignore pinned values.
392  if(local_eqn >= 0)
393  {
394  // Store the dof lookup in temporary pair: First entry in pair
395  // is the global equation number; second entry is the LOCAL dof type.
396  // It is assumed that the velocity dof types comes first. We do not
397  // re-classify the pressure dof types as they are not constrained.
398  dof_lookup.first = Bulk_element_pt->eqn_number(local_eqn);
399  dof_lookup.second = velocity_i;
400  dof_lookup_list.push_front(dof_lookup);
401  } // ignore pinned nodes "if(local-eqn>=0)"
402  } // for loop over the velocity components
403 
404  // Now we classify the dof types of this face element.
405  // Cast to a boundary node
406  BoundaryNodeBase *bnod_pt =
407  dynamic_cast<BoundaryNodeBase*>(node_pt(node_i));
408 
409  // Loop over directions in this Face(!)Element
410  for(unsigned dim_i=0;dim_i<dim_el;dim_i++)
411  {
412  // Local eqn number:
413  int local_eqn
415  (node_i,
417 
418  if (local_eqn>=0)
419  {
420  // Store dof lookup in temporary pair: First entry in pair
421  // is global equation number; second entry is dof type.
422  // We assume that the first 2(3) types are fluid dof types in
423  // 2(3) spatial dimensions (classified above):
424  //
425  // 0 1 2 3 4
426  // up vp wp L1 L2
427  //
428  // Thus the offset is the dim_bulk.
429  dof_lookup.first = this->eqn_number(local_eqn);
430  dof_lookup.second = dim_i + dim_bulk;
431 
432  // Add to list
433  dof_lookup_list.push_front(dof_lookup);
434  } // if local_eqn > 0
435  } // for: loop over the directions in the face element.
436  } // for: loop over the nodes in the face element.
437 
438  } // get_dof_numbers_for_unknowns
439 
440  }; // End of ImposeParallelOutflowElement class
441 
442 }
443 
444 #endif
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
void continuous_tangent_and_outer_unit_normal(const Vector< double > &s, Vector< Vector< double > > &tang_vec, Vector< double > &unit_normal) const
Compute the tangent vector(s) and the outer unit normal vector at the specified local coordinate...
Definition: elements.cc:5211
double * Pressure_pt
pointer to imposed pressure – if null then no pressure imposed.
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
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...
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
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
ImposeParallelOutflowElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor takes a "bulk" element, the index that identifies which face the ImposeParallelOutflowEle...
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
ImposeParallelOutflowElement are elements that coincide with the faces of higher-dimensional "bulk" e...
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
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.
void output(std::ostream &outfile)
Overload the output function.
FaceElement()
Constructor: Initialise all appropriate member data.
Definition: elements.h:4201
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 ...
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
unsigned ndof_types() const
The number of degrees of freedom types that this element is divided into: The ImposeParallelOutflowEl...
double *& pressure_pt()
Access function for the pressure.
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
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from 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