immersed_rigid_body_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 Rigid Body Elements that are immersed in an external fluid
31 #ifndef OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
32 #define OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 #include "../generic/elements.h"
41 #include "../generic/triangle_mesh.h"
42 #include "../generic/fsi.h"
43 
44 namespace oomph
45 {
46 
47 //=====================================================================
48 // DOXYERROR: I removed a \f from the formula \f$\mbox{\boldmath$F$} = \mbox{\boldmath$F$}^{*}/\mu U\f$ below because it broke doxygen. Someone who knows the maths should check that this still makes sense.
49 /// Class that solves the equations of motion for a general
50 /// two-dimensional rigid body subject to a particular imposed force and torque
51 /// distribution and immersed within an external fluid.
52 /// The body's position is entirely specified by the location of its
53 /// centre of mass, \f$\mbox{\boldmath$X$}\f$,
54 /// and a single angle, \f$\phi\f$, that represents a possible rotation.
55 /// The equations of motion are then simply Newton's second law for the
56 /// conservation of linear momentum in two directions and angular
57 /// momentum about the single possible axis of rotation.
58 ///
59 /// The non-dimensionalisation is based on the viscous scales of
60 /// the surrounding fluid, in which case, the governing equations are
61 /// \f[ Re St^{2} M \frac{\mbox{d}^{2} \mbox{\boldmath$X$}}{\mbox{d} t^{2}}
62 /// = \mbox{\boldmath$F$} + \oint \mbox{\boldmath$t$} \mbox{d} s
63 /// \quad\mbox{and}\quad Re St^{2} I \frac{\mbox{d}^{2}\Phi}{\mbox{d} t^{2}}
64 /// = \mbox{\boldmath$T$} + \oint \mbox{\boldmath$q$} \mbox{d} s,
65 /// \f]
66 /// where \f$\mbox{\boldmath$F$} = \mbox{\boldmath$F$}^{*}/\mu U\f$ is the
67 /// external force per unit length;
68 /// \f$\mbox{\boldmath$t$}\f$ is the net force per unit length
69 /// applied on the rigid body
70 /// by the surrounding fluid;
71 /// \f$\mbox{\boldmath$T$} = \mbox{\boldmath$T$}^{*}/(\mu UL)\f$
72 /// is the external torque per unit length; and \f$\mbox{\boldmath$q$}\f$ is
73 /// the net torque per unit length applied on the body by the fluid.
74 /// \f$M\f$ is a scaled mass the density ratio multiplied by a shape parameter
75 /// and \f$I\f$ is a scaled moment of inertia. Finally, \f$Re\f$ and \f$St\f$
76 /// are the Reynolds and Strouhal numbers of the surrounding fluid.
77 /// Note that these equations may be used without the external fluid,
78 /// in which case the non-dimensionalisation doesn't make a lot of sense,
79 /// but may be re-interpreted accordingly.
80 ///
81 /// A Data object whose three values
82 /// represent the x and y displacements of the body's centre of mass and
83 /// its rotation about the centre of mass may be passed in as external
84 /// data or, if not, it will be constructed internally.
85 ///
86 /// For general usage, an underlying geometric object must passed to the
87 /// constructor and the position will be determined by applying rigid body
88 /// motions to the underlying object based on the initial_centre_of_mass
89 /// and initial_phi (angle), which default to zero. If these defaults are
90 /// not suitable, the values must be set externally. In addition a mass
91 /// and moment of inertia should also be set externally.
92 ///
93 /// If added to a mesh in the Problem (in its incarnation as a
94 /// GeneralisedElement) the displacement/rotation of the body
95 /// is computed in response to (i) user-specifiable applied forces
96 /// and a torque and (ii) the net drag (and associated torque) from
97 /// a mesh of elements that can exert a drag onto the body (typically
98 /// Navier-Stokes FaceElements that apply a viscous drag to an
99 /// immersed body, represented by the body.)
100 //=====================================================================
102  {
103 
104  public:
105 
106  /// \short Function pointer to function that specifies
107  /// external force
108  typedef void (*ExternalForceFctPt)(const double& time,
109  Vector<double>& external_force);
110 
111  /// \short Function pointer to function that specifies
112  /// external torque
113  typedef void (*ExternalTorqueFctPt)(const double& time,
114  double& external_torque);
115 
116  protected:
117 
118  /// \short Default constructor that intialises everything to
119  /// zero. This is expected to be called only from derived clases
120  /// such as the ImmersedRigidBodyTriangleMeshPolygon that can provided
121  /// their own position() functions.
123  Data* const &centre_displacement_data_pt=0) :
124  Initial_Phi(0.0), Mass(0.0), Moment_of_inertia(0.0),
126  Geom_object_pt(0),
134  {
135  this->initialise(time_stepper_pt);
136  }
137 
138  public:
139 
140  /// \short Constructor that takes an underlying geometric object:
141  /// and timestepper.
142  ImmersedRigidBodyElement(GeomObject* const &geom_object_pt,
144  Data* const &centre_displacement_data_pt=0) :
145  Initial_Phi(0.0), Mass(0.0), Moment_of_inertia(0.0),
147  Geom_object_pt(geom_object_pt),
155  {
156  this->initialise(time_stepper_pt);
157  }
158 
159  /// Set the rotation of the object to be included
161 
162  /// \short Set the rotation of the object to be ignored (only really
163  /// useful if you have a circular shape)
165 
166  ///Access function for the initial angle
167  double &initial_phi() {return Initial_Phi;}
168 
169  ///Access function for the initial centre of mass
170  double &initial_centre_of_mass(const unsigned &i)
171  {return Initial_centre_of_mass[i];}
172 
173  ///Access function for the initial centre of mass (const version)
174  const double &initial_centre_of_mass(const unsigned &i) const
175  {return Initial_centre_of_mass[i];}
176 
177  ///Overload the position to apply the rotation and translation
178  void position(const Vector<double> &xi, Vector<double> &r) const
179  {
180  Vector<double> initial_x(2);
181  Geom_object_pt->position(xi,initial_x);
182  this->apply_rigid_body_motion(0,initial_x,r);
183  }
184 
185  ///Overload to include the time history of the motion of the object
186  void position(const unsigned& t, const Vector<double>& xi,
187  Vector<double>& r) const
188  {
189  Vector<double> initial_x(2);
190  Geom_object_pt->position(xi,initial_x);
191  this->apply_rigid_body_motion(t,initial_x,r);
192  }
193 
194 
195  ///Work out the position derivative, including rigid body motion
196  void dposition_dt(const Vector<double> &zeta, const unsigned &j,
197  Vector<double> &drdt);
198 
199  /// \short Destuctor: Cleanup if required
201  {
203  {
205  //Null out the pointer stored as internal data
206  this->internal_data_pt(0) = 0;
207  }
208  }
209 
210  /// Access to dimensionless "mass" shape parameter that must be set by hand
211  /// for non polygonal shapes
212  double &mass_shape() {return Mass;}
213 
214  /// Access to dimensionless polar "moment of inertia" shape parameter
216 
217  /// \short Pointer to Data for centre of gravity displacement.
218  /// Values: 0: x-displ; 1: y-displ; 2: rotation angle.
221 
222  /// x-displacement of centre of mass
224  {return *(Centre_displacement_data_pt->value_pt(0));}
225 
226  /// y-displacement of centre of mass
228  {return *(Centre_displacement_data_pt->value_pt(1));}
229 
230  /// rotation of centre of mass
232  {return *(Centre_displacement_data_pt->value_pt(2));}
233 
234  /// Get current centre of gravity
236  {
237  Vector<double> cog(2);
238  for(unsigned i=0;i<2;i++)
239  {
240  cog[i] = Initial_centre_of_mass[i] +
242  }
243  return cog;
244  }
245 
246  /// Pin the i-th coordinate of the centre of mass
247  void pin_centre_of_mass_coordinate(const unsigned& i)
248  {
250  }
251 
252  /// Unpin the i-th coordinate of the centre of mass
253  void unpin_centre_of_mass_coordinate(const unsigned& i)
254  {
256  }
257 
258  /// Pin the rotation angle
260  {
262  }
263 
264  /// Unpin the rotation angle
266  {
268  }
269 
270  /// Output position velocity and acceleration of centre of gravity
271  void output_centre_of_gravity(std::ostream& outfile);
272 
273  /// Get the contribution to the residuals
275  {
276  // Get generic function without jacobian terms
279  }
280 
281 
282  /// Get residuals including contribution to jacobian
284  DenseMatrix<double>& jacobian)
285  {
286  // Get generic function, but don't bother to get jacobian terms
287  bool flag = false;
288  get_residuals_rigid_body_generic(residuals,jacobian,flag);
289  //Get the internal effects by FD, because a change in internal
290  //data can change the fluid loading, so this is one of those
291  //subtle interactions
292  this->fill_in_jacobian_from_internal_by_fd(residuals,jacobian);
293  //Get the effect of the fluid loading on the rigid body
294  this->fill_in_jacobian_from_external_by_fd(residuals,jacobian);
295  }
296 
297  /// \short Update the positions of the nodes in fluid elements
298  /// adjacent to the rigid body, defined as being elements in the
299  /// drag mesh.
301  {
302  //If there is no drag mesh then do nothing
303  if (Drag_mesh_pt==0) {return;}
304  //Otherwise update the elements adjacent to the rigid body
305  //(all the elements adjacent to those in the drag mesh)
306  else
307  {
308  unsigned nel=Drag_mesh_pt->nelement();
309  for (unsigned e=0;e<nel;e++)
310  {
311  dynamic_cast<FaceElement*>(Drag_mesh_pt->element_pt(e))->
312  bulk_element_pt()->node_update();
313  }
314  }
315  }
316 
317 
318  /// \short After an external data change, update the nodal positions
319  inline void update_in_external_fd(const unsigned &i)
321 
322  ///\short Do nothing to reset within finite-differencing of external data
323  inline void reset_in_external_fd(const unsigned &i) { }
324 
325  ///\short After all external data finite-differencing, update nodal positions
326  inline void reset_after_external_fd()
328 
329  ///\short After an internal data change, update the nodal positions
330  inline void update_in_internal_fd(const unsigned &i)
332 
333  ///\short Do nothing to reset within finite-differencing of internal data
334  inline void reset_in_internal_fd(const unsigned &i) { }
335 
336  ///\short After all internal data finite-differencing, update nodal positions
337  inline void reset_after_internal_fd()
339 
340  /// \short Get force and torque from specified fct pointers and
341  /// drag mesh
342  void get_force_and_torque(const double& time,
343  Vector<double>& force,
344  double& torque);
345 
346  /// \short Access to function pointer to function that specifies
347  /// external force
349 
350  /// \short Access to function pointer to function that specifies
351  /// external torque
353 
354  /// \short Access fct to mesh containing face elements that allow
355  /// the computation of the drag on the body
356  Mesh* const &drag_mesh_pt() {return Drag_mesh_pt;}
357 
358  /// \short Function to set the drag mesh and add the appropriate load
359  /// and geometric data as external data to the Rigid Body
360  void set_drag_mesh(Mesh* const &drag_mesh_pt);
361 
362  /// \short Function to clear the drag mesh and all associated external data
364  {
365  //Delete the hijacked data created as the load
367  //Flush the external data
368  this->flush_external_data();
369  //Set the Drag_mesh pointer to null
370  Drag_mesh_pt = 0;
371  }
372 
373  /// The position of the object depends on one data item
374  unsigned ngeom_data() const {return 1;}
375 
376  /// \short Return pointer to the j-th (only) Data item that the object's
377  /// shape depends on.
378  Data* geom_data_pt(const unsigned& j)
379  {return this->Centre_displacement_data_pt;}
380 
381  /// \short Access function to the direction of gravity
382  Vector<double>* &g_pt() {return G_pt;}
383 
384  /// \short Access function for gravity
385  const Vector<double> &g() const {return *G_pt;}
386 
387  /// \short Access function for the pointer to the fluid Reynolds number
388  double* &re_pt() {return Re_pt;}
389 
390  ///Access function for the fluid Reynolds number
391  const double &re() const{return *Re_pt;}
392 
393  /// \short Access function for the pointer to the fluid Strouhal number
394  double* &st_pt() {return St_pt;}
395 
396  ///Access function for the fluid Strouhal number
397  const double &st() const {return *St_pt;}
398 
399  /// \short Access function for pointer to the fluid inverse Froude number
400  /// (dimensionless gravitational loading)
401  double* &re_invfr_pt() {return ReInvFr_pt;}
402 
403  /// Access to the fluid inverse Froude number
404  const double &re_invfr() {return *ReInvFr_pt;}
405 
406  /// \short Access function for the pointer to the density ratio
407  double* &density_ratio_pt() {return Density_ratio_pt;}
408 
409  ///Access function for the the density ratio
410  const double &density_ratio() const {return *Density_ratio_pt;}
411 
412  protected:
413 
414  /// \short Helper function to adjust the position in
415  /// response to changes in position and angle of the solid
416  /// about the centre of mass
417  inline void apply_rigid_body_motion(const unsigned &t,
418  const Vector<double> &initial_x,
419  Vector<double> &r) const
420  {
421  //Scale relative to the centre of mass
422  double X = initial_x[0] - Initial_centre_of_mass[0];
423  double Y = initial_x[1] - Initial_centre_of_mass[1];
424 
425  //Find the original angle and radius
426  double phi_orig=atan2(Y,X);
427  double r_orig=sqrt(X*X+Y*Y);
428 
429  // Updated position vector
430  r[0] = Initial_centre_of_mass[0] +
432 
433  r[1] = Initial_centre_of_mass[1] +
435 
436  //Add in the rotation terms if there are to be included
438  {
439  r[0] +=
440  r_orig*cos(phi_orig + this->Centre_displacement_data_pt->value(t,2));
441  r[1] +=
442  r_orig*sin(phi_orig + this->Centre_displacement_data_pt->value(t,2));
443  }
444  else
445  {
446  r[0] += r_orig*cos(phi_orig);
447  r[1] += r_orig*sin(phi_orig);
448  }
449  }
450 
451 
452  private:
453 
454  /// \short Return the equation number associated with the i-th
455  /// centre of gravity displacment
456  /// 0: x-displ; 1: y-displ; 2: rotation angle.
457  inline int centre_displacement_local_eqn(const unsigned &i)
458  {
460  {
462  }
463  else
464  {
466  }
467  }
468 
469  /// \short Initialisation function
470  void initialise(TimeStepper* const &time_stepper_pt);
471 
472  /// Get residuals and/or Jacobian
474  DenseMatrix<double> &jacobian,
475  const bool& flag);
476 
477  /// Storage for the external data that is formed from hijacked data
478  /// that must be deleted by this element
479  std::list<unsigned> List_of_external_hijacked_data;
480 
481  /// Delete the storage for the external data formed from hijacked data
483  {
484  for(std::list<unsigned>::iterator it =
486  it!=List_of_external_hijacked_data.end();++it)
487  {
488  //Delete the associated external data
489  delete this->external_data_pt(*it);
490  this->external_data_pt(*it) =0;
491  }
492  //Now clear the list
494  }
495 
496  protected:
497 
498  /// X-coordinate of initial centre of gravity
500 
501  /// Original rotation angle
502  double Initial_Phi;
503 
504  // Mass of body
505  double Mass;
506 
507  /// Polar moment of inertia of body
509 
510  /// \short Data for centre of gravity displacement.
511  /// Values: 0: x-displ; 1: y-displ; 2: rotation angle.
513 
514  private:
515 
516  /// Underlying geometric object
518 
519  /// \short Function pointer to function that specifies
520  /// external force
522 
523  /// \short Function pointer to function that specifies
524  /// external torque
526 
527  /// \short Mesh containing face elements that allow the computation of
528  /// the drag on the body
530 
531  /// The direction of gravity
533 
534  /// Reynolds number of external fluid
535  double* Re_pt;
536 
537  /// Strouhal number of external fluid
538  double *St_pt;
539 
540  /// Reynolds number divided by Froude number of external fluid
541  double *ReInvFr_pt;
542 
543  /// Density ratio of the solid to the external fluid
545 
546  /// Static default value for physical constants
548 
549  /// Static default value for physical ratios
551 
552  /// Static default value for gravity
554 
555  /// \short Index for the data (internal or external) that contains the
556  /// centre-of-gravity displacement
558 
559  /// Boolean flag to indicate whether data is internal
561 
562  /// Boolean to indicate that the rotation variable does not affect the boundary
563  /// shape
565 
566 };
567 
568 
569 
570 //=====================================================================
571 /// Class upgrading a TriangleMeshPolygon to a "hole" for use during
572 /// triangle mesh generation. For mesh generation purposes, the main (and only)
573 /// addition to the base class is the provision of the coordinates
574 /// of a hole inside the polygon. To faciliate the movement of the "hole"
575 /// through the domain we also provide a Data object whose three values
576 /// represent the x and y displacements of its centre of gravity and
577 /// the polygon's rotation about its centre of gravity.
578 /// If added to a mesh in the Problem (in its incarnation as a
579 /// GeneralisedElement) the displacement/rotation of the polygon
580 /// is computed in response to (i) user-specifiable applied forces
581 /// and a torque and (ii) the net drag (and associated torque) from
582 /// a mesh of elements that can exert a drag onto the polygon (typically
583 /// Navier-Stokes FaceElements that apply a viscous drag to an
584 /// immersed body, represented by the polygon.)
585 //=====================================================================
587  public TriangleMeshPolygon,
589 {
590 
591 public:
592 
593 
594  /// \short Constructor: Specify coordinates of a point inside the hole
595  /// and a vector of pointers to TriangleMeshPolyLines
596  /// that define the boundary segments of the polygon.
597  /// Each TriangleMeshPolyLine has its own boundary ID and can contain
598  /// multiple (straight-line) segments. The optional final argument
599  /// is a pointer to a Data object whose three values represent
600  /// the two displacements of and the rotation angle about the polygon's
601  /// centre of mass.
604  boundary_polyline_pt,
607 
608  /// \short Empty Destuctor
610  {
611 
612  }
613 
614  ///Overload (again) the position to apply the rotation and translation
615  void position(const Vector<double> &xi, Vector<double> &r) const
616  {
617  Vector<double> initial_x(2);
618  this->get_initial_position(xi,initial_x);
619  this->apply_rigid_body_motion(0,initial_x,r);
620  }
621 
622 
623  ///Overload (again) the position to apply the rotation and translation
624  void position(const unsigned &t,
625  const Vector<double> &xi, Vector<double> &r) const
626  {
627  Vector<double> initial_x(2);
628  this->get_initial_position(xi,initial_x);
629  this->apply_rigid_body_motion(t,initial_x,r);
630  }
631 
632 
633 
634  /// \short Update the reference configuration by re-setting the original
635  /// position of the vertices to their current ones, re-set the
636  /// original position of the centre of mass, and the displacements
637  /// and rotations relative to it
639 
640 private:
641 
642  /// \short Get the initial position of the polygon
644  {
645  //Find the number of polylines (boundaries)
646  unsigned n_poly = this->npolyline();
647 
648  //The boundary coordinate will be contiguous from polyline to
649  //polyline and each polyline will have the scaled arclength coordinate
650  //in the range 0->1.
651 
652  //Find the maximum coordinate
653  double zeta_max = Zeta_vertex[n_poly-1].back();
654 
655  //If we are above the maximum complain
656  if(xi[0] > zeta_max)
657  {
658  std::ostringstream error_message;
659  error_message << "Value of intrinsic coordinate " << xi[0] <<
660  "greater than maximum " << zeta_max << "\n";
661  throw
662  OomphLibError(error_message.str(),
663  "TriangleMeshPolygon::position()",
664  OOMPH_EXCEPTION_LOCATION);
665  }
666 
667  //If equal to the maximum, return the final vertex
668  if(xi[0] == zeta_max)
669  {
670  unsigned n_vertex = this->polyline_pt(n_poly-1)->nvertex();
671  r = this->polyline_pt(n_poly-1)->vertex_coordinate(n_vertex-1);
672  return;
673  }
674 
675  //Otherwise
676 
677  //Find out which polyline we are in
678  //If we've got here this should be less than n_poly
679  unsigned p = static_cast<unsigned> (floor(xi[0]));
680 
681  //If we are "above" the last polyline then throw an error
682  //This should have been caught by the trap above
683  if(p>=n_poly)
684  {
685  std::ostringstream error_message;
686  error_message
687  << "Something has gone wrong.\n"
688  << "The integer part of the input intrinsic coordinate is "
689  << p <<
690  "\nwhich is equal to or greater than the number of polylines, "
691  << n_poly << ".\n"
692  << "This should have triggered an earlier error\n";
693 
694 
695  throw OomphLibError(error_message.str(),
696  OOMPH_CURRENT_FUNCTION,
697  OOMPH_EXCEPTION_LOCATION);
698  }
699 
700  //Cache the appropriate polyline
701  TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
702 
703  //If we are at the first vertex in the polyline, return it
704  if(xi[0] == Zeta_vertex[p][0])
705  {
706  r = line_pt->vertex_coordinate(0);
707  return;
708  }
709 
710  //Otherwise loop over the other points to find the appropriate
711  //segment
712 
713  //Find the number of vertices in the chosen polyline
714  unsigned n_vertex = line_pt->nvertex();
715  //Now start from the first node in the appropriate polyline and loop up
716  for(unsigned v=1;v<n_vertex;v++)
717  {
718  //First time that zeta falls below the vertex coordinate
719  //we have something
720  if(xi[0] < Zeta_vertex[p][v])
721  {
722  double fraction = (xi[0] - Zeta_vertex[p][v-1])/
723  (Zeta_vertex[p][v] - Zeta_vertex[p][v-1]);
724  Vector<double> first = line_pt->vertex_coordinate(v-1);
725  Vector<double> last = line_pt->vertex_coordinate(v);
726  r.resize(2);
727  for(unsigned i=0;i<2;i++)
728  {
729  r[i] = first[i] + fraction*(last[i] - first[i]);
730  }
731  return;
732  }
733  }
734  }
735 
736 
737  /// Helper function to assign the values of the (scaled) arc-length
738  /// to each node of each polyline. The direction will be the natural
739  /// order of the vertices within the polyline.
740  void assign_zeta()
741  {
742  //Find the number of polylines
743  unsigned n_poly = this->npolyline();
744 
745  //Allocate appropriate storage for the zeta values
746  Zeta_vertex.resize(n_poly);
747 
748  //Temporary storage for the vertex coordinates
749  Vector<double> vertex_coord_first;
750  Vector<double> vertex_coord_next;
751 
752  //Set the initial value of zeta
753  double zeta_offset = 0.0;
754 
755  //Loop over the polylines
756  for(unsigned p=0;p<n_poly;++p)
757  {
758  //Cache the pointer to the polyline
759  TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
760 
761  //Find the number of vertices in the polyline
762  unsigned n_vertex = line_pt->nvertex();
763 
764  //Allocate storage and set initial value
765  Zeta_vertex[p].resize(n_vertex);
766  Zeta_vertex[p][0] = 0.0;
767 
768  //Loop over the vertices in the polyline and calculate the length
769  //between each for use as the intrinsic coordinate
770  vertex_coord_first = line_pt->vertex_coordinate(0);
771  for(unsigned v=1;v<n_vertex;v++)
772  {
773  vertex_coord_next = line_pt->vertex_coordinate(v);
774  double length =
775  sqrt(pow(vertex_coord_next[0] - vertex_coord_first[0],2.0)
776  + pow(vertex_coord_next[1] - vertex_coord_first[1],2.0));
777  Zeta_vertex[p][v] = Zeta_vertex[p][v-1] + length;
778  vertex_coord_first = vertex_coord_next;
779  }
780 
781  //Now scale the length to unity and add the offset
782  Zeta_vertex[p][0] += zeta_offset;
783  for(unsigned v=1;v<n_vertex;v++)
784  {
785  Zeta_vertex[p][v] /= Zeta_vertex[p][n_vertex-1];
786  Zeta_vertex[p][v] += zeta_offset;
787  }
788  zeta_offset += 1.0;
789  } //End of loop over polylines
790  }
791 
792  /// \short Vector of intrisic coordinate values at the nodes
794 
795 };
796 
797 }
798 #endif
GeomObject * Geom_object_pt
Underlying geometric object.
A Generalised Element class.
Definition: elements.h:76
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
int centre_displacement_local_eqn(const unsigned &i)
Return the equation number associated with the i-th centre of gravity displacment 0: x-displ; 1: y-di...
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
void pin_centre_of_mass_coordinate(const unsigned &i)
Pin the i-th coordinate of the centre of mass.
Vector< double > centre_of_gravity()
Get current centre of gravity.
void delete_external_hijacked_data()
Delete the storage for the external data formed from hijacked data.
void reset_after_internal_fd()
After all internal data finite-differencing, update nodal positions.
double Initial_Phi
Original rotation angle.
double * ReInvFr_pt
Reynolds number divided by Froude number of external fluid.
Data *& centre_displacement_data_pt()
Pointer to Data for centre of gravity displacement. Values: 0: x-displ; 1: y-displ; 2: rotation angle...
void dposition_dt(const Vector< double > &zeta, const unsigned &j, Vector< double > &drdt)
Work out the position derivative, including rigid body motion.
void set_drag_mesh(Mesh *const &drag_mesh_pt)
Function to set the drag mesh and add the appropriate load and geometric data as external data to the...
Vector< double > Initial_centre_of_mass
X-coordinate of initial centre of gravity.
double & moment_of_inertia_shape()
Access to dimensionless polar "moment of inertia" shape parameter.
void reset_reference_configuration()
Update the reference configuration by re-setting the original position of the vertices to their curre...
Mesh *const & drag_mesh_pt()
Access fct to mesh containing face elements that allow the computation of the drag on the body...
void set_geometric_rotation()
Set the rotation of the object to be included.
void fill_in_jacobian_from_internal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the internal degrees of freedom using finite differe...
Definition: elements.cc:1064
void reset_in_internal_fd(const unsigned &i)
Do nothing to reset within finite-differencing of internal data.
ImmersedRigidBodyTriangleMeshPolygon(const Vector< double > &hole_center, const Vector< TriangleMeshCurveSection * > &boundary_polyline_pt, TimeStepper *const &time_stepper_pt, Data *const &centre_displacement_data_pt=0)
Constructor: Specify coordinates of a point inside the hole and a vector of pointers to TriangleMeshP...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Get the contribution to the residuals.
cstr elem_len * i
Definition: cfortran.h:607
void unpin_rotation_angle()
Unpin the rotation angle.
void position(const Vector< double > &xi, Vector< double > &r) const
Overload (again) the position to apply the rotation and translation.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Get residuals including contribution to jacobian.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Overload to include the time history of the motion of the object.
const Vector< double > & g() const
Access function for gravity.
void get_residuals_rigid_body_generic(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &flag)
Get residuals and/or Jacobian.
Vector< Vector< double > > Zeta_vertex
Vector of intrisic coordinate values at the nodes.
static double Default_Physical_Ratio_Value
Static default value for physical ratios.
void flush_drag_mesh()
Function to clear the drag mesh and all associated external data.
double & centre_rotation_angle()
rotation of centre of mass
const double & re_invfr()
Access to the fluid inverse Froude number.
char t
Definition: cfortran.h:572
void unpin_centre_of_mass_coordinate(const unsigned &i)
Unpin the i-th coordinate of the centre of mass.
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
unsigned Index_for_centre_displacement
Index for the data (internal or external) that contains the centre-of-gravity displacement.
static Vector< double > Default_Gravity_vector
Static default value for gravity.
e
Definition: cfortran.h:575
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void fill_in_jacobian_from_external_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the external degrees of freedom using finite differe...
Definition: elements.cc:1158
void flush_external_data()
Flush all external data.
Definition: elements.cc:365
double * Density_ratio_pt
Density ratio of the solid to the external fluid.
unsigned nvertex() const
Number of vertices.
const double & re() const
Access function for the fluid Reynolds number.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
const double & initial_centre_of_mass(const unsigned &i) const
Access function for the initial centre of mass (const version)
double & centre_y_displacement()
y-displacement of centre of mass
void reset_after_external_fd()
After all external data finite-differencing, update nodal positions.
void node_update_adjacent_fluid_elements()
Update the positions of the nodes in fluid elements adjacent to the rigid body, defined as being elem...
void get_initial_position(const Vector< double > &xi, Vector< double > &r) const
Get the initial position of the polygon.
void get_force_and_torque(const double &time, Vector< double > &force, double &torque)
Get force and torque from specified fct pointers and drag mesh.
double & centre_x_displacement()
x-displacement of centre of mass
Class defining a polyline for use in Triangle Mesh generation.
double * St_pt
Strouhal number of external fluid.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
void unset_geometric_rotation()
Set the rotation of the object to be ignored (only really useful if you have a circular shape) ...
double Moment_of_inertia
Polar moment of inertia of body.
const double & st() const
Access function for the fluid Strouhal number.
Vector< double > vertex_coordinate(const unsigned &i) const
Coordinate vector of i-th vertex (const version)
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void update_in_internal_fd(const unsigned &i)
After an internal data change, update the nodal positions.
ExternalForceFctPt External_force_fct_pt
Function pointer to function that specifies external force.
Data * Centre_displacement_data_pt
Data for centre of gravity displacement. Values: 0: x-displ; 1: y-displ; 2: rotation angle...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
ExternalTorqueFctPt & external_torque_fct_pt()
Access to function pointer to function that specifies external torque.
ExternalTorqueFctPt External_torque_fct_pt
Function pointer to function that specifies external torque.
ExternalForceFctPt & external_force_fct_pt()
Access to function pointer to function that specifies external force.
double *& density_ratio_pt()
Access function for the pointer to the density ratio.
ImmersedRigidBodyElement(GeomObject *const &geom_object_pt, TimeStepper *const &time_stepper_pt, Data *const &centre_displacement_data_pt=0)
Constructor that takes an underlying geometric object: and timestepper.
void(* ExternalTorqueFctPt)(const double &time, double &external_torque)
Function pointer to function that specifies external torque.
ImmersedRigidBodyElement(TimeStepper *const &time_stepper_pt, Data *const &centre_displacement_data_pt=0)
Default constructor that intialises everything to zero. This is expected to be called only from deriv...
unsigned npolyline() const
Number of constituent polylines.
static double Default_Physical_Constant_Value
Static default value for physical constants.
bool Displacement_data_is_internal
Boolean flag to indicate whether data is internal.
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:662
void initialise(TimeStepper *const &time_stepper_pt)
Initialisation function.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
Vector< double > *& g_pt()
Access function to the direction of gravity.
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th (only) Data item that the object's shape depends on.
double *& re_invfr_pt()
Access function for pointer to the fluid inverse Froude number (dimensionless gravitational loading) ...
double *& st_pt()
Access function for the pointer to the fluid Strouhal number.
void reset_in_external_fd(const unsigned &i)
Do nothing to reset within finite-differencing of external data.
const double & density_ratio() const
Access function for the the density ratio.
Mesh * Drag_mesh_pt
Mesh containing face elements that allow the computation of the drag on the body. ...
void apply_rigid_body_motion(const unsigned &t, const Vector< double > &initial_x, Vector< double > &r) const
Helper function to adjust the position in response to changes in position and angle of the solid abou...
unsigned ngeom_data() const
The position of the object depends on one data item.
double & initial_centre_of_mass(const unsigned &i)
Access function for the initial centre of mass.
void pin_rotation_angle()
Pin the rotation angle.
double & initial_phi()
Access function for the initial angle.
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
void output_centre_of_gravity(std::ostream &outfile)
Output position velocity and acceleration of centre of gravity.
void(* ExternalForceFctPt)(const double &time, Vector< double > &external_force)
Function pointer to function that specifies external force.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
Vector< double > * G_pt
The direction of gravity.
~ImmersedRigidBodyElement()
Destuctor: Cleanup if required.
double *& re_pt()
Access function for the pointer to the fluid Reynolds number.
A general mesh class.
Definition: mesh.h:74
void update_in_external_fd(const unsigned &i)
After an external data change, update the nodal positions.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Overload (again) the position to apply the rotation and translation.
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data...
Definition: elements.h:313
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
double * Re_pt
Reynolds number of external fluid.
void position(const Vector< double > &xi, Vector< double > &r) const
Overload the position to apply the rotation and translation.
virtual void node_update()
Update the positions of all nodes in the element using each node update function. The default impleme...
Definition: elements.cc:4865
Class defining a closed polygon for the Triangle mesh generation.