polar_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_POLAR_FLUID_TRACTION_ELEMENTS_HEADER
35 #define OOMPH_POLAR_FLUID_TRACTION_ELEMENTS_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41 
42 
43 //OOMPH-LIB headers
44 #include "../generic/Qelements.h"
45 
46 namespace oomph
47 {
48 
49 //======================================================================
50 ///A class for elements that allow the imposition of an applied traction
51 ///to the Navier--Stokes equations
52 ///The geometrical information can be read from the FaceGeometery<ELEMENT>
53 ///class and and thus, we can be generic enough without the need to have
54 ///a separate equations class
55 //======================================================================
56 template <class ELEMENT>
57 class PolarNavierStokesTractionElement : public virtual FaceGeometry<ELEMENT>,
58 public virtual FaceElement
59 {
60 
61 private:
62 
63  ///Pointer to an imposed traction function
64  void (*Traction_fct_pt)(const double& time, const Vector<double> &x,
65  Vector<double> &result);
66 
67  ///The highest dimension of the problem
68  unsigned Dim;
69 
70 protected:
71 
72 
73  /// \short Access function that returns the local equation numbers
74  /// for velocity components.
75  /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
76  /// The default is to asssume that n is the local node number
77  /// and the i-th velocity component is the i-th unknown stored at the node.
78  virtual inline int u_local_eqn(const unsigned &n, const unsigned &i)
79  {return nodal_local_eqn(n,i);}
80 
81  ///\short Function to compute the shape and test functions and to return
82  ///the Jacobian of mapping
83  inline double shape_and_test_at_knot(const unsigned &ipt,
84  Shape &psi, Shape &test)
85  const
86  {
87  //Find number of nodes
88  unsigned n_node = nnode();
89  //Calculate the shape functions
90  shape_at_knot(ipt,psi);
91  //Set the test functions to be the same as the shape functions
92  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
93  //Return the value of the jacobian
94  return J_eulerian_at_knot(ipt);
95  }
96 
97 
98  ///Function to calculate the traction applied to the fluid
99  void get_traction(double time, const Vector<double> &x,
100  Vector<double> &result)
101  {
102  //If the function pointer is zero return zero
103  if(Traction_fct_pt == 0)
104  {
105  //Loop over dimensions and set body forces to zero
106  for(unsigned i=0;i<Dim;i++) {result[i] = 0.0;}
107  }
108  //Otherwise call the function
109  else
110  {
111  (*Traction_fct_pt)(time,x,result);
112  }
113  }
114 
115  ///\short This function returns the residuals for the
116  /// traction function.
117  ///flag=1(or 0): do (or don't) compute the Jacobian as well.
119  DenseMatrix<double> &jacobian,
120  DenseMatrix<double> &mass_matrix,
121  unsigned flag);
122  /// Pointer to the angle alpha
123  double *Alpha_pt;
124 
125  /// \short Pointer to the Data item that stores the external pressure
127 
128  /// \short The Data that contains the traded pressure is stored
129  /// as external Data for the element. Which external Data item is it?
131 
132  //Traction elements need to know whether they're at the inlet or outlet
133  //as the unit outward normal has a differing sign dependent on
134  //the boundary
135  // -1=inlet, 1=outlet
136  int Boundary;
137 
138  //Pointer to homotopy parameter
139  double Eta;
140 
141 public:
142 
143  /// Alpha
144  const double &alpha() const {return *Alpha_pt;}
145 
146  /// Pointer to Alpha
147  double* &alpha_pt() {return Alpha_pt;}
148 
149  ///Function for setting up external pressure
150  void set_external_pressure_data(Data* pext_data_pt)
151  {
152  //Set external pressure pointer
153  Pext_data_pt=pext_data_pt;
154 
155  // Add to the element's external data so it gets included
156  // in the black-box local equation numbering scheme
158  this->add_external_data(pext_data_pt);
159  }
160 
161  /// Boundary
162  const int boundary() const {return Boundary;}
163 
164  /// Function to set boundary
165  void set_boundary(int bound) {Boundary=bound;}
166 
167  /// Eta
168  const double get_eta() const {return Eta;}
169 
170  /// Function to set Eta
171  void set_eta(double eta) {Eta=eta;}
172 
173  ///Constructor, which takes a "bulk" element and the value of the index
174  ///and its limit
176  const int &face_index) :
177  FaceGeometry<ELEMENT>(), FaceElement()
178  {
179 
180  //Attach the geometrical information to the element. N.B. This function
181  //also assigns nbulk_value from the required_nvalue of the bulk element
182  element_pt->build_face_element(face_index,this);
183 
184 #ifdef PARANOID
185  {
186  //Check that the element is not a refineable 3d element
187  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
188  //If it's three-d
189  if(elem_pt->dim()==3)
190  {
191  //Is it refineable
192  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
193  if(ref_el_pt!=0)
194  {
195  if (this->has_hanging_nodes())
196  {
197  throw OomphLibError(
198  "This flux element will not work correctly if nodes are hanging\n",
199  OOMPH_CURRENT_FUNCTION,
200  OOMPH_EXCEPTION_LOCATION);
201  }
202  }
203  }
204  }
205 #endif
206 
207  //Set the body force function pointer to zero
208  Traction_fct_pt = 0;
209 
210  //Set the external pressure pointer to be zero
211  Pext_data_pt = 0;
212 
213  //Set the dimension from the dimension of the first node
214  Dim = this->node_pt(0)->ndim();
215 
216  //Set Eta to one by default
217  Eta=1.0;
218 
219  }
220 
221  /// Destructor should not delete anything
223 
224  //Access function for the imposed traction pointer
225  void (*&traction_fct_pt())(const double& t, const Vector<double>& x,
226  Vector<double>& result)
227  {return Traction_fct_pt;}
228 
229  ///This function returns just the residuals
231  {
232  //Call the generic residuals function with flag set to 0
233  //using a dummy matrix argument
236  }
237 
238  ///This function returns the residuals and the jacobian
240  DenseMatrix<double> &jacobian)
241  {
242  //Call the generic routine with the flag set to 1
244  }
245 
246  ///\short Compute the element's residual Vector and the jacobian matrix
247  /// Plus the mass matrix especially for eigenvalue problems
249  Vector<double> &residuals,
250  DenseMatrix<double> &jacobian,DenseMatrix<double> &mass_matrix)
251  {
252  //Call the generic routine with the flag set to 2
254  }
255 
256  ///Overload the output function
257  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
258 
259 ///Output function: x,y,[z],u,v,[w],p in tecplot format
260 void output(std::ostream &outfile, const unsigned &nplot)
261  {FiniteElement::output(outfile,nplot);}
262 
263  /// local velocities
264  double u(const unsigned &l, const unsigned &i )
265  {return nodal_value(l,i);}
266 
267  /// local position
268  double x(const unsigned &l, const unsigned &i )
269  {return nodal_position(l,i);}
270 
271 };
272 
273 
274 
275 ///////////////////////////////////////////////////////////////////////
276 ///////////////////////////////////////////////////////////////////////
277 ///////////////////////////////////////////////////////////////////////
278 
279 
280 
281 //============================================================================
282 /// Function that returns the residuals for the imposed traction Navier_Stokes
283 /// equations
284 //============================================================================
285 template<class ELEMENT>
288  DenseMatrix<double> &jacobian,
289  DenseMatrix<double> &mass_matrix,
290  unsigned flag)
291 {
292  //Find out how many nodes there are
293  unsigned n_node = nnode();
294 
295  // Get continuous time from timestepper of first node
296  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
297 
298  //Set up memory for the shape and test functions
299  Shape psif(n_node), testf(n_node);
300 
301  //Set the value of n_intpt
302  unsigned n_intpt = integral_pt()->nweight();
303 
304  //Get Alpha
305  const double Alpha = alpha();
306 
307  //Storage for external pressure
308  double pext=0.0;
309 
310  //Get boundary multiplier
311  //This is necessary because the sign of the traction is
312  //dependent on the boundary
313  const int multiplier = boundary();
314 
315  //Get the homotopy parameter (if necessary)
316  const double eta=get_eta();
317 
318  //Integers to store local equation numbers
319  int local_eqn=0,local_unknown=0,pext_local_eqn=0,pext_local_unknown=0;
320 
321  ////////////////////////////////////////NEW//////////////////////////////////////////
322 
323  //Get local equation number of external pressure
324  //Note that if we have not passed an external pressure pointer to this element
325  //(and at the same time added it to the element's external data)
326  //than this will be -1 to indicate that it is not a degree of freedom here
327  if(Pext_data_pt==0)
328  {
329  pext_local_eqn=-1;
330  }
331  else
332  {
333  //If at a non-zero degree of freedom add in the entry
334  pext_local_eqn = external_local_eqn(External_data_number_of_Pext,0);
335 
336  //Get external pressure
337  pext=Pext_data_pt->value(0);
338  }
339 
340  //The local unkown number of pext will be the same
341  pext_local_unknown=pext_local_eqn;
342 
343  ////////////////////////////////////////NEW//////////////////////////////////////////
344 
345  //Loop over the integration points
346  for(unsigned ipt=0;ipt<n_intpt;ipt++)
347  {
348  //Get the integral weight
349  double w = integral_pt()->weight(ipt);
350 
351  //Find the shape and test functions and return the Jacobian
352  //of the mapping
353  double J = shape_and_test_at_knot(ipt,psif,testf);
354 
355  //Premultiply the weights and the Jacobian
356  double W = w*J;
357 
358  //Need to find position to feed into Traction function
359  Vector<double> interpolated_x(Dim);
360  Vector<double> interpolated_u(Dim);
361 
362  //Initialise to zero
363  for(unsigned i=0;i<Dim;i++)
364  {
365  interpolated_x[i] = 0.0;
366  interpolated_u[i] = 0.0;
367  }
368 
369  //Calculate velocities and derivatives:
370  // Loop over nodes
371  for(unsigned l=0;l<n_node;l++)
372  {
373  //Loop over directions
374  for(unsigned i=0;i<Dim;i++)
375  {
376  //Get the nodal value
377  interpolated_u[i] += this->nodal_value(l,i)*psif[l];
378  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
379  }
380  }
381 
382  //Get the user-defined traction terms
383  Vector<double> traction(Dim);
384  get_traction(time,interpolated_x,traction);
385 
386  //Now add to the appropriate equations
387 
388  //Loop over the test functions
389  for(unsigned l=0;l<n_node;l++)
390  {
391  //Only alter u velocity component
392  {
393  unsigned i=0;
394  local_eqn = u_local_eqn(l,i);
395  /*IF it's not a boundary condition*/
396  if(local_eqn >= 0)
397  {
398  //Add the user-defined traction terms
399  residuals[local_eqn] -= multiplier*eta*3.0*(interpolated_u[i]/interpolated_x[0])
400  *testf[l]*interpolated_x[0]*Alpha*W;
401 
402  ////////////////////////////////////////NEW//////////////////////////////////////////
403 
404  //Plus additional external pressure contribution at inlet
405  //This is zero if we haven't passed a Pext_data_pt to the element
406  residuals[local_eqn] += pext*testf[l]*interpolated_x[0]*Alpha*W;
407 
408  ////////////////////////////////////////NEW//////////////////////////////////////////
409 
410  //CALCULATE THE JACOBIAN
411  if(flag)
412  {
413  //Loop over the velocity shape functions again
414  for(unsigned l2=0;l2<n_node;l2++)
415  {
416  //We only have an i2=0 contribution
417  unsigned i2=0;
418  {
419  //If at a non-zero degree of freedom add in the entry
420  local_unknown = u_local_eqn(l2,i2);
421  if(local_unknown >= 0)
422  {
423  //Add contribution to Elemental Matrix
424  jacobian(local_eqn,local_unknown)
425  -= multiplier*eta*3.0*(psif[l2]/interpolated_x[0])
426  *testf[l]*interpolated_x[0]*Alpha*W;
427 
428  } //End of (Jacobian's) if not boundary condition statement
429  } //End of i2=0 section
430  } //End of l2 loop
431 
432  ////////////////////////////////////////NEW//////////////////////////////////////////
433  //Add pext's contribution to these residuals
434  //This only needs to be done once hence why it is outside the l2 loop
435  if(pext_local_unknown >= 0)
436  {
437  //Add contribution to Elemental Matrix
438  jacobian(local_eqn,pext_local_unknown) += testf[l]*interpolated_x[0]*Alpha*W;
439  }
440  ////////////////////////////////////////NEW//////////////////////////////////////////
441 
442  } /*End of Jacobian calculation*/
443 
444  } //end of if not boundary condition statement
445  } //End of i=0 section
446 
447  } //End of loop over shape functions
448 
449  ////////////////////////////////////////NEW//////////////////////////////////////////
450 
451  /// The additional residual for the mass flux
452  /// (ie. the extra equation for pext)
453  /// This is an integral equation along the whole boundary
454  /// It lies outside the loop over shape functions above
455  {
456  /*IF it's not a boundary condition*/
457  if(pext_local_eqn >= 0)
458  {
459  //Add the user-defined traction terms
460  residuals[pext_local_eqn] += interpolated_u[0]*interpolated_x[0]*Alpha*W;
461 
462  //No longer necessary due to my FluxCosntraint element
463  //Now take off a fraction of the desired mass flux
464  //Divided by number of elements and number of int points in each
465  //HACK
466  //residuals[pext_local_eqn] -= (1.0/(30.*3.));
467  //HACK
468 
469  //CALCULATE THE JACOBIAN
470  if(flag)
471  {
472  //Loop over the velocity shape functions again
473  for(unsigned l2=0;l2<n_node;l2++)
474  {
475  //We only have an i2=0 contribution
476  unsigned i2=0;
477  {
478  //If at a non-zero degree of freedom add in the entry
479  local_unknown = u_local_eqn(l2,i2);
480  if(local_unknown >= 0)
481  {
482  //Add contribution to Elemental Matrix
483  jacobian(pext_local_eqn,local_unknown)
484  += psif[l2]*interpolated_x[0]*Alpha*W;
485 
486  } //End of (Jacobian's) if not boundary condition statement
487  } //End of i2=0 section
488 
489  } //End of l2 loop
490  } /*End of Jacobian calculation*/
491 
492  } //end of if not boundary condition statement
493  } //End of additional residual for the mass flux
494 
495  ////////////////////////////////////////NEW//////////////////////////////////////////
496 
497  } //End of loop over integration points
498 
499 } //End of fill_in_generic_residual_contribution
500 
501 } //End of namespace oomph
502 
503 #endif
void set_external_pressure_data(Data *pext_data_pt)
Function for setting up external pressure.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to an imposed traction function.
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
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eig...
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
Data * Pext_data_pt
Pointer to the Data item that stores the external pressure.
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
void(*&)(const double &t, const Vector< double > &x, Vector< double > &result) traction_fct_pt()
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
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
char t
Definition: cfortran.h:572
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
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.
unsigned External_data_number_of_Pext
The Data that contains the traded pressure is stored as external Data for the element. Which external Data item is it?
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
PolarNavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index)
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
double x(const unsigned &l, const unsigned &i)
local position
unsigned Dim
The highest dimension of the problem.
void get_traction(double time, const Vector< double > &x, Vector< double > &result)
Function to calculate the traction applied to the fluid.
double u(const unsigned &l, const unsigned &i)
local velocities
~PolarNavierStokesTractionElement()
Destructor should not delete anything.
double * Alpha_pt
Pointer to the angle alpha.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
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
void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
void set_boundary(int bound)
Function to set boundary.
void output(std::ostream &outfile)
Overload the output function.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void set_eta(double eta)
Function to set Eta.
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_residuals(Vector< double > &residuals)
This function returns just the residuals.
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
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.