navier_stokes_surface_drag_torque_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 specific surface elements
31 
32 #ifndef OOMPH_NAVIER_STOKES_SURFACE_DRAG_TORQUE_ELEMENTS_HEADER
33 #define OOMPH_NAVIER_STOKES_SURFACE_DRAG_TORQUE_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 namespace oomph
41 {
42 
43 //======================================================================
44 /// A class of elements that allow the determination of the
45 /// drag and toque, relative to a given centre of rotation, along
46 /// a domain boundary.
47 /// The element operates as a FaceElement and attaches itself
48 /// to a bulk element of the type specified by the template
49 /// argument.
50 //======================================================================
51 template <class ELEMENT>
53  public virtual FaceGeometry<ELEMENT>,
54  public virtual FaceElement,
55  public virtual ElementWithDragFunction
56  {
57 
58  public:
59 
60  ///\short Constructor, which takes a "bulk" element and the value of an index
61  ///describing to which face the element should be attached.
63  const int &face_index) :
64  FaceGeometry<ELEMENT>(), FaceElement()
65  {
66  //Attach the geometrical information to the element. N.B. This function
67  //also assigns nbulk_value from the required_nvalue of the bulk element
68  element_pt->build_face_element(face_index,this);
69 
70  //Set the dimension from the dimension of the first node
71  this->Dim = node_pt(0)->ndim();
72 
73  //Default centre of rotation is the origin
74  this->Centre_of_rotation.resize(this->Dim,0.0);
75  }
76 
77  /// \short Set the translation and rotation of the rigid object
78  /// as external data
79  void set_translation_and_rotation(Data* const &object_data_pt)
80  {
81  this->Translation_index = this->add_external_data(object_data_pt);
82  }
83 
84 
85  /// \short Access function for the centre of rotation
86  double &centre_of_rotation(const unsigned &i)
87  {return this->Centre_of_rotation[i];}
88 
89 
90  /// \short Function that specifies the drag force and the torque about
91  /// the origin
92  virtual void get_drag_and_torque(Vector<double>& drag_force,
93  Vector<double>& drag_torque)
94  {
95  // Spatial dimension of element
96  unsigned ndim=dim();
97 
98  // Initialise force
99  for(unsigned i=0;i<ndim+1;i++) {drag_force[i] = 0.0;}
100 
101  //Now there will be one torque component in 2D (1D surface element)
102  //and three in 3D (2D surface element)
103  for(unsigned i=0;i<2*ndim-1;i++) {drag_torque[i] = 0.0;}
104 
105  //Vector of local coordinates in face element
106  Vector<double> s(ndim);
107 
108  // Vector for global Eulerian coordinates
109  Vector<double> x(ndim+1);
110 
111  // Vector for local coordinates in bulk element
112  Vector<double> s_bulk(ndim+1);
113 
114  //Set the value of n_intpt
115  unsigned n_intpt = integral_pt()->nweight();
116 
117  // Get pointer to assocated bulk element
118  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
119 
120  // Loop over the integration points
121  for (unsigned ipt=0;ipt<n_intpt;ipt++)
122  {
123 
124  //Assign values of s in FaceElement and local coordinates in bulk element
125  for(unsigned i=0;i<ndim;i++)
126  {
127  s[i] = integral_pt()->knot(ipt,i);
128  }
129 
130  //Get the bulk coordinates
131  this->get_local_coordinate_in_bulk(s,s_bulk);
132 
133  //Get the integral weight
134  double w = integral_pt()->weight(ipt);
135 
136  // Get jacobian of mapping
137  double J=J_eulerian(s);
138 
139  //Premultiply the weights and the Jacobian
140  double W = w*J;
141 
142  // Get x position as Vector
143  interpolated_x(s,x);
144 
145 #ifdef PARANOID
146 
147  // Get x position as Vector from bulk element
148  Vector<double> x_bulk(ndim+1);
149  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
150 
151  double max_legal_error=1.0e-5;
152  double error=0.0;
153  for(unsigned i=0;i<ndim+1;i++)
154  {
155  error+=fabs(x[i]- x_bulk[i]);
156  }
157  if (error>max_legal_error)
158  {
159  std::ostringstream error_stream;
160  error_stream << "difference in Eulerian posn from bulk and face: "
161  << error << " exceeds threshold " << max_legal_error
162  << std::endl;
163  throw OomphLibError(
164  error_stream.str(),
165  OOMPH_CURRENT_FUNCTION,
166  OOMPH_EXCEPTION_LOCATION);
167  }
168 #endif
169 
170  // Outer unit normal of the fluid
171  Vector<double> normal(ndim+1);
172  outer_unit_normal(s,normal);
173 
174  // Get velocity from bulk element
175  Vector<double> veloc(ndim+1);
176  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
177 
178  // Get traction from bulk element this is directed into the fluid
179  //so we need to reverse the sign
180  Vector<double> traction(ndim+1);
181  bulk_el_pt->get_traction(s_bulk,normal,traction);
182 
183  // Integrate (note the minus sign)
184  for (unsigned i=0;i<ndim+1;i++)
185  {
186  drag_force[i]-=traction[i]*W;
187  }
188 
189  //Now Calculate the torque which is r x F
190  //Scale X relative to the centre of rotation
191  Vector<double> X(ndim+1);
192  for(unsigned i=0;i<ndim+1;i++)
193  {X[i] = x[i] - Centre_of_rotation[i]
195 
196  //In 2D it's just a single scalar
197  //again note the minus sign
198  if(ndim==1)
199  {
200  drag_torque[0] -= (X[0]*traction[1] - X[1]*traction[0])*W;
201  }
202  else if(ndim==2)
203  {
204  drag_torque[0] -= (X[1]*traction[2] - X[2]*traction[1])*W;
205  drag_torque[1] -= (X[2]*traction[0] - X[0]*traction[2])*W;
206  drag_torque[2] -= (X[0]*traction[1] - X[1]*traction[0])*W;
207  }
208  else
209  {
210  throw OomphLibError(
211  "Dimension of a surface element must be 1 or 2\n",
212  OOMPH_CURRENT_FUNCTION,
213  OOMPH_EXCEPTION_LOCATION);
214  }
215  }
216 
217  //std::cout << "Forces " << drag_force[0] << " " << drag_force[1] << " "
218  // << drag_torque[0] << "\n";
219  }
220 
221 
222  /// \short Specify the value of nodal zeta from the face geometry
223  /// The "global" intrinsic coordinate of the element when
224  /// viewed as part of a geometric object should be given by
225  /// the FaceElement representation, by default (needed to break
226  /// indeterminacy if bulk element is SolidElement)
227  double zeta_nodal(const unsigned &n, const unsigned &k,
228  const unsigned &i) const
229  {return FaceElement::zeta_nodal(n,k,i);}
230 
231 
232 
233  /// \short Output function
234  void output(std::ostream &outfile, const unsigned &n_plot)
235  {
236  // Get pointer to assocated bulk element
237  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
238 
239  // Elemental dimension
240  unsigned dim_el=dim();
241 
242  //Local coordinates
243  Vector<double> s(dim_el);
244 
245  Vector<double> s_bulk(dim_el+1);
246 
247  //Calculate the Eulerian coordinates and Lagrange multiplier
248  Vector<double> x(dim_el+1,0.0);
249 
250  // Outer unit normal of the fluid
251  Vector<double> normal(dim_el+1);
252 
253  // Velocity from bulk element
254  Vector<double> veloc(dim_el+1);
255 
256  // Tractions (from bulk element)
257  Vector<double> traction(dim_el+1);
258 
259  // Tecplot header info
260  outfile << this->tecplot_zone_string(n_plot);
261 
262  // Loop over plot points
263  unsigned num_plot_points=this->nplot_points(n_plot);
264  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
265  {
266  // Get local coordinates of plot point
267  this->get_s_plot(iplot,n_plot,s);
268 
269  this->get_local_coordinate_in_bulk(s,s_bulk);
270 
271  //Get x position from bulk
272  bulk_el_pt->interpolated_x(s_bulk,x);
273  //Get outer unit normal
274  outer_unit_normal(s,normal);
275  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
276  // Get traction from bulk element this is directed into the fluid
277  //so we need to reverse the sign
278  bulk_el_pt->get_traction(s_bulk,normal,traction);
279 
280  //Now Calculate the torque which is r x F
281  //Scale X relative to the centre of rotation
282  Vector<double> X(dim_el+1);
283  for(unsigned i=0;i<dim_el+1;i++)
284  {X[i] = x[i] - Centre_of_rotation[i]
286 
287  for(unsigned i=0;i<dim_el+1;i++)
288  {
289  outfile << x[i] << " ";
290  }
291 
292  for(unsigned i=0;i<dim_el+1;i++)
293  {
294  outfile << normal[i] << " ";
295  }
296 
297  for(unsigned i=0;i<dim_el+1;i++)
298  {
299  outfile << veloc[i] << " ";
300  }
301 
302  for(unsigned i=0;i<dim_el+1;i++)
303  {
304  outfile << -1.0*traction[i] << " ";
305  }
306  outfile << -(X[0]*traction[1] - X[1]*traction[0]);
307 
308  outfile << std::endl;
309 
310  }
311 
312  this->write_tecplot_zone_footer(outfile,n_plot);
313  }
314 
315 
316 private:
317 
318  ///The highest dimension of the problem
319  unsigned Dim;
320 
321  ///The centre of rotation for the torque calculation
323 
324  /// The index of where the translation and rotation data is stored
326 
327 
328 };
329 
330 
331 }
332 
333 #endif
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
virtual void get_drag_and_torque(Vector< double > &drag_force, Vector< double > &drag_torque)
Function that specifies the drag force and the torque about the origin.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:185
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
double & centre_of_rotation(const unsigned &i)
Access function for the centre of rotation.
A general Finite Element class.
Definition: elements.h:1271
unsigned Translation_index
The index of where the translation and rotation data is stored.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
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.
Vector< double > Centre_of_rotation
The centre of rotation for the torque calculation.
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6064
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
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
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:662
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
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
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
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 set_translation_and_rotation(Data *const &object_data_pt)
Set the translation and rotation of the rigid object as external data.
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
NavierStokesSurfaceDragTorqueElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of an index describing to which face the elem...