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// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for Rigid Body Elements that are immersed in an external fluid
27 #ifndef OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
28 #define OOMPH_IMMERSED_RIGID_BODY_ELEMENTS_HEADER
29 
30 
31 // Config header generated by autoconfig
32 #ifdef HAVE_CONFIG_H
33 #include <oomph-lib-config.h>
34 #endif
35 
36 #include "../generic/elements.h"
37 #include "../generic/triangle_mesh.h"
38 #include "../generic/fsi.h"
39 
40 namespace oomph
41 {
42  //=====================================================================
43  // DOXYERROR: I removed a \f from the formula \f$\mbox{\boldmath$F$} = \mbox{\boldmath$F$}^{*}/\mu U\f$ below because it broke doxygen. Someone
44  // who knows the maths should check that this still makes sense.
45  /// Class that solves the equations of motion for a general
46  /// two-dimensional rigid body subject to a particular imposed force and
47  /// torque distribution and immersed within an external fluid. The body's
48  /// position is entirely specified by the location of its centre of mass,
49  /// \f$\mbox{\boldmath$X$}\f$, and a single angle, \f$\phi\f$, that represents
50  /// a possible rotation. The equations of motion are then simply Newton's
51  /// second law for the conservation of linear momentum in two directions and
52  /// angular momentum about the single possible axis of rotation.
53  ///
54  /// The non-dimensionalisation is based on the viscous scales of
55  /// the surrounding fluid, in which case, the governing equations are
56  /// \f[ Re St^{2} M \frac{\mbox{d}^{2} \mbox{\boldmath$X$}}{\mbox{d} t^{2}} = \mbox{\boldmath$F$} + \oint \mbox{\boldmath$t$} \mbox{d} s \quad\mbox{and}\quad Re St^{2} I \frac{\mbox{d}^{2}\Phi}{\mbox{d} t^{2}} = \mbox{\boldmath$T$} + \oint \mbox{\boldmath$q$} \mbox{d} s, \f]
57  /// where \f$\mbox{\boldmath$F$} = \mbox{\boldmath$F$}^{*}/\mu U\f$ is the
58  /// external force per unit length;
59  /// \f$\mbox{\boldmath$t$}\f$ is the net force per unit length
60  /// applied on the rigid body
61  /// by the surrounding fluid;
62  /// \f$\mbox{\boldmath$T$} = \mbox{\boldmath$T$}^{*}/(\mu UL)\f$
63  /// is the external torque per unit length; and \f$\mbox{\boldmath$q$}\f$ is
64  /// the net torque per unit length applied on the body by the fluid.
65  /// \f$M\f$ is a scaled mass the density ratio multiplied by a shape parameter
66  /// and \f$I\f$ is a scaled moment of inertia. Finally, \f$Re\f$ and \f$St\f$
67  /// are the Reynolds and Strouhal numbers of the surrounding fluid.
68  /// Note that these equations may be used without the external fluid,
69  /// in which case the non-dimensionalisation doesn't make a lot of sense,
70  /// but may be re-interpreted accordingly.
71  ///
72  /// A Data object whose three values
73  /// represent the x and y displacements of the body's centre of mass and
74  /// its rotation about the centre of mass may be passed in as external
75  /// data or, if not, it will be constructed internally.
76  ///
77  /// For general usage, an underlying geometric object must passed to the
78  /// constructor and the position will be determined by applying rigid body
79  /// motions to the underlying object based on the initial_centre_of_mass
80  /// and initial_phi (angle), which default to zero. If these defaults are
81  /// not suitable, the values must be set externally. In addition a mass
82  /// and moment of inertia should also be set externally.
83  ///
84  /// If added to a mesh in the Problem (in its incarnation as a
85  /// GeneralisedElement) the displacement/rotation of the body
86  /// is computed in response to (i) user-specifiable applied forces
87  /// and a torque and (ii) the net drag (and associated torque) from
88  /// a mesh of elements that can exert a drag onto the body (typically
89  /// Navier-Stokes FaceElements that apply a viscous drag to an
90  /// immersed body, represented by the body.)
91  //=====================================================================
93  {
94  public:
95  /// Function pointer to function that specifies
96  /// external force
97  typedef void (*ExternalForceFctPt)(const double& time,
98  Vector<double>& external_force);
99 
100  /// Function pointer to function that specifies
101  /// external torque
102  typedef void (*ExternalTorqueFctPt)(const double& time,
103  double& external_torque);
104 
105  protected:
106  /// Default constructor that intialises everything to
107  /// zero. This is expected to be called only from derived clases
108  /// such as the ImmersedRigidBodyTriangleMeshPolygon that can provided
109  /// their own position() functions.
111  Data* const& centre_displacement_data_pt = 0)
112  : Initial_Phi(0.0),
113  Mass(0.0),
114  Moment_of_inertia(0.0),
116  Geom_object_pt(0),
119  Drag_mesh_pt(0),
126  {
127  this->initialise(time_stepper_pt);
128  }
129 
130  public:
131  /// Constructor that takes an underlying geometric object:
132  /// and timestepper.
133  ImmersedRigidBodyElement(GeomObject* const& geom_object_pt,
135  Data* const& centre_displacement_data_pt = 0)
136  : Initial_Phi(0.0),
137  Mass(0.0),
138  Moment_of_inertia(0.0),
140  Geom_object_pt(geom_object_pt),
143  Drag_mesh_pt(0),
150  {
151  this->initialise(time_stepper_pt);
152  }
153 
154  /// Set the rotation of the object to be included
156  {
158  }
159 
160  /// Set the rotation of the object to be ignored (only really
161  /// useful if you have a circular shape)
163  {
165  }
166 
167  /// Access function for the initial angle
168  double& initial_phi()
169  {
170  return Initial_Phi;
171  }
172 
173  /// Access function for the initial centre of mass
174  double& initial_centre_of_mass(const unsigned& i)
175  {
176  return Initial_centre_of_mass[i];
177  }
178 
179  /// Access function for the initial centre of mass (const version)
180  const double& initial_centre_of_mass(const unsigned& i) const
181  {
182  return Initial_centre_of_mass[i];
183  }
184 
185  /// Overload the position to apply the rotation and translation
186  void position(const Vector<double>& xi, Vector<double>& r) const
187  {
188  Vector<double> initial_x(2);
189  Geom_object_pt->position(xi, initial_x);
190  this->apply_rigid_body_motion(0, initial_x, r);
191  }
192 
193  /// Overload to include the time history of the motion of the object
194  void position(const unsigned& t,
195  const Vector<double>& xi,
196  Vector<double>& r) const
197  {
198  Vector<double> initial_x(2);
199  Geom_object_pt->position(xi, initial_x);
200  this->apply_rigid_body_motion(t, initial_x, r);
201  }
202 
203 
204  /// Work out the position derivative, including rigid body motion
205  void dposition_dt(const Vector<double>& zeta,
206  const unsigned& j,
207  Vector<double>& drdt);
208 
209  /// Destuctor: Cleanup if required
211  {
213  {
215  // Null out the pointer stored as internal data
216  this->internal_data_pt(0) = 0;
217  }
218  }
219 
220  /// Access to dimensionless "mass" shape parameter that must be set by hand
221  /// for non polygonal shapes
222  double& mass_shape()
223  {
224  return Mass;
225  }
226 
227  /// Access to dimensionless polar "moment of inertia" shape parameter
229  {
230  return Moment_of_inertia;
231  }
232 
233  /// Pointer to Data for centre of gravity displacement.
234  /// Values: 0: x-displ; 1: y-displ; 2: rotation angle.
236  {
238  }
239 
240  /// x-displacement of centre of mass
242  {
244  }
245 
246  /// y-displacement of centre of mass
248  {
250  }
251 
252  /// rotation of centre of mass
254  {
256  }
257 
258  /// Get current centre of gravity
260  {
261  Vector<double> cog(2);
262  for (unsigned i = 0; i < 2; i++)
263  {
264  cog[i] = Initial_centre_of_mass[i] +
266  }
267  return cog;
268  }
269 
270  /// Pin the i-th coordinate of the centre of mass
271  void pin_centre_of_mass_coordinate(const unsigned& i)
272  {
274  }
275 
276  /// Unpin the i-th coordinate of the centre of mass
277  void unpin_centre_of_mass_coordinate(const unsigned& i)
278  {
280  }
281 
282  /// Pin the rotation angle
284  {
286  }
287 
288  /// Unpin the rotation angle
290  {
292  }
293 
294  /// Output position velocity and acceleration of centre of gravity
295  void output_centre_of_gravity(std::ostream& outfile);
296 
297  /// Get the contribution to the residuals
299  {
300  // Get generic function without jacobian terms
302  residuals, GeneralisedElement::Dummy_matrix, false);
303  }
304 
305 
306  /// Get residuals including contribution to jacobian
308  DenseMatrix<double>& jacobian)
309  {
310  // Get generic function, but don't bother to get jacobian terms
311  bool flag = false;
312  get_residuals_rigid_body_generic(residuals, jacobian, flag);
313  // Get the internal effects by FD, because a change in internal
314  // data can change the fluid loading, so this is one of those
315  // subtle interactions
316  this->fill_in_jacobian_from_internal_by_fd(residuals, jacobian);
317  // Get the effect of the fluid loading on the rigid body
318  this->fill_in_jacobian_from_external_by_fd(residuals, jacobian);
319  }
320 
321  /// Update the positions of the nodes in fluid elements
322  /// adjacent to the rigid body, defined as being elements in the
323  /// drag mesh.
325  {
326  // If there is no drag mesh then do nothing
327  if (Drag_mesh_pt == 0)
328  {
329  return;
330  }
331  // Otherwise update the elements adjacent to the rigid body
332  //(all the elements adjacent to those in the drag mesh)
333  else
334  {
335  unsigned nel = Drag_mesh_pt->nelement();
336  for (unsigned e = 0; e < nel; e++)
337  {
338  dynamic_cast<FaceElement*>(Drag_mesh_pt->element_pt(e))
339  ->bulk_element_pt()
340  ->node_update();
341  }
342  }
343  }
344 
345 
346  /// After an external data change, update the nodal positions
347  inline void update_in_external_fd(const unsigned& i)
348  {
350  }
351 
352  /// Do nothing to reset within finite-differencing of external data
353  inline void reset_in_external_fd(const unsigned& i) {}
354 
355  /// After all external data finite-differencing, update nodal
356  /// positions
358  {
360  }
361 
362  /// After an internal data change, update the nodal positions
363  inline void update_in_internal_fd(const unsigned& i)
364  {
366  }
367 
368  /// Do nothing to reset within finite-differencing of internal data
369  inline void reset_in_internal_fd(const unsigned& i) {}
370 
371  /// After all internal data finite-differencing, update nodal
372  /// positions
374  {
376  }
377 
378  /// Get force and torque from specified fct pointers and
379  /// drag mesh
380  void get_force_and_torque(const double& time,
381  Vector<double>& force,
382  double& torque);
383 
384  /// Access to function pointer to function that specifies
385  /// external force
387  {
388  return External_force_fct_pt;
389  }
390 
391  /// Access to function pointer to function that specifies
392  /// external torque
394  {
395  return External_torque_fct_pt;
396  }
397 
398  /// Access fct to mesh containing face elements that allow
399  /// the computation of the drag on the body
400  Mesh* const& drag_mesh_pt()
401  {
402  return Drag_mesh_pt;
403  }
404 
405  /// Function to set the drag mesh and add the appropriate load
406  /// and geometric data as external data to the Rigid Body
407  void set_drag_mesh(Mesh* const& drag_mesh_pt);
408 
409  /// Function to clear the drag mesh and all associated external data
411  {
412  // Delete the hijacked data created as the load
414  // Flush the external data
415  this->flush_external_data();
416  // Set the Drag_mesh pointer to null
417  Drag_mesh_pt = 0;
418  }
419 
420  /// The position of the object depends on one data item
421  unsigned ngeom_data() const
422  {
423  return 1;
424  }
425 
426  /// Return pointer to the j-th (only) Data item that the object's
427  /// shape depends on.
428  Data* geom_data_pt(const unsigned& j)
429  {
430  return this->Centre_displacement_data_pt;
431  }
432 
433  /// Access function to the direction of gravity
435  {
436  return G_pt;
437  }
438 
439  /// Access function for gravity
440  const Vector<double>& g() const
441  {
442  return *G_pt;
443  }
444 
445  /// Access function for the pointer to the fluid Reynolds number
446  double*& re_pt()
447  {
448  return Re_pt;
449  }
450 
451  /// Access function for the fluid Reynolds number
452  const double& re() const
453  {
454  return *Re_pt;
455  }
456 
457  /// Access function for the pointer to the fluid Strouhal number
458  double*& st_pt()
459  {
460  return St_pt;
461  }
462 
463  /// Access function for the fluid Strouhal number
464  const double& st() const
465  {
466  return *St_pt;
467  }
468 
469  /// Access function for pointer to the fluid inverse Froude number
470  /// (dimensionless gravitational loading)
471  double*& re_invfr_pt()
472  {
473  return ReInvFr_pt;
474  }
475 
476  /// Access to the fluid inverse Froude number
477  const double& re_invfr()
478  {
479  return *ReInvFr_pt;
480  }
481 
482  /// Access function for the pointer to the density ratio
483  double*& density_ratio_pt()
484  {
485  return Density_ratio_pt;
486  }
487 
488  /// Access function for the the density ratio
489  const double& density_ratio() const
490  {
491  return *Density_ratio_pt;
492  }
493 
494  protected:
495  /// Helper function to adjust the position in
496  /// response to changes in position and angle of the solid
497  /// about the centre of mass
498  inline void apply_rigid_body_motion(const unsigned& t,
499  const Vector<double>& initial_x,
500  Vector<double>& r) const
501  {
502  // Scale relative to the centre of mass
503  double X = initial_x[0] - Initial_centre_of_mass[0];
504  double Y = initial_x[1] - Initial_centre_of_mass[1];
505 
506  // Find the original angle and radius
507  double phi_orig = atan2(Y, X);
508  double r_orig = sqrt(X * X + Y * Y);
509 
510  // Updated position vector
511  r[0] = Initial_centre_of_mass[0] +
512  this->Centre_displacement_data_pt->value(t, 0);
513 
514  r[1] = Initial_centre_of_mass[1] +
515  this->Centre_displacement_data_pt->value(t, 1);
516 
517  // Add in the rotation terms if there are to be included
519  {
520  r[0] += r_orig *
521  cos(phi_orig + this->Centre_displacement_data_pt->value(t, 2));
522  r[1] += r_orig *
523  sin(phi_orig + this->Centre_displacement_data_pt->value(t, 2));
524  }
525  else
526  {
527  r[0] += r_orig * cos(phi_orig);
528  r[1] += r_orig * sin(phi_orig);
529  }
530  }
531 
532 
533  private:
534  /// Return the equation number associated with the i-th
535  /// centre of gravity displacment
536  /// 0: x-displ; 1: y-displ; 2: rotation angle.
537  inline int centre_displacement_local_eqn(const unsigned& i)
538  {
540  {
542  }
543  else
544  {
546  }
547  }
548 
549  /// Initialisation function
550  void initialise(TimeStepper* const& time_stepper_pt);
551 
552  /// Get residuals and/or Jacobian
554  DenseMatrix<double>& jacobian,
555  const bool& flag);
556 
557  /// Storage for the external data that is formed from hijacked data
558  /// that must be deleted by this element
559  std::list<unsigned> List_of_external_hijacked_data;
560 
561  /// Delete the storage for the external data formed from hijacked data
563  {
564  for (std::list<unsigned>::iterator it =
566  it != List_of_external_hijacked_data.end();
567  ++it)
568  {
569  // Delete the associated external data
570  delete this->external_data_pt(*it);
571  this->external_data_pt(*it) = 0;
572  }
573  // Now clear the list
575  }
576 
577  protected:
578  /// X-coordinate of initial centre of gravity
580 
581  /// Original rotation angle
582  double Initial_Phi;
583 
584  // Mass of body
585  double Mass;
586 
587  /// Polar moment of inertia of body
589 
590  /// Data for centre of gravity displacement.
591  /// Values: 0: x-displ; 1: y-displ; 2: rotation angle.
593 
594  private:
595  /// Underlying geometric object
597 
598  /// Function pointer to function that specifies
599  /// external force
601 
602  /// Function pointer to function that specifies
603  /// external torque
605 
606  /// Mesh containing face elements that allow the computation of
607  /// the drag on the body
609 
610  /// The direction of gravity
612 
613  /// Reynolds number of external fluid
614  double* Re_pt;
615 
616  /// Strouhal number of external fluid
617  double* St_pt;
618 
619  /// Reynolds number divided by Froude number of external fluid
620  double* ReInvFr_pt;
621 
622  /// Density ratio of the solid to the external fluid
624 
625  /// Static default value for physical constants
627 
628  /// Static default value for physical ratios
630 
631  /// Static default value for gravity
633 
634  /// Index for the data (internal or external) that contains the
635  /// centre-of-gravity displacement
637 
638  /// Boolean flag to indicate whether data is internal
640 
641  /// Boolean to indicate that the rotation variable does not affect the
642  /// boundary shape
644  };
645 
646 
647  //=====================================================================
648  /// Class upgrading a TriangleMeshPolygon to a "hole" for use during
649  /// triangle mesh generation. For mesh generation purposes, the main (and
650  /// only) addition to the base class is the provision of the coordinates of a
651  /// hole inside the polygon. To faciliate the movement of the "hole" through
652  /// the domain we also provide a Data object whose three values represent the
653  /// x and y displacements of its centre of gravity and the polygon's rotation
654  /// about its centre of gravity. If added to a mesh in the Problem (in its
655  /// incarnation as a GeneralisedElement) the displacement/rotation of the
656  /// polygon is computed in response to (i) user-specifiable applied forces and
657  /// a torque and (ii) the net drag (and associated torque) from a mesh of
658  /// elements that can exert a drag onto the polygon (typically Navier-Stokes
659  /// FaceElements that apply a viscous drag to an immersed body, represented by
660  /// the polygon.)
661  //=====================================================================
664  {
665  public:
666  /// Constructor: Specify coordinates of a point inside the hole
667  /// and a vector of pointers to TriangleMeshPolyLines
668  /// that define the boundary segments of the polygon.
669  /// Each TriangleMeshPolyLine has its own boundary ID and can contain
670  /// multiple (straight-line) segments. The optional final argument
671  /// is a pointer to a Data object whose three values represent
672  /// the two displacements of and the rotation angle about the polygon's
673  /// centre of mass.
675  const Vector<double>& hole_center,
676  const Vector<TriangleMeshCurveSection*>& boundary_polyline_pt,
678  Data* const& centre_displacement_data_pt = 0);
679 
680  /// Empty Destuctor
682 
683  /// Overload (again) the position to apply the rotation and translation
684  void position(const Vector<double>& xi, Vector<double>& r) const
685  {
686  Vector<double> initial_x(2);
687  this->get_initial_position(xi, initial_x);
688  this->apply_rigid_body_motion(0, initial_x, r);
689  }
690 
691 
692  /// Overload (again) the position to apply the rotation and translation
693  void position(const unsigned& t,
694  const Vector<double>& xi,
695  Vector<double>& r) const
696  {
697  Vector<double> initial_x(2);
698  this->get_initial_position(xi, initial_x);
699  this->apply_rigid_body_motion(t, initial_x, r);
700  }
701 
702 
703  /// Update the reference configuration by re-setting the original
704  /// position of the vertices to their current ones, re-set the
705  /// original position of the centre of mass, and the displacements
706  /// and rotations relative to it
708 
709  private:
710  /// Get the initial position of the polygon
712  {
713  // Find the number of polylines (boundaries)
714  unsigned n_poly = this->npolyline();
715 
716  // The boundary coordinate will be contiguous from polyline to
717  // polyline and each polyline will have the scaled arclength coordinate
718  // in the range 0->1.
719 
720  // Find the maximum coordinate
721  double zeta_max = Zeta_vertex[n_poly - 1].back();
722 
723  // If we are above the maximum complain
724  if (xi[0] > zeta_max)
725  {
726  std::ostringstream error_message;
727  error_message << "Value of intrinsic coordinate " << xi[0]
728  << "greater than maximum " << zeta_max << "\n";
729  throw OomphLibError(error_message.str(),
730  "TriangleMeshPolygon::position()",
731  OOMPH_EXCEPTION_LOCATION);
732  }
733 
734  // If equal to the maximum, return the final vertex
735  if (xi[0] == zeta_max)
736  {
737  unsigned n_vertex = this->polyline_pt(n_poly - 1)->nvertex();
738  r = this->polyline_pt(n_poly - 1)->vertex_coordinate(n_vertex - 1);
739  return;
740  }
741 
742  // Otherwise
743 
744  // Find out which polyline we are in
745  // If we've got here this should be less than n_poly
746  unsigned p = static_cast<unsigned>(floor(xi[0]));
747 
748  // If we are "above" the last polyline then throw an error
749  // This should have been caught by the trap above
750  if (p >= n_poly)
751  {
752  std::ostringstream error_message;
753  error_message
754  << "Something has gone wrong.\n"
755  << "The integer part of the input intrinsic coordinate is " << p
756  << "\nwhich is equal to or greater than the number of polylines, "
757  << n_poly << ".\n"
758  << "This should have triggered an earlier error\n";
759 
760 
761  throw OomphLibError(error_message.str(),
762  OOMPH_CURRENT_FUNCTION,
763  OOMPH_EXCEPTION_LOCATION);
764  }
765 
766  // Cache the appropriate polyline
767  TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
768 
769  // If we are at the first vertex in the polyline, return it
770  if (xi[0] == Zeta_vertex[p][0])
771  {
772  r = line_pt->vertex_coordinate(0);
773  return;
774  }
775 
776  // Otherwise loop over the other points to find the appropriate
777  // segment
778 
779  // Find the number of vertices in the chosen polyline
780  unsigned n_vertex = line_pt->nvertex();
781  // Now start from the first node in the appropriate polyline and loop up
782  for (unsigned v = 1; v < n_vertex; v++)
783  {
784  // First time that zeta falls below the vertex coordinate
785  // we have something
786  if (xi[0] < Zeta_vertex[p][v])
787  {
788  double fraction = (xi[0] - Zeta_vertex[p][v - 1]) /
789  (Zeta_vertex[p][v] - Zeta_vertex[p][v - 1]);
790  Vector<double> first = line_pt->vertex_coordinate(v - 1);
791  Vector<double> last = line_pt->vertex_coordinate(v);
792  r.resize(2);
793  for (unsigned i = 0; i < 2; i++)
794  {
795  r[i] = first[i] + fraction * (last[i] - first[i]);
796  }
797  return;
798  }
799  }
800  }
801 
802 
803  /// Helper function to assign the values of the (scaled) arc-length
804  /// to each node of each polyline. The direction will be the natural
805  /// order of the vertices within the polyline.
806  void assign_zeta()
807  {
808  // Find the number of polylines
809  unsigned n_poly = this->npolyline();
810 
811  // Allocate appropriate storage for the zeta values
812  Zeta_vertex.resize(n_poly);
813 
814  // Temporary storage for the vertex coordinates
815  Vector<double> vertex_coord_first;
816  Vector<double> vertex_coord_next;
817 
818  // Set the initial value of zeta
819  double zeta_offset = 0.0;
820 
821  // Loop over the polylines
822  for (unsigned p = 0; p < n_poly; ++p)
823  {
824  // Cache the pointer to the polyline
825  TriangleMeshPolyLine* const line_pt = this->polyline_pt(p);
826 
827  // Find the number of vertices in the polyline
828  unsigned n_vertex = line_pt->nvertex();
829 
830  // Allocate storage and set initial value
831  Zeta_vertex[p].resize(n_vertex);
832  Zeta_vertex[p][0] = 0.0;
833 
834  // Loop over the vertices in the polyline and calculate the length
835  // between each for use as the intrinsic coordinate
836  vertex_coord_first = line_pt->vertex_coordinate(0);
837  for (unsigned v = 1; v < n_vertex; v++)
838  {
839  vertex_coord_next = line_pt->vertex_coordinate(v);
840  double length =
841  sqrt(pow(vertex_coord_next[0] - vertex_coord_first[0], 2.0) +
842  pow(vertex_coord_next[1] - vertex_coord_first[1], 2.0));
843  Zeta_vertex[p][v] = Zeta_vertex[p][v - 1] + length;
844  vertex_coord_first = vertex_coord_next;
845  }
846 
847  // Now scale the length to unity and add the offset
848  Zeta_vertex[p][0] += zeta_offset;
849  for (unsigned v = 1; v < n_vertex; v++)
850  {
851  Zeta_vertex[p][v] /= Zeta_vertex[p][n_vertex - 1];
852  Zeta_vertex[p][v] += zeta_offset;
853  }
854  zeta_offset += 1.0;
855  } // End of loop over polylines
856  }
857 
858  /// Vector of intrisic coordinate values at the nodes
860  };
861 
862 } // namespace oomph
863 #endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
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:324
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:293
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4342
virtual void node_update()
Update the positions of all nodes in the element using each node update function. The default impleme...
Definition: elements.cc:5102
A Generalised Element class.
Definition: elements.h:73
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:1227
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:1130
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:659
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
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:227
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:267
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:311
void flush_external_data()
Flush all external data.
Definition: elements.cc:392
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
Class that solves the equations of motion for a general two-dimensional rigid body subject to a parti...
const Vector< double > & g() const
Access function for gravity.
ExternalTorqueFctPt & external_torque_fct_pt()
Access to 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...
Vector< double > centre_of_gravity()
Get current centre of gravity.
double & mass_shape()
Access to dimensionless "mass" shape parameter that must be set by hand for non polygonal shapes.
bool Include_geometric_rotation
Boolean to indicate that the rotation variable does not affect the boundary shape.
void position(const Vector< double > &xi, Vector< double > &r) const
Overload the position to apply the rotation and translation.
double * St_pt
Strouhal number of external fluid.
Vector< double > *& g_pt()
Access function to the direction of gravity.
void pin_centre_of_mass_coordinate(const unsigned &i)
Pin the i-th coordinate of the centre of mass.
static double Default_Physical_Constant_Value
Static default value for physical constants.
void flush_drag_mesh()
Function to clear the drag mesh and all associated external data.
void initialise(TimeStepper *const &time_stepper_pt)
Initialisation function.
void dposition_dt(const Vector< double > &zeta, const unsigned &j, Vector< double > &drdt)
Work out the position derivative, including rigid body motion.
double & initial_centre_of_mass(const unsigned &i)
Access function for the initial centre of mass.
void reset_after_internal_fd()
After all internal data finite-differencing, update nodal positions.
ExternalTorqueFctPt External_torque_fct_pt
Function pointer to function that specifies external torque.
void(* ExternalTorqueFctPt)(const double &time, double &external_torque)
Function pointer to function that specifies external torque.
double *& re_invfr_pt()
Access function for pointer to the fluid inverse Froude number (dimensionless gravitational loading)
unsigned ngeom_data() const
The position of the object depends on one data item.
void update_in_external_fd(const unsigned &i)
After an external data change, update the nodal positions.
void reset_in_external_fd(const unsigned &i)
Do nothing to reset within finite-differencing of external data.
double * ReInvFr_pt
Reynolds number divided by Froude number of external fluid.
void delete_external_hijacked_data()
Delete the storage for the external data formed from hijacked data.
~ImmersedRigidBodyElement()
Destuctor: Cleanup if required.
double & centre_y_displacement()
y-displacement of centre of mass
double *& re_pt()
Access function for the pointer to the fluid Reynolds number.
double Initial_Phi
Original rotation angle.
void reset_in_internal_fd(const unsigned &i)
Do nothing to reset within finite-differencing of internal data.
bool Displacement_data_is_internal
Boolean flag to indicate whether data is internal.
void reset_after_external_fd()
After all external data finite-differencing, update nodal positions.
void get_residuals_rigid_body_generic(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &flag)
Get residuals and/or Jacobian.
void(* ExternalForceFctPt)(const double &time, Vector< double > &external_force)
Function pointer to function that specifies external force.
void output_centre_of_gravity(std::ostream &outfile)
Output position velocity and acceleration of centre of gravity.
const double & st() const
Access function for the fluid Strouhal number.
double *& density_ratio_pt()
Access function for the pointer to the density ratio.
Data * Centre_displacement_data_pt
Data for centre of gravity displacement. Values: 0: x-displ; 1: y-displ; 2: rotation angle.
void unpin_rotation_angle()
Unpin the rotation angle.
unsigned Index_for_centre_displacement
Index for the data (internal or external) that contains the centre-of-gravity displacement.
double * Re_pt
Reynolds number of external fluid.
static Vector< double > Default_Gravity_vector
Static default value for gravity.
void update_in_internal_fd(const unsigned &i)
After an internal data change, update the nodal positions.
GeomObject * Geom_object_pt
Underlying geometric object.
void pin_rotation_angle()
Pin the rotation angle.
double & initial_phi()
Access function for the initial angle.
ExternalForceFctPt External_force_fct_pt
Function pointer to function that specifies external force.
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...
const double & re() const
Access function for the fluid Reynolds number.
double Moment_of_inertia
Polar moment of inertia of body.
Mesh *const & drag_mesh_pt()
Access fct to mesh containing face elements that allow the computation of the drag on the body.
ExternalForceFctPt & external_force_fct_pt()
Access to function pointer to function that specifies external force.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Get the contribution to the residuals.
void unset_geometric_rotation()
Set the rotation of the object to be ignored (only really useful if you have a circular shape)
Data *& centre_displacement_data_pt()
Pointer to Data for centre of gravity displacement. Values: 0: x-displ; 1: y-displ; 2: rotation angle...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Get residuals including contribution to jacobian.
double & moment_of_inertia_shape()
Access to dimensionless polar "moment of inertia" shape parameter.
double & centre_rotation_angle()
rotation of centre of mass
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 double & re_invfr()
Access to the fluid inverse Froude number.
Vector< double > * G_pt
The direction of gravity.
double & centre_x_displacement()
x-displacement of centre of mass
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 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...
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th (only) Data item that the object's shape depends on.
Mesh * Drag_mesh_pt
Mesh containing face elements that allow the computation of the drag on the body.
void get_force_and_torque(const double &time, Vector< double > &force, double &torque)
Get force and torque from specified fct pointers and drag mesh.
const double & density_ratio() const
Access function for the the density ratio.
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...
const double & initial_centre_of_mass(const unsigned &i) const
Access function for the initial centre of mass (const version)
double *& st_pt()
Access function for the pointer to the fluid Strouhal number.
static double Default_Physical_Ratio_Value
Static default value for physical ratios.
Vector< double > Initial_centre_of_mass
X-coordinate of initial centre of gravity.
void node_update_adjacent_fluid_elements()
Update the positions of the nodes in fluid elements adjacent to the rigid body, defined as being elem...
std::list< unsigned > List_of_external_hijacked_data
Storage for the external data that is formed from hijacked data that must be deleted by this element.
double * Density_ratio_pt
Density ratio of the solid to the external fluid.
void set_geometric_rotation()
Set the rotation of the object to be included.
void unpin_centre_of_mass_coordinate(const unsigned &i)
Unpin the i-th coordinate of the centre of mass.
Class upgrading a TriangleMeshPolygon to a "hole" for use during triangle mesh generation....
void assign_zeta()
Helper function to assign the values of the (scaled) arc-length to each node of each polyline....
void position(const Vector< double > &xi, Vector< double > &r) const
Overload (again) the position to apply the rotation and translation.
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...
Vector< Vector< double > > Zeta_vertex
Vector of intrisic coordinate values at the nodes.
void reset_reference_configuration()
Update the reference configuration by re-setting the original position of the vertices to their curre...
void get_initial_position(const Vector< double > &xi, Vector< double > &r) const
Get the initial position of the polygon.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Overload (again) the position to apply the rotation and translation.
A general mesh class.
Definition: mesh.h:67
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
Class defining a polyline for use in Triangle Mesh generation.
unsigned nvertex() const
Number of vertices.
Vector< double > vertex_coordinate(const unsigned &i) const
Coordinate vector of i-th vertex (const version)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
unsigned npolyline() const
Number of constituent polylines.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...