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-2022 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
40namespace 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.
112 : Initial_Phi(0.0),
113 Mass(0.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,
136 : Initial_Phi(0.0),
137 Mass(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 {
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 {
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 {
389 }
390
391 /// Access to function pointer to function that specifies
392 /// external torque
394 {
396 }
397
398 /// Access fct to mesh containing face elements that allow
399 /// the computation of the drag on the body
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
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] +
513
514 r[1] = Initial_centre_of_mass[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
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
560
561 /// Delete the storage for the external data formed from hijacked data
563 {
564 for (std::list<unsigned>::iterator it =
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
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,
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.
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(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
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
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
virtual void node_update()
Update the positions of all nodes in the element using each node update function. The default impleme...
Definition: elements.cc:5072
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:1199
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:1102
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:387
/////////////////////////////////////////////////////////////////////
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...
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...
double & centre_y_displacement()
y-displacement of centre of mass
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.
const double & density_ratio() const
Access function for the the density ratio.
Mesh *const & drag_mesh_pt()
Access fct to mesh containing face elements that allow the computation of the drag on the body.
double & initial_phi()
Access function for the initial angle.
const double & st() const
Access function for the fluid Strouhal number.
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.
double & moment_of_inertia_shape()
Access to dimensionless polar "moment of inertia" shape parameter.
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.
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.
const Vector< double > & g() const
Access function for gravity.
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.
double & mass_shape()
Access to dimensionless "mass" shape parameter that must be set by hand for non polygonal shapes.
double & initial_centre_of_mass(const unsigned &i)
Access function for the initial centre of mass.
void delete_external_hijacked_data()
Delete the storage for the external data formed from hijacked data.
~ImmersedRigidBodyElement()
Destuctor: Cleanup if required.
double Initial_Phi
Original rotation angle.
double *& st_pt()
Access function for the pointer to the fluid Strouhal number.
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.
double *& re_invfr_pt()
Access function for pointer to the fluid inverse Froude number (dimensionless gravitational loading)
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th (only) Data item that the object's shape depends on.
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.
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.
Vector< double > *& g_pt()
Access function to the direction of gravity.
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.
ExternalForceFctPt & external_force_fct_pt()
Access to function pointer to function that specifies external force.
const double & re_invfr()
Access to the fluid inverse Froude number.
Data *& centre_displacement_data_pt()
Pointer to Data for centre of gravity displacement. Values: 0: x-displ; 1: y-displ; 2: rotation angle...
GeomObject * Geom_object_pt
Underlying geometric object.
void pin_rotation_angle()
Pin the rotation 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...
Vector< double > centre_of_gravity()
Get current centre of gravity.
double Moment_of_inertia
Polar moment of inertia of body.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Get the contribution to the residuals.
ExternalTorqueFctPt & external_torque_fct_pt()
Access to function pointer to function that specifies external torque.
void unset_geometric_rotation()
Set the rotation of the object to be ignored (only really useful if you have a circular shape)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Get residuals including contribution to jacobian.
const double & initial_centre_of_mass(const unsigned &i) const
Access function for the initial centre of mass (const version)
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.
double *& re_pt()
Access function for the pointer to the fluid Reynolds number.
Vector< double > * G_pt
The direction of gravity.
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...
Mesh * Drag_mesh_pt
Mesh containing face elements that allow the computation of the drag on the body.
double & centre_x_displacement()
x-displacement of centre of mass
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 & re() const
Access function for the fluid Reynolds number.
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...
static double Default_Physical_Ratio_Value
Static default value for physical ratios.
double *& density_ratio_pt()
Access function for the pointer to the density ratio.
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.
Vector< double > vertex_coordinate(const unsigned &i) const
Coordinate vector of i-th vertex (const version)
unsigned nvertex() const
Number of vertices.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
unsigned npolyline() const
Number of constituent polylines.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...