fluid_traction_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 elements that are used to integrate fluid tractions
31 //This includes the guts (i.e. equations) because we want to inline them
32 //for faster operation, although it slows down the compilation!
33 
34 #ifndef OOMPH_FLUID_TRACTION_ELEMENTS_HEADER
35 #define OOMPH_FLUID_TRACTION_ELEMENTS_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39 #include <oomph-lib-config.h>
40 #endif
41 
42 
43 //OOMPH-LIB headers
44 #include "../generic/Qelements.h"
45 #include "../generic/Telements.h"
46 
47 namespace oomph
48 {
49 
50 //======================================================================
51 ///A class for elements that allow the imposition of an applied traction
52 ///to the Navier--Stokes equations
53 ///The geometrical information can be read from the FaceGeometry<ELEMENT>
54 ///class and and thus, we can be generic enough without the need to have
55 ///a separate equations class
56 //======================================================================
57 template <class ELEMENT>
58 class NavierStokesTractionElement : public virtual FaceGeometry<ELEMENT>,
59  public virtual FaceElement
60 {
61 
62 private:
63 
64  ///Pointer to an imposed traction function
65  void (*Traction_fct_pt)(const double& time, const Vector<double> &x,
66  const Vector<double> &n, Vector<double> &result);
67 
68 protected:
69 
70 
71  /// \short The "global" intrinsic coordinate of the element when
72  /// viewed as part of a geometric object should be given by
73  /// the FaceElement representation, by default
74  double zeta_nodal(const unsigned &n, const unsigned &k,
75  const unsigned &i) const
76  {return FaceElement::zeta_nodal(n,k,i);}
77 
78  /// \short Access function that returns the local equation numbers
79  /// for velocity components.
80  /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
81  /// The default is to asssume that n is the local node number
82  /// and the i-th velocity component is the i-th unknown stored at the node.
83  virtual inline int u_local_eqn(const unsigned &n, const unsigned &i)
84  {return nodal_local_eqn(n,i);}
85 
86  ///\short Function to compute the shape and test functions and to return
87  ///the Jacobian of mapping
88  inline double shape_and_test_at_knot(const unsigned &ipt,
89  Shape &psi, Shape &test)
90  const
91  {
92  //Find number of nodes
93  unsigned n_node = nnode();
94  //Calculate the shape functions
95  shape_at_knot(ipt,psi);
96  //Set the test functions to be the same as the shape functions
97  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
98  //Return the value of the jacobian
99  return J_eulerian_at_knot(ipt);
100  }
101 
102 
103  ///Function to calculate the traction applied to the fluid
104  void get_traction(const double& time, const Vector<double> &x,
105  const Vector<double> &n, Vector<double> &result)
106  {
107  //If the function pointer is zero return zero
108  if(Traction_fct_pt == 0)
109  {
110  //Loop over dimensions and set body forces to zero
111  for(unsigned i=0;i<Dim;i++) {result[i] = 0.0;}
112  }
113  //Otherwise call the function
114  else
115  {
116  (*Traction_fct_pt)(time,x,n,result);
117  }
118  }
119 
120 
121  ///\short This function returns the residuals for the
122  /// traction function.
123  ///flag=1(or 0): do (or don't) compute the Jacobian as well.
125  Vector<double> &residuals,
126  DenseMatrix<double> &jacobian,
127  unsigned flag);
128 
129  ///The highest dimension of the problem
130  unsigned Dim;
131 
132 public:
133 
134  ///Constructor, which takes a "bulk" element and the value of the index
135  ///and its limit
137  const int &face_index,
138  const bool&
139  called_from_refineable_constructor=false) :
140  FaceGeometry<ELEMENT>(), FaceElement()
141  {
142 
143  //Attach the geometrical information to the element. N.B. This function
144  //also assigns nbulk_value from the required_nvalue of the bulk element
145  element_pt->build_face_element(face_index,this);
146 
147 #ifdef PARANOID
148  {
149  //Check that the element is not a refineable 3d element
150  if (!called_from_refineable_constructor)
151  {
152  //If it's three-d
153  if(element_pt->dim()==3)
154  {
155  //Is it refineable
156  RefineableElement* ref_el_pt=
157  dynamic_cast<RefineableElement*>(element_pt);
158  if(ref_el_pt!=0)
159  {
160  if (this->has_hanging_nodes())
161  {
162  throw OomphLibError(
163  "This flux element will not work correctly if nodes are hanging\n",
164  OOMPH_CURRENT_FUNCTION,
165  OOMPH_EXCEPTION_LOCATION);
166  }
167  }
168  }
169  }
170  }
171 #endif
172 
173  //Set the body force function pointer to zero
174  Traction_fct_pt = 0;
175 
176  //Set the dimension from the dimension of the first node
177  Dim = this->node_pt(0)->ndim();
178  }
179 
180  /// Destructor should not delete anything
182 
183  //Access function for the imposed traction pointer
184  void (*&traction_fct_pt())(const double& t, const Vector<double>& x,
185  const Vector<double> &n, Vector<double>& result)
186  {return Traction_fct_pt;}
187 
188  ///This function returns just the residuals
190  {
191  //Call the generic residuals function with flag set to 0
192  //using a dummy matrix argument
195  }
196 
197  ///This function returns the residuals and the jacobian
199  DenseMatrix<double> &jacobian)
200  {
201  //Call the generic routine with the flag set to 1
203  }
204 
205  ///Overload the output function
206  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
207 
208 ///Output function: x,y,[z],u,v,[w],p in tecplot format
209 void output(std::ostream &outfile, const unsigned &nplot)
210  {FiniteElement::output(outfile,nplot);}
211 
212 
213 };
214 
215 
216 
217 ///////////////////////////////////////////////////////////////////////
218 ///////////////////////////////////////////////////////////////////////
219 ///////////////////////////////////////////////////////////////////////
220 
221 
222 
223 //============================================================================
224 /// Function that returns the residuals for the imposed traction Navier_Stokes
225 /// equations
226 //============================================================================
227 template<class ELEMENT>
230  Vector<double> &residuals,
231  DenseMatrix<double> &jacobian,
232  unsigned flag)
233 {
234  //Find out how many nodes there are
235  unsigned n_node = nnode();
236 
237  // Get continuous time from timestepper of first node
238  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
239 
240  //Set up memory for the shape and test functions
241  Shape psif(n_node), testf(n_node);
242 
243  //Set the value of n_intpt
244  unsigned n_intpt = integral_pt()->nweight();
245 
246  //Integers to store local equation numbers
247  int local_eqn=0;
248 
249  //Loop over the integration points
250  for(unsigned ipt=0;ipt<n_intpt;ipt++)
251  {
252  //Get the integral weight
253  double w = integral_pt()->weight(ipt);
254 
255  //Find the shape and test functions and return the Jacobian
256  //of the mapping
257  double J = shape_and_test_at_knot(ipt,psif,testf);
258 
259  //Premultiply the weights and the Jacobian
260  double W = w*J;
261 
262  //Need to find position to feed into Traction function
263  Vector<double> interpolated_x(Dim);
264 
265  //Initialise to zero
266  for(unsigned i=0;i<Dim;i++) {interpolated_x[i] = 0.0;}
267 
268  //Calculate velocities and derivatives
269  for(unsigned l=0;l<n_node;l++)
270  {
271  //Loop over velocity components
272  for(unsigned i=0;i<Dim;i++) {interpolated_x[i] +=
273  nodal_position(l,i)*psif[l];}
274  }
275 
276  // Get the outer unit normal
277  Vector<double> interpolated_n(Dim);
278  outer_unit_normal(ipt,interpolated_n);
279 
280  //Get the user-defined traction terms
281  Vector<double> traction(Dim);
282  get_traction(time,interpolated_x,interpolated_n,traction);
283 
284  //Now add to the appropriate equations
285 
286  //Loop over the test functions
287  for(unsigned l=0;l<n_node;l++)
288  {
289  //Loop over the velocity components
290  for(unsigned i=0;i<Dim;i++)
291  {
292  local_eqn = u_local_eqn(l,i);
293  /*IF it's not a boundary condition*/
294  if(local_eqn >= 0)
295  {
296  //Add the user-defined traction terms
297  residuals[local_eqn] += traction[i]*testf[l]*W;
298 
299  //Assuming the the traction DOES NOT depend upon velocities
300  //or pressures, the jacobian is always zero, so no jacobian
301  //terms are required
302  }
303  } //End of loop over dimension
304  } //End of loop over shape functions
305 
306  }
307 
308 }
309 
310 
311 
312 /////////////////////////////////////////////////////////////////////////
313 /////////////////////////////////////////////////////////////////////////
314 /////////////////////////////////////////////////////////////////////////
315 
316 
317 
318 //======================================================================
319 /// A class for elements that allow the imposition of an applied traction
320 /// to the Navier--Stokes equations
321 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
322 /// class and and thus, we can be generic enough without the need to have
323 /// a separate equations class.
324 ///
325 /// THIS IS THE REFINEABLE VERSION.
326 //======================================================================
327 template <class ELEMENT>
329 public virtual NavierStokesTractionElement<ELEMENT>,
331 {
332  public:
333 
334  ///Constructor, which takes a "bulk" element and the face index
336  const int &face_index) :
337  // we're calling this from the constructor of the refineable version.
338  NavierStokesTractionElement<ELEMENT>(element_pt, face_index,true)
339  {}
340 
341  /// Destructor should not delete anything
343 
344 
345  /// \short Number of continuously interpolated values are the
346  /// same as those in the bulk element.
347  unsigned ncont_interpolated_values() const
348  {
349  return dynamic_cast<ELEMENT*>(this->bulk_element_pt())->
351  }
352 
353  ///This function returns just the residuals
355  {
356  //Call the generic residuals function using a dummy matrix argument
359  }
360 
361  ///\short This function returns the residuals and the Jacobian
363  DenseMatrix<double> &jacobian)
364  {
365  //Call the generic routine
367  jacobian,1);
368  }
369 
370 
371  protected:
372 
373  ///\short This function returns the residuals for the
374  /// traction function.
375  ///flag=1(or 0): do (or don't) compute the Jacobian as well.
377  Vector<double> &residuals,
378  DenseMatrix<double> &jacobian,
379  unsigned flag);
380 
381 };
382 
383 
384 ///////////////////////////////////////////////////////////////////////
385 ///////////////////////////////////////////////////////////////////////
386 ///////////////////////////////////////////////////////////////////////
387 
388 
389 
390 //============================================================================
391 /// Function that returns the residuals for the imposed traction Navier_Stokes
392 /// equations
393 //============================================================================
394 template<class ELEMENT>
397  Vector<double> &residuals,
398  DenseMatrix<double> &jacobian,
399  unsigned flag)
400 {
401 
402  // Get the indices at which the velocity components are stored
403  unsigned u_nodal_index[this->Dim];
404  for(unsigned i=0;i<this->Dim;i++)
405  {
406  u_nodal_index[i] = dynamic_cast<ELEMENT*>(
407  this->bulk_element_pt())->u_index_nst(i);
408  }
409 
410  //Find out how many nodes there are
411  unsigned n_node = nnode();
412 
413  // Get continuous time from timestepper of first node
414  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
415 
416  //Set up memory for the shape and test functions
417  Shape psif(n_node), testf(n_node);
418 
419  //Set the value of n_intpt
420  unsigned n_intpt = integral_pt()->nweight();
421 
422  //Integers to store local equation numbers
423  int local_eqn=0;
424 
425  //Loop over the integration points
426  for(unsigned ipt=0;ipt<n_intpt;ipt++)
427  {
428  //Get the integral weight
429  double w = integral_pt()->weight(ipt);
430 
431  //Find the shape and test functions and return the Jacobian
432  //of the mapping
433  double J = this->shape_and_test_at_knot(ipt,psif,testf);
434 
435  //Premultiply the weights and the Jacobian
436  double W = w*J;
437 
438  //Need to find position to feed into Traction function
439  Vector<double> interpolated_x(this->Dim);
440 
441  //Initialise to zero
442  for(unsigned i=0;i<this->Dim;i++) {interpolated_x[i] = 0.0;}
443 
444  //Calculate velocities and derivatives
445  for(unsigned l=0;l<n_node;l++)
446  {
447  //Loop over velocity components
448  for(unsigned i=0;i<this->Dim;i++)
449  {
450  interpolated_x[i] += nodal_position(l,i)*psif[l];
451  }
452  }
453 
454  // Get the outer unit normal
455  Vector<double> interpolated_n(this->Dim);
456  this->outer_unit_normal(ipt,interpolated_n);
457 
458  //Get the user-defined traction terms
459  Vector<double> traction(this->Dim);
460  this->get_traction(time,interpolated_x,interpolated_n,traction);
461 
462  //Now add to the appropriate equations
463 
464  //Number of master nodes and storage for the weight of the shape function
465  unsigned n_master=1; double hang_weight=1.0;
466 
467  //Pointer to hang info object
468  HangInfo* hang_info_pt=0;
469 
470  //Loop over the nodes for the test functions/equations
471  //----------------------------------------------------
472  for(unsigned l=0;l<n_node;l++)
473  {
474  //Local boolean to indicate whether the node is hanging
475  bool is_node_hanging = this->node_pt(l)->is_hanging();
476 
477  //If the node is hanging
478  if(is_node_hanging)
479  {
480  hang_info_pt = this->node_pt(l)->hanging_pt();
481 
482  //Read out number of master nodes from hanging data
483  n_master = hang_info_pt->nmaster();
484  }
485  //Otherwise the node is its own master
486  else
487  {
488  n_master = 1;
489  }
490 
491  //Loop over the master nodes
492  for(unsigned m=0;m<n_master;m++)
493  {
494  // Loop over velocity components for equations
495  for(unsigned i=0;i<this->Dim;i++)
496  {
497  //Get the equation number
498  //If the node is hanging
499  if(is_node_hanging)
500  {
501  //Get the equation number from the master node
502  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
503  u_nodal_index[i]);
504  //Get the hang weight from the master node
505  hang_weight = hang_info_pt->master_weight(m);
506  }
507  //If the node is not hanging
508  else
509  {
510  // Local equation number
511  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
512 
513  // Node contributes with full weight
514  hang_weight = 1.0;
515  }
516 
517  //If it's not a boundary condition...
518  if(local_eqn>= 0)
519  {
520 
521 /* //Loop over the test functions */
522 /* for(unsigned l=0;l<n_node;l++) */
523 /* { */
524 /* //Loop over the velocity components */
525 /* for(unsigned i=0;i<Dim;i++) */
526 /* { */
527 /* local_eqn = u_local_eqn(l,i); */
528 /* /\*IF it's not a boundary condition*\/ */
529 /* if(local_eqn >= 0) */
530 /* { */
531 
532 
533 
534  //Add the user-defined traction terms
535  residuals[local_eqn] += traction[i]*testf[l]*W*hang_weight;
536 
537  //Assuming the the traction DOES NOT depend upon velocities
538  //or pressures, the jacobian is always zero, so no jacobian
539  //terms are required
540  }
541  }
542  } //End of loop over dimension
543  } //End of loop over shape functions
544 
545  }
546 }
547 
548 
549 }
550 
551 #endif
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation numbers for velocity components. u_local_eqn(n,i) = local equation number or < 0 if pinned. The default is to asssume that n is the local node number and the i-th velocity component is the i-th unknown stored at the node.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping.
void fill_in_generic_residual_contribution_fluid_traction(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
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
void refineable_fill_in_generic_residual_contribution_fluid_traction(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
void(*&)(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &result) traction_fct_pt()
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
NavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index, const bool &called_from_refineable_constructor=false)
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
unsigned ncont_interpolated_values() const
Number of continuously interpolated values are the same as those in the bulk element.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
~RefineableNavierStokesTractionElement()
Destructor should not delete anything.
void output(std::ostream &outfile)
Overload the output function.
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
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
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
Class that contains data for hanging nodes.
Definition: nodes.h:684
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
unsigned Dim
The highest dimension of the problem.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the Jacobian.
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 ...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
~NavierStokesTractionElement()
Destructor should not delete anything.
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
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
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
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and 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
void get_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Function to calculate the traction applied to the fluid.
RefineableNavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the face index.