polar_stress_integral_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 the shear stress
31 //along either side wall
32 
33 #ifndef OOMPH_POLAR_STRESS_INTEGRAL_ELEMENTS_HEADER
34 #define OOMPH_POLAR_STRESS_INTEGRAL_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <config.h>
39 #endif
40 
41 
42 //OOMPH-LIB headers
43 #include "../generic/Qelements.h"
44 
45 namespace oomph
46 {
47 
48 //======================================================================
49 ///A class for elements that allow the imposition of an applied traction
50 ///to the Navier--Stokes equations
51 ///The geometrical information can be read from the FaceGeometery<ELEMENT>
52 ///class and and thus, we can be generic enough without the need to have
53 ///a separate equations class
54 //======================================================================
55 template <class ELEMENT>
56 class PolarStressIntegralElement : public virtual FaceGeometry<ELEMENT>,
57 public virtual FaceElement
58 {
59 
60 private:
61 
62  ///Pointer to an imposed traction function
63  void (*Traction_fct_pt)(const double& time, const Vector<double> &x,
64  Vector<double> &result);
65 
66  ///The highest dimension of the problem
67  unsigned Dim;
68 
69 protected:
70 
71 
72  /// \short Access function that returns the local equation numbers
73  /// for velocity components.
74  /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
75  /// The default is to asssume that n is the local node number
76  /// and the i-th velocity component is the i-th unknown stored at the node.
77  virtual inline int u_local_eqn(const unsigned &n, const unsigned &i)
78  {return nodal_local_eqn(n,i);}
79 
80  ///\short Function to compute the shape and test functions and to return
81  ///the Jacobian of mapping
82  inline double shape_and_test_at_knot(const unsigned &ipt,
83  Shape &psi, Shape &test)
84  const
85  {
86  //Find number of nodes
87  unsigned n_node = nnode();
88  //Calculate the shape functions
89  shape_at_knot(ipt,psi);
90  //Set the test functions to be the same as the shape functions
91  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
92  //Return the value of the jacobian
93  return J_eulerian_at_knot(ipt);
94  }
95 
96  /// Pointer to the angle alpha
97  double *Alpha_pt;
98 
99  //Traction elements need to know whether they're at the inlet or outlet
100  //as the unit outward normal has a differing sign dependent on
101  //the boundary
102  // phi=-1, phi=1
103  int Boundary;
104 
105 public:
106 
107  /// Alpha
108  const double &alpha() const {return *Alpha_pt;}
109 
110  /// Pointer to Alpha
111  double* &alpha_pt() {return Alpha_pt;}
112 
113  /// Boundary
114  const int boundary() const {return Boundary;}
115 
116  /// Function to set boundary
117  void set_boundary(int bound) {Boundary=bound;}
118 
119  ///Function to calculate the shear stress along boundary
120  double get_shear_stress();
121 
122  ///Constructor, which takes a "bulk" element and the value of the index
123  ///and its limit
125  const int &face_index) :
126  FaceGeometry<ELEMENT>(), FaceElement()
127  {
128  //Attach the geometrical information to the element. N.B. This function
129  //also assigns nbulk_value from the required_nvalue of the bulk element
130  element_pt->build_face_element(face_index,this);
131 
132 #ifdef PARANOID
133  {
134  //Check that the element is not a refineable 3d element
135  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
136 
137  //If it's three-d
138  if(elem_pt->dim()==3)
139  {
140  //Is it refineable
141  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
142  if(ref_el_pt!=0)
143  {
144  if (this->has_hanging_nodes())
145  {
146  throw OomphLibError(
147  "This flux element will not work correctly if nodes are hanging\n",
148  OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151  }
152  }
153  }
154 #endif
155 
156  //Set the dimension from the dimension of the first node
157  Dim = this->node_pt(0)->ndim();
158 
159  }
160 
161  /// Destructor should not delete anything
163 
164  ///This function returns just the residuals
166  {
167  //Do nothing
168  }
169 
170  ///This function returns the residuals and the jacobian
172  DenseMatrix<double> &jacobian)
173  {
174  //Do nothing
175  }
176 
177  ///\short Compute the element's residual Vector and the jacobian matrix
178  /// Plus the mass matrix especially for eigenvalue problems
180  Vector<double> &residuals,
181  DenseMatrix<double> &jacobian,DenseMatrix<double> &mass_matrix)
182  {
183  //Do nothing
184  }
185 
186  ///Overload the output function
187  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
188 
189 ///Output function: x,y,[z],u,v,[w],p in tecplot format
190 void output(std::ostream &outfile, const unsigned &nplot)
191  {FiniteElement::output(outfile,nplot);}
192 
193  /// local velocities
194  double u(const unsigned &l, const unsigned &i )
195  {return nodal_value(l,i);}
196 
197  /// local position
198  double x(const unsigned &l, const unsigned &i )
199  {return nodal_position(l,i);}
200 
201 };
202 
203 //============================================================================
204 /// Function that returns the shear stress
205 //============================================================================
206 template<class ELEMENT>
209 {
210  //Storage for shear stress
211  double dudphi,shear_contribution=0.0;
212 
213  //Set the value of n_intpt
214  unsigned n_intpt = integral_pt()->nweight();
215 
216  //Storage for local coordinate
217  Vector<double> s_local(1);
218  //Storage for local coordinate in bulk
219  Vector<double> s_bulk(2);
220 
221  //Find out how many nodes there are
222  unsigned n_node = nnode();
223 
224  //Set up memory for the shape and test functions
225  Shape psif(n_node), testf(n_node);
226 
227  //Loop over the integration points
228  for(unsigned ipt=0;ipt<n_intpt;ipt++)
229  {
230  //Get the integral weight
231  double w = integral_pt()->weight(ipt);
232 
233  //Find the shape and test functions and return the Jacobian
234  //of the mapping
235  double J = shape_and_test_at_knot(ipt,psif,testf);
236 
237  //Premultiply the weights and the Jacobian
238  double W = w*J;
239 
240  //Need to find position to feed into Traction function
241  Vector<double> interpolated_x(Dim);
242 
243  //Initialise to zero
244  for(unsigned i=0;i<Dim;i++)
245  {
246  interpolated_x[i] = 0.0;
247  }
248 
249  //Calculate velocities and derivatives
250  for(unsigned l=0;l<n_node;l++)
251  {
252  //Loop over velocity components
253  for(unsigned i=0;i<Dim;i++)
254  {
255  interpolated_x[i] += nodal_position(l,i)*psif[l];
256  }
257  }
258 
259  //Get the local coordinate
260  s_local[0] = integral_pt()->knot(ipt,0);
261 
262  //Get bulk coordinate
263  s_bulk = this->local_coordinate_in_bulk(s_local);
264 
265  // Upcast from GeneralisedElement to the present element
266  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->Bulk_element_pt);
267 
268  //Get du_dphi from bulk element
269  dudphi = el_pt->interpolated_dudx_pnst(s_bulk,0,1);
270 
271  //The contribution to the unweighted shear stress
272  shear_contribution += dudphi*W;
273  //dudphi*interpolated_x[0]*W;
274 
275  } //End of loop over integration points
276 
277  return shear_contribution;
278 
279 } //End of get_shear_stress
280 
281 } //End of namespace oomph
282 
283 #endif
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 fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
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 set_boundary(int bound)
Function to set boundary.
double x(const unsigned &l, const unsigned &i)
local position
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
~PolarStressIntegralElement()
Destructor should not delete anything.
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_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...
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 u(const unsigned &l, const unsigned &i)
local velocities
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.
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
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.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to an imposed traction function.
unsigned Dim
The highest dimension of the problem.
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 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
double get_shear_stress()
Function to calculate the shear stress along boundary.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void output(std::ostream &outfile)
Overload the output function.
double * Alpha_pt
Pointer to the angle alpha.
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
PolarStressIntegralElement(FiniteElement *const &element_pt, const int &face_index)