jeffery_orbit.cc
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-2023 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 // Driver to demonstrate the interaction of a fluid flow with a rigid body.
27 // The driver solves a classic problem of the motion of an rigid ellipse in
28 // a shear flow. The problem imposes a smooth, but rapid, start-up from
29 // no flow to a shear imposed on the boundaries of a rectangular domain
30 // and the ellipse soon settles into approximate Jeffery orbits for small
31 // enough Reynolds number.
32 
33 //Generic routines
34 #include "generic.h"
35 
36 // The equations
37 #include "navier_stokes.h"
38 #include "solid.h"
39 #include "constitutive.h"
40 #include "rigid_body.h"
41 
42 // The mesh
43 #include "meshes/triangle_mesh.h"
44 
45 //Thin wrapper to "custom" TaylorHoodElements that overload output functions
47 
48 #include <algorithm>
49 
50 using namespace std;
51 using namespace oomph;
52 
53 //==start_of_namespace==============================
54 /// Namespace for Problem Parameters
55 //==================================================
57  {
58  /// Reynolds number
59  double Re=1.0;
60 
61  /// Strouhal number
62  double St = 1.0;
63 
64  /// Density ratio (Solid density / Fluid density)
65  double Density_ratio = 1.0;
66 
67  /// Initial axis of the elliptical solid in x-direction
68  double A = 0.25;
69 
70  /// Initial axis of the elliptical solid in y-direction
71  /// (N.B. 2B = 1 is the reference length scale)
72  double B = 0.5;
73 
74  /// Pseudo-solid (mesh) Poisson ratio
75  double Nu=0.3;
76 
77  /// Pseudo-solid (mesh) "density"
78  /// Set to zero because we don't want inertia in the node update!
79  double Lambda_sq=0.0;
80 
81  /// Constitutive law used to determine the mesh deformation
82  ConstitutiveLaw *Constitutive_law_pt=
83  new GeneralisedHookean(&Problem_Parameter::Nu);
84 
85  } // end_of_namespace
86 
87 
88 //=======================================================================
89 /// Exact solution for the rotation of an ellipse in unbounded shear flow
90 /// as computed by Jeffery (1922)
91 //=======================================================================
93 {
94  /// Null function
95  double null(const double &t) {return 0.0;}
96 
97  /// Angular position as a function of time t
98  double angle(const double &t)
99  {
100  const double a = Problem_Parameter::A;
101  const double b = Problem_Parameter::B;
102 
103  return atan((b/a)*tan((a*b*t)/(b*b + a*a)));
104  }
105 
106  /// Angular velocity as function of time t
107  double velocity(const double &t)
108  {
109  const double a = Problem_Parameter::A;
110  const double b = Problem_Parameter::B;
111 
112  //Get the angle
113  double chi = angle(t);
114 
115  //Now return the velocity
116  return (a*a*sin(chi)*sin(chi) + b*b*cos(chi)*cos(chi))/(a*a + b*b);
117  }
118 
119  /// Angular acceleration as a function of time t (should always be zero)
120  double acceleration(const double &t)
121  {
122  const double a = Problem_Parameter::A;
123  const double b = Problem_Parameter::A;
124 
125  //Get the angle and velocity
126  double chi = angle(t);
127  double chi_dot = velocity(t);
128 
129  //Now return the acceleration
130  return 2.0*(a*a - b*b)*(sin(chi)*cos(chi))*chi_dot/(a*a + b*b);
131  }
132 }
133 
134 /// ////////////////////////////////////////////////////////
135 /// ////////////////////////////////////////////////////////
136 /// ////////////////////////////////////////////////////////
137 
138 
139 //===================start_of_general_ellipse=================================
140 /// A geometric object for an ellipse with initial centre of mass at
141 /// (centre_x, centre_y) with axis in the x direction given by 2a
142 /// and in the y-direction given by 2b. The boundary of the ellipse is
143 /// parametrised by its angle.
144 //============================================================================
145 class GeneralEllipse : public GeomObject
146 {
147 private:
148 
149  //Storage for the centre of mass and semi-major and semi-minor axes
150  double Centre_x, Centre_y, A, B;
151 
152 public:
153 
154  /// Simple Constructor that transfers appropriate geometric
155  /// parameters into internal data
156  GeneralEllipse(const double &centre_x, const double &centre_y,
157  const double &a, const double &b)
158  : GeomObject(1,2), Centre_x(centre_x), Centre_y(centre_y), A(a), B(b)
159  {}
160 
161  /// Empty Destructor
163 
164  /// Return the position of the ellipse boundary as a function of
165  /// the angle xi[0]
166  void position(const Vector<double> &xi, Vector<double> &r) const
167  {
168  r[0] = Centre_x + A*cos(xi[0]);
169  r[1] = Centre_y + B*sin(xi[0]);
170  }
171 
172  //Return the position which is always fixed
173  void position(const unsigned &t,
174  const Vector<double> &xi, Vector<double> &r) const
175  {
176  return position(xi,r);
177  }
178 
179 };
180 //end_of_general_ellipse
181 
182 
183 /// ////////////////////////////////////////////////////////
184 /// ////////////////////////////////////////////////////////
185 /// ////////////////////////////////////////////////////////
186 
187 
188 //==start_of_problem_class============================================
189 /// Unstructured Navier-Stokes ALE Problem for a rigid ellipse
190 /// immersed within a viscous fluid
191 //====================================================================
192 template<class ELEMENT>
194 {
195 
196 public:
197 
198  /// Constructor
200 
201  /// Destructor
203 
204  /// Reset the boundary conditions when timestepping
206  {
207  this->set_boundary_velocity();
208  }
209 
210  /// Wipe the meshes of Lagrange multiplier and drag elements
211  void actions_before_adapt();
212 
213  /// Rebuild the meshes of Lagrange multiplier and drag elements
214  void actions_after_adapt();
215 
216  /// Re-apply the no slip condition (imposed indirectly via dependent
217  /// velocities)
219  {
220  // Update mesh -- this applies the auxiliary node update function
221  Fluid_mesh_pt->node_update();
222  }
223 
224  /// Set boundary condition, assign auxiliary node update fct.
225  /// Complete the build of all elements, attach power elements that allow
226  /// computation of drag vector
227  void complete_problem_setup();
228 
229  /// Set the boundary velocity
230  void set_boundary_velocity();
231 
232  /// Function that solves a simplified problem to ensure that
233  /// the positions of the boundary nodes are initially consistent with
234  /// the lagrange multiplier formulation
235  void solve_for_consistent_nodal_positions();
236 
237  /// Doc the solution
238  void doc_solution(const bool& project=false);
239 
240  /// Output the exact solution
241  void output_exact_solution(std::ofstream &output_file);
242 
243 private:
244 
245  /// Create elements that enforce prescribed boundary motion
246  /// for the pseudo-solid fluid mesh by Lagrange multipliers
247  void create_lagrange_multiplier_elements();
248 
249  /// Delete elements that impose the prescribed boundary displacement
250  /// and wipe the associated mesh
251  void delete_lagrange_multiplier_elements();
252 
253  /// Create elements that calculate the drag and torque on
254  /// the boundaries
255  void create_drag_elements();
256 
257  /// Delete elements that calculate the drag and torque on the
258  /// boundaries
259  void delete_drag_elements();
260 
261  /// Pin the degrees of freedom associated with the solid bodies
262  void pin_rigid_body();
263 
264  /// Unpin the degrees of freedom associated with the solid bodies
265  void unpin_rigid_body();
266 
267  /// Pointers to mesh of Lagrange multiplier elements
269 
270  /// Pointer to Fluid_mesh
271  RefineableSolidTriangleMesh<ELEMENT>* Fluid_mesh_pt;
272 
273  /// Triangle mesh polygon for outer boundary
274  TriangleMeshPolygon* Outer_boundary_polygon_pt;
275 
276  /// Mesh of drag elements
277  Vector<Mesh*> Drag_mesh_pt;
278 
279  /// Mesh of the generalised elements for the rigid bodies
281 
282  /// Storage for the geom object
283  Vector<GeomObject*> Rigid_body_pt;
284 
285  /// Internal DocInfo object
286  DocInfo Doc_info;
287 
288  /// File to document the norm of the solution (for validation purposes)
289  ofstream Norm_file;
290 
291  /// File to document the motion of the centre of gravity
292  ofstream Cog_file;
293 
294  /// File to document the exact motion of the centre of gravity
295  ofstream Cog_exact_file;
296 
297 }; // end_of_problem_class
298 
299 
300 //==start_constructor=====================================================
301 /// Constructor: Open output files, construct time steppers, build
302 /// fluid mesh, immersed rigid body and combine to form the problem
303 //========================================================================
304 template<class ELEMENT>
307 {
308  // Output directory
309  this->Doc_info.set_directory("RESLT");
310 
311  // Open norm file
312  this->Norm_file.open("RESLT/norm.dat");
313 
314  // Open file to trace the centre of gravity
315  this->Cog_file.open("RESLT/cog_trace.dat");
316 
317  // Open file to document the exact motion of the centre of gravity
318  this->Cog_exact_file.open("RESLT/cog_exact_trace.dat");
319 
320  // Allocate the timestepper -- this constructs the Problem's
321  // time object with a sufficient amount of storage to store the
322  // previous timsteps.
323  this->add_time_stepper_pt(new BDF<2>);
324 
325  // Allocate a timestepper for the rigid body
326  this->add_time_stepper_pt(new Newmark<2>);
327 
328  // Define the boundaries: Polyline with 4 different
329  // boundaries for the outer boundary and 1 internal elliptical hole
330 
331  // Build the boundary segments for outer boundary, consisting of
332  //--------------------------------------------------------------
333  // four separate polyline segments
334  //---------------------------------
335  Vector<TriangleMeshCurveSection*> boundary_segment_pt(4);
336 
337  //Set the length of the channel
338  double half_length = 5.0;
339  double half_height = 2.5;
340 
341  // Initialize boundary segment
342  Vector<Vector<double> > bound_seg(2);
343  for(unsigned i=0;i<2;i++) {bound_seg[i].resize(2);}
344 
345  // First boundary segment
346  bound_seg[0][0]=-half_length;
347  bound_seg[0][1]=-half_height;
348  bound_seg[1][0]=-half_length;
349  bound_seg[1][1]=half_height;
350 
351  // Specify 1st boundary id
352  unsigned bound_id = 0;
353 
354  // Build the 1st boundary segment
355  boundary_segment_pt[0] = new TriangleMeshPolyLine(bound_seg,bound_id);
356 
357  // Second boundary segment
358  bound_seg[0][0]=-half_length;
359  bound_seg[0][1]=half_height;
360  bound_seg[1][0]=half_length;
361  bound_seg[1][1]=half_height;
362 
363  // Specify 2nd boundary id
364  bound_id = 1;
365 
366  // Build the 2nd boundary segment
367  boundary_segment_pt[1] = new TriangleMeshPolyLine(bound_seg,bound_id);
368 
369  // Third boundary segment
370  bound_seg[0][0]=half_length;
371  bound_seg[0][1]=half_height;
372  bound_seg[1][0]=half_length;
373  bound_seg[1][1]=-half_height;
374 
375  // Specify 3rd boundary id
376  bound_id = 2;
377 
378  // Build the 3rd boundary segment
379  boundary_segment_pt[2] = new TriangleMeshPolyLine(bound_seg,bound_id);
380 
381  // Fourth boundary segment
382  bound_seg[0][0]=half_length;
383  bound_seg[0][1]=-half_height;
384  bound_seg[1][0]=-half_length;
385  bound_seg[1][1]=-half_height;
386 
387  // Specify 4th boundary id
388  bound_id = 3;
389 
390  // Build the 4th boundary segment
391  boundary_segment_pt[3] = new TriangleMeshPolyLine(bound_seg,bound_id);
392 
393  // Create the triangle mesh polygon for outer boundary using boundary segment
394  Outer_boundary_polygon_pt = new TriangleMeshPolygon(boundary_segment_pt);
395 
396  // Now build the moving rigid body
397  //-------------------------------------
398 
399  // We have one rigid body
400  Rigid_body_pt.resize(1);
401  Vector<TriangleMeshClosedCurve*> hole_pt(1);
402 
403  // Build Rigid Body
404  //-----------------
405  double x_center = 0.0;
406  double y_center = 0.0;
407  double A = Problem_Parameter::A;
408  double B = Problem_Parameter::B;
409  GeomObject* temp_hole_pt = new GeneralEllipse(x_center,y_center,A,B);
410  Rigid_body_pt[0] = new ImmersedRigidBodyElement(temp_hole_pt,
411  this->time_stepper_pt(1));
412 
413  // Build the two parts of the curvilinear boundary from the rigid body
414  Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt(2);
415 
416  //First section (boundary 4)
417  double zeta_start=0.0;
418  double zeta_end=MathematicalConstants::Pi;
419  unsigned nsegment=8;
420  unsigned boundary_id=4;
421  curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
422  Rigid_body_pt[0],zeta_start,zeta_end,nsegment,boundary_id);
423 
424  //Second section (boundary 5)
425  zeta_start=MathematicalConstants::Pi;
426  zeta_end=2.0*MathematicalConstants::Pi;
427  nsegment=8;
428  boundary_id=5;
429  curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
430  Rigid_body_pt[0],zeta_start,zeta_end,
431  nsegment,boundary_id);
432 
433  // Combine to form a hole in the fluid mesh
434  Vector<double> hole_coords(2);
435  hole_coords[0]=0.0;
436  hole_coords[1]=0.0;
437  Vector<TriangleMeshClosedCurve*> curvilinear_hole_pt(1);
438  hole_pt[0]=
439  new TriangleMeshClosedCurve(
440  curvilinear_boundary_pt,hole_coords);
441 
442  // Now build the mesh, based on the boundaries specified by
443  //---------------------------------------------------------
444  // polygons just created
445  //----------------------
446 
447  TriangleMeshClosedCurve* closed_curve_pt=Outer_boundary_polygon_pt;
448 
449  double uniform_element_area=1.0;
450 
451  // Use the TriangleMeshParameters object for gathering all
452  // the necessary arguments for the TriangleMesh object
453  TriangleMeshParameters triangle_mesh_parameters(
454  closed_curve_pt);
455 
456  // Define the holes on the domain
457  triangle_mesh_parameters.internal_closed_curve_pt() =
458  hole_pt;
459 
460  // Define the maximum element area
461  triangle_mesh_parameters.element_area() =
462  uniform_element_area;
463 
464  // Create the mesh
465  Fluid_mesh_pt =
466  new RefineableSolidTriangleMesh<ELEMENT>(
467  triangle_mesh_parameters, this->time_stepper_pt());
468 
469  // Set error estimator for bulk mesh
470  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
471  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
472 
473  // Set targets for spatial adaptivity
474  Fluid_mesh_pt->max_permitted_error()=0.005;
475  Fluid_mesh_pt->min_permitted_error()=0.001;
476  Fluid_mesh_pt->max_element_size()=1.0;
477  Fluid_mesh_pt->min_element_size()=0.001;
478 
479  // Use coarser mesh during validation
480  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
481  {
482  Fluid_mesh_pt->min_element_size()=0.01;
483  }
484 
485  // Set boundary condition, assign auxiliary node update fct,
486  // complete the build of all elements, attach power elements that allow
487  // computation of drag vector
488  complete_problem_setup();
489 
490  //Set the parameters of the rigid body elements
491  ImmersedRigidBodyElement* rigid_element1_pt =
492  dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0]);
493  rigid_element1_pt->initial_centre_of_mass(0) = x_center;
494  rigid_element1_pt->initial_centre_of_mass(1) = y_center;
495  rigid_element1_pt->mass_shape() = MathematicalConstants::Pi*A*B;
496  rigid_element1_pt->moment_of_inertia_shape() =
497  0.25*MathematicalConstants::Pi*A*B*(A*A + B*B);
498  rigid_element1_pt->re_pt() = &Problem_Parameter::Re;
499  rigid_element1_pt->st_pt() = &Problem_Parameter::St;
500  rigid_element1_pt->density_ratio_pt() = &Problem_Parameter::Density_ratio;
501 
502  //Pin the position of the centre of mass
503  rigid_element1_pt->pin_centre_of_mass_coordinate(0);
504  rigid_element1_pt->pin_centre_of_mass_coordinate(1);
505 
506  // Create the mesh for the rigid bodies
507  Rigid_body_mesh_pt = new Mesh;
508  Rigid_body_mesh_pt->add_element_pt(rigid_element1_pt);
509 
510  // Create the drag mesh for the rigid bodies
511  Drag_mesh_pt.resize(1);
512  for(unsigned m=0;m<1;m++) {Drag_mesh_pt[m] = new Mesh;}
513  this->create_drag_elements();
514 
515  //Add the drag mesh to the appropriate rigid bodies
516  rigid_element1_pt->set_drag_mesh(Drag_mesh_pt[0]);
517 
518 
519 
520  // Create Lagrange multiplier mesh for boundary motion
521  //----------------------------------------------------
522  // Construct the mesh of elements that enforce prescribed boundary motion
523  // of pseudo-solid fluid mesh by Lagrange multipliers
524  Lagrange_multiplier_mesh_pt=new SolidMesh;
525  create_lagrange_multiplier_elements();
526 
527 
528  // Combine meshes
529  //---------------
530 
531  // Add Fluid_mesh_pt sub meshes
532  this->add_sub_mesh(Fluid_mesh_pt);
533 
534  // Add Lagrange_multiplier sub meshes
535  this->add_sub_mesh(this->Lagrange_multiplier_mesh_pt);
536 
537  this->add_sub_mesh(this->Rigid_body_mesh_pt);
538 
539  // Build global mesh
540  this->build_global_mesh();
541 
542  // Setup equation numbering scheme
543  cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
544 
545 } // end_of_constructor
546 
547 
548 //========================================================================
549 /// Destructor that cleans up memory and closes files
550 //========================================================================
551 template<class ELEMENT>
554 {
555  // Delete Fluid timestepper
556  delete this->time_stepper_pt(0);
557 
558  //Delete solid timestepper
559  delete this->time_stepper_pt(1);
560 
561  // Kill data associated with outer boundary
562  unsigned n=Outer_boundary_polygon_pt->npolyline();
563  for (unsigned j=0;j<n;j++)
564  {
565  delete Outer_boundary_polygon_pt->polyline_pt(j);
566  }
567  delete Outer_boundary_polygon_pt;
568 
569  //Flush the drag mesh from the rigid body
570  dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
571  flush_drag_mesh();
572 
573  // Flush Lagrange multiplier mesh
574  delete_lagrange_multiplier_elements();
575  delete Lagrange_multiplier_mesh_pt;
576 
577  // Delete error estimator
578  delete Fluid_mesh_pt->spatial_error_estimator_pt();
579 
580  // Delete fluid mesh
581  delete Fluid_mesh_pt;
582 
583  // Kill const eqn
585 
586  // Close norm and trace files
587  this->Cog_exact_file.close();
588  this->Cog_file.close();
589  this->Norm_file.close();
590 }
591 
592 
593 //==========start_actions_before_adapt==================================
594 /// Actions before adapt: Wipe the mesh of Lagrange multiplier elements
595 //=======================================================================
596 template<class ELEMENT>
598 {
599  //Flush the drag mesh from the rigid body
600  dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
601  flush_drag_mesh();
602 
603  //Kill the drag element
604  delete_drag_elements();
605 
606  // Kill the elements and wipe surface mesh
607  delete_lagrange_multiplier_elements();
608 
609  // Rebuild the Problem's global mesh from its various sub-meshes
610  this->rebuild_global_mesh();
611 } // end of actions_before_adapt
612 
613 
614 //============start_actions_after_adapt=====================================
615 /// Actions after adapt: Rebuild the mesh of Lagrange multiplier elements
616 //==========================================================================
617 template<class ELEMENT>
619 {
620  //Reset the lagrangian coordinates for the solid mechanics
621  //an updated lagrangian approach
622  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
623 
624  // Create the elements that impose the displacement constraint
625  create_lagrange_multiplier_elements();
626 
627  // Create the drag elements anew
628  create_drag_elements();
629 
630  // Add the drag elements to the rigid body
631  dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
632  set_drag_mesh(this->Drag_mesh_pt[0]);
633 
634  // Rebuild the Problem's global mesh from its various sub-meshes
635  this->rebuild_global_mesh();
636 
637  // Setup the problem again -- remember that fluid mesh has been
638  // completely rebuilt and its element's don't have any
639  // pointers to Re etc. yet
640  complete_problem_setup();
641 
642  // Output solution after adaptation/projection
643  bool doc_projection=true;
644  doc_solution(doc_projection);
645 }// end of actions_after_adapt
646 
647 
648 //============start_complete_problem_setup=================================
649 /// Set boundary condition, assign auxiliary node update fct.
650 /// Complete the build of all elements, attach power elements that allow
651 /// computation of drag vector
652 //=========================================================================
653 template<class ELEMENT>
655 {
656  // Set the boundary conditions for fluid problem: All nodes are
657  // free by default -- just pin the ones that have Dirichlet conditions
658  // here.
659  unsigned nbound=Fluid_mesh_pt->nboundary();
660  for(unsigned ibound=0;ibound<nbound;ibound++)
661  {
662  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
663  for (unsigned inod=0;inod<num_nod;inod++)
664  {
665  // Cache pointer to node
666  Node* const nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
667 
668  //Pin x-velocity unless on inlet (0) and outlet (2) boundaries
669  //of the external rectangular box
670  if((ibound!=0) && (ibound!=2)) {nod_pt->pin(0);}
671  //Pin the y-velocity on all boundaries
672  nod_pt->pin(1);
673 
674  // Pin pseudo-solid positions apart from on the
675  // ellipse boundary that is allowed to move
676  // Cache cast pointer to solid node
677  SolidNode* const solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
678 
679  //Pin the solid positions on all external boundaries
680  if(ibound < 4)
681  {
682  solid_node_pt->pin_position(0);
683  solid_node_pt->pin_position(1);
684  }
685  // Unpin the position of all the nodes on hole boundaries:
686  // since they will be moved using Lagrange Multiplier
687  else
688  {
689  solid_node_pt->unpin_position(0);
690  solid_node_pt->unpin_position(1);
691 
692  // Assign auxiliary node update fct, which determines the
693  // velocity on the moving boundary using the position history
694  // values
695  // A more accurate version may be obtained by using velocity
696  // based on the actual position of the geometric object,
697  // but this introduces additional dependencies between the
698  // Data of the rigid body and the fluid elements.
699  nod_pt->set_auxiliary_node_update_fct_pt(
700  FSI_functions::apply_no_slip_on_moving_wall);
701  }
702  } //End of loop over boundary nodes
703  } // End loop over boundaries
704 
705  // Complete the build of all elements so they are fully functional
706  unsigned n_element = Fluid_mesh_pt->nelement();
707  for(unsigned e=0;e<n_element;e++)
708  {
709  // Upcast from GeneralisedElement to the present element
710  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
711 
712  // Set the Reynolds number
713  el_pt->re_pt() = &Problem_Parameter::Re;
714 
715  // Set the Womersley number (same as Re since St=1)
716  el_pt->re_st_pt() = &Problem_Parameter::Re;
717 
718  // Set the constitutive law for pseudo-elastic mesh deformation
719  el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
720 
721  // Set the "density" for pseudo-elastic mesh deformation
722  el_pt->lambda_sq_pt()=&Problem_Parameter::Lambda_sq;
723  }
724 
725  // Re-apply Dirichlet boundary conditions for current and history values
726  // (projection ignores boundary conditions!)
727  this->set_boundary_velocity();
728 
729 } //end_of_complete_problem_setup
730 
731 
732 //=================start_set_boundary_velocity===========================
733 /// Set the boundary velocity for current and history values
734 //=======================================================================
735 template<class ELEMENT>
738 {
739  //Loop over top and bottom walls, inlet and outlet
740  for(unsigned ibound=0;ibound<4;ibound++)
741  {
742  //If we are on the inlet or outlet only zero the
743  //y velocity, leave x alone
744  if((ibound==0) || (ibound==2))
745  {
746  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
747  for (unsigned inod=0;inod<num_nod;inod++)
748  {
749  // Get node
750  Node* const nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
751 
752  // Get number of previous (history) values
753  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
754 
755  //Zero all current and previous y-velocities
756  for(unsigned t=0;t<=n_prev;t++)
757  {
758  nod_pt->set_value(t,1,0.0);
759  }
760  }
761  }
762  //Otherwise on the top and bottom walls set a smooth x-velocity
763  //and zero y-velocity
764  else
765  {
766  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
767  for (unsigned inod=0;inod<num_nod;inod++)
768  {
769  // Get node
770  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
771 
772  // Get number of previous (history) values
773  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
774 
775  //Now set the boundary velocity
776  double y = nod_pt->x(1);
777  //Get the previous time
778  for(unsigned t=0;t<=n_prev;t++)
779  {
780  //Get the time
781  double time_ = this->time_pt()->time(t);
782 
783  //Get the velocity ramp
784  //Initially zero (nothing at all is going on)
785  double ramp = 0.0;
786  double delta = 5.0;
787 
788  double e1 = exp(-delta);
789  double a1 = 1.0 - (1.0 + delta + 0.5*delta*delta)*e1;
790  double b1 = (3.0 + 3.0*delta + delta*delta)*e1 - 3.0;
791  double c1 = 3.0 - (3.0 + 2.0*delta + 0.5*delta*delta)*e1;
792  //Smooth start
793  if((time_ >= 0.0) & (time_ <= 1.0))
794  {
795  ramp = a1*time_*time_*time_
796  + b1*time_*time_
797  + c1*time_;
798  }
799  //Coupled to exponential levelling
800  else if (time_ > 1.0)
801  {
802  ramp = 1.0 - exp(-delta*time_);
803  }
804 
805  nod_pt->set_value(t,0,-y*ramp);
806  }
807 
808  // Zero all current and previous veloc values
809  // for the v-velocity
810  for (unsigned t=0;t<=n_prev;t++)
811  {
812  nod_pt->set_value(t,1,0.0);
813  }
814  }
815  }
816  }
817 }
818 
819 
820 //====================start_of_pin_rigid_body=====================
821 /// Pin the degrees of freedom associated with the solid bodies
822 //================================================================
823 template<class ELEMENT>
825 {
826  unsigned n_rigid_body = Rigid_body_pt.size();
827  for(unsigned i=0;i<n_rigid_body;++i)
828  {
829  unsigned n_geom_data = Rigid_body_pt[i]->ngeom_data();
830  for(unsigned j=0;j<n_geom_data;j++)
831  {
832  Rigid_body_pt[i]->geom_data_pt(j)->pin_all();
833  }
834  }
835 }
836 
837 
838 //=================start_unpin_rigid_body==========================
839 /// Unpin the degrees of freedom associated with the solid bodies
840 //=================================================================
841 template<class ELEMENT>
843 {
844  unsigned n_rigid_body = Rigid_body_pt.size();
845  for(unsigned i=0;i<n_rigid_body;++i)
846  {
847  unsigned n_geom_data = Rigid_body_pt[i]->ngeom_data();
848  for(unsigned j=0;j<n_geom_data;j++)
849  {
850  Rigid_body_pt[i]->geom_data_pt(j)->unpin_all();
851  }
852  }
853 }
854 
855 
856 //==========start_solve_for_consistent_nodal_positions================
857 /// Assemble and solve a simplified problem that ensures that the
858 /// positions of the boundary nodes are consistent with the weak
859 /// imposition of the displacement boundary conditions on the surface
860 /// of the ellipse.
861 //===================================================================
862 template<class ELEMENT>
865 {
866  //First pin all degrees of freedom in the rigid body
867  this->pin_rigid_body();
868 
869  //Must reassign equation numbrs
870  this->assign_eqn_numbers();
871 
872  //Do a steady solve to map the nodes to the boundary of the ellipse
873  this->steady_newton_solve();
874 
875  //Now unpin the rigid body...
876  this->unpin_rigid_body();
877 
878  //...and then repin the position of the centre of mass
879  ImmersedRigidBodyElement* rigid_element1_pt =
880  dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0]);
881  rigid_element1_pt->pin_centre_of_mass_coordinate(0);
882  rigid_element1_pt->pin_centre_of_mass_coordinate(1);
883 
884  //and then reassign equation numbers
885  this->assign_eqn_numbers();
886 
887 } //end_solve_for_consistent_nodal_positions
888 
889 
890 
891 
892 //============start_of_create_lagrange_multiplier_elements===============
893 /// Create elements that impose the prescribed boundary displacement
894 /// for the pseudo-solid fluid mesh
895 //=======================================================================
896 template<class ELEMENT>
899 {
900  // The idea is to apply Lagrange multipliers to the boundaries in
901  // the mesh that have associated geometric objects
902 
903  //Find the number of boundaries
904  unsigned n_boundary = Fluid_mesh_pt->nboundary();
905 
906  // Loop over the boundaries
907  for(unsigned b=0;b<n_boundary;b++)
908  {
909  //Get the geometric object associated with the boundary
910  GeomObject* boundary_geom_obj_pt =
911  Fluid_mesh_pt->boundary_geom_object_pt(b);
912 
913  //Only bother to do anything if there is a geometric object
914  if(boundary_geom_obj_pt!=0)
915  {
916  // How many bulk fluid elements are adjacent to boundary b?
917  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
918 
919  // Loop over the bulk fluid elements adjacent to boundary b?
920  for(unsigned e=0;e<n_element;e++)
921  {
922  // Get pointer to the bulk fluid element that is
923  // adjacent to boundary b
924  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
925  Fluid_mesh_pt->boundary_element_pt(b,e));
926 
927  //Find the index of the face of element e along boundary b
928  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
929 
930  // Create new element. Note that we use different Lagrange
931  // multiplier fields for each distinct boundary (here indicated
932  // by b.
933  ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>* el_pt =
934  new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
935  bulk_elem_pt,face_index,b);
936 
937  // Add it to the mesh
938  Lagrange_multiplier_mesh_pt->add_element_pt(el_pt);
939 
940  // Set the GeomObject that defines the boundary shape and set
941  // which bulk boundary we are attached to (needed to extract
942  // the boundary coordinate from the bulk nodes)
943  el_pt->set_boundary_shape_geom_object_pt(
944  boundary_geom_obj_pt,b);
945 
946  // Loop over the nodes to pin Lagrange multiplier
947  unsigned nnod=el_pt->nnode();
948  for(unsigned j=0;j<nnod;j++)
949  {
950  Node* nod_pt = el_pt->node_pt(j);
951 
952  // How many nodal values were used by the "bulk" element
953  // that originally created this node?
954  unsigned n_bulk_value=el_pt->nbulk_value(j);
955 
956  // Pin two of the four Lagrange multipliers at vertices
957  // This is not totally robust, but will work in this application
958  unsigned nval=nod_pt->nvalue();
959  if (nval==7)
960  {
961  for (unsigned i=0;i<2;i++)
962  {
963  // Pin lagrangian multipliers
964  nod_pt->pin(n_bulk_value+2+i);
965  }
966  }
967  }
968  } // end loop over the element
969  } //End of case if there is a geometric object
970  } //End of loop over boundaries
971 }
972 // end of create_lagrange_multiplier_elements
973 
974 
975 //===============start_delete_lagrange_multiplier_elements==================
976 /// Delete elements that impose the prescribed boundary displacement
977 /// and wipe the associated mesh
978 //==========================================================================
979 template<class ELEMENT>
982 {
983  // How many surface elements are in the surface mesh
984  unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
985 
986  // Loop over the surface elements
987  for(unsigned e=0;e<n_element;e++)
988  {
989  // Kill surface element
990  delete Lagrange_multiplier_mesh_pt->element_pt(e);
991  }
992 
993  // Wipe the mesh
994  Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
995 
996 } // end of delete_lagrange_multiplier_elements
997 
998 
999 
1000 
1001 //============start_of_create_drag_elements===============
1002 /// Create elements that calculate the drag and torque on
1003 /// the obstacles in the fluid mesh
1004 //=======================================================================
1005 template<class ELEMENT>
1007 {
1008  //The idea is only to attach drag elements to those
1009  //boundaries associated with each particular rigid body
1010  //The loop is slightly inefficient, but should work in general
1011 
1012  // Get the number of rigid bodies
1013  unsigned n_rigid = Rigid_body_pt.size();
1014 
1015  // Get the number of boundaries in the mesh
1016  unsigned n_boundary = Fluid_mesh_pt->nboundary();
1017 
1018  //Loop over the rigid bodies
1019  for(unsigned r=0;r<n_rigid;r++)
1020  {
1021  //Get the rigid boundary geometric object
1022  ImmersedRigidBodyElement* rigid_el_pt =
1023  dynamic_cast<ImmersedRigidBodyElement*>(this->Rigid_body_pt[r]);
1024 
1025  // Loop over all boundaries
1026  for(unsigned b=0;b<n_boundary;b++)
1027  {
1028  //Does the boundary correspond to the current rigid body
1029  if(dynamic_cast<ImmersedRigidBodyElement*>
1030  (Fluid_mesh_pt->boundary_geom_object_pt(b)) == rigid_el_pt)
1031  {
1032  // How many bulk fluid elements are adjacent to boundary b?
1033  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
1034 
1035  // Loop over the bulk fluid elements adjacent to boundary b?
1036  for(unsigned e=0;e<n_element;e++)
1037  {
1038  // Get pointer to the bulk fluid element that is
1039  // adjacent to boundary b
1040  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1041  Fluid_mesh_pt->boundary_element_pt(b,e));
1042 
1043  //Find the index of the face of element e along boundary b
1044  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
1045 
1046  // Create new element. Note that we use different Lagrange
1047  // multiplier fields for each distinct boundary (here indicated
1048  // by b.
1049  NavierStokesSurfaceDragTorqueElement<ELEMENT>* el_pt =
1050  new NavierStokesSurfaceDragTorqueElement<ELEMENT>(
1051  bulk_elem_pt,face_index);
1052 
1053  // Add it to the mesh
1054  Drag_mesh_pt[r]->add_element_pt(el_pt);
1055 
1056  //Set the original centre of rotation
1057  //from the rigid body in the drag element
1058  for(unsigned i=0;i<2;i++)
1059  {el_pt->centre_of_rotation(i) =
1060  rigid_el_pt->initial_centre_of_mass(i);}
1061 
1062  //Set the data to the translation and rotation data
1063  //as well because these data will enter the Jacobian terms
1064  //in the DragTorqueElement
1065  el_pt->set_translation_and_rotation(rigid_el_pt->geom_data_pt(0));
1066  } // end loop over the element
1067  }
1068  } //End of loop over boundaries
1069  } // end loop over the rigid bodies
1070 }
1071 // end of create_drag_elements
1072 
1073 
1074 //=======================================================================
1075 /// Delete elements that calculate the drag and torque on the
1076 /// boundaries
1077 //=======================================================================
1078 template<class ELEMENT>
1080 {
1081  unsigned n_bodies = Drag_mesh_pt.size();
1082  for(unsigned n=0;n<n_bodies;n++)
1083  {
1084  // How many surface elements are in the surface mesh
1085  unsigned n_element = Drag_mesh_pt[n]->nelement();
1086 
1087  // Loop over the surface elements
1088  for(unsigned e=0;e<n_element;e++)
1089  {
1090  // Kill surface element
1091  delete Drag_mesh_pt[n]->element_pt(e);
1092  }
1093 
1094  // Wipe the mesh
1095  Drag_mesh_pt[n]->flush_element_and_node_storage();
1096  }
1097 } // end of delete_drag_elements
1098 
1099 
1100 
1101 
1102 //==start_of_doc_solution=================================================
1103 /// Doc the solution
1104 //========================================================================
1105 template<class ELEMENT>
1107  const bool& project)
1108 {
1109 
1110  oomph_info << "Docing step: " << this->Doc_info.number()
1111  << std::endl;
1112 
1113  ofstream some_file;
1114  char filename[100];
1115 
1116  // Number of plot points
1117  unsigned npts;
1118  npts=5;
1119 
1120 
1121  // Output solution and projection files
1122  if(!project)
1123  {
1124  sprintf(filename,"%s/soln%i.dat",
1125  this->Doc_info.directory().c_str(),
1126  this->Doc_info.number());
1127  }
1128  else
1129  {
1130  sprintf(filename,"%s/proj%i.dat",
1131  this->Doc_info.directory().c_str(),
1132  this->Doc_info.number()-1);
1133  }
1134 
1135  // Assemble square of L2 norm
1136  double square_of_l2_norm=0.0;
1137  unsigned nel=Fluid_mesh_pt->nelement();
1138  for (unsigned e=0;e<nel;e++)
1139  {
1140  square_of_l2_norm+=
1141  dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(e))->
1142  square_of_l2_norm();
1143  }
1144 
1145  std::cout << "Checking " << sqrt(square_of_l2_norm) << "\n";
1146  this->Norm_file << sqrt(square_of_l2_norm) << "\n";
1147 
1148  some_file.open(filename);
1149  some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1150  ->variable_identifier();
1151  this->Fluid_mesh_pt->output(some_file,npts);
1152  some_file.close();
1153 
1154  // No trace file writing after projection
1155  if(project) return;
1156 
1157  //Output the boundary only if not projecting
1158  sprintf(filename,"%s/int%i.dat",
1159  this->Doc_info.directory().c_str(),
1160  this->Doc_info.number());
1161  some_file.open(filename);
1162  this->Lagrange_multiplier_mesh_pt->output(some_file);
1163  some_file.close();
1164 
1165  // Increment the doc_info number
1166  this->Doc_info.number()++;
1167 
1168  //Output the motion of the centre of gravity
1169  dynamic_cast<ImmersedRigidBodyElement*>(this->Rigid_body_pt[0])->
1170  output_centre_of_gravity(this->Cog_file);
1171 
1172  //Output the exact solution
1173  this->output_exact_solution(this->Cog_exact_file);
1174 
1175 }
1176 
1177 
1178 //=====================================================================
1179 /// Output the exact solution
1180 //=====================================================================
1181 template<class ELEMENT>
1183 output_exact_solution(std::ofstream &output_file)
1184  {
1185  //Get the current time
1186  double time = this->time();
1187  //Output the exact solution
1188  output_file << time << " " << Jeffery_Solution::angle(time)
1189  << " " << Jeffery_Solution::velocity(time)
1190  << " " << Jeffery_Solution::acceleration(time) << std::endl;
1191  }
1192 
1193 
1194 
1195 //==========start_of_main======================================
1196 /// Driver code for immersed ellipse problem
1197 //============================================================
1198 int main(int argc, char **argv)
1199 {
1200  // Store command line arguments
1201  CommandLineArgs::setup(argc,argv);
1202 
1203  // Define possible command line arguments and parse the ones that
1204  // were actually specified
1205 
1206  // Validation?
1207  CommandLineArgs::specify_command_line_flag("--validation");
1208 
1209  // Parse command line
1210  CommandLineArgs::parse_and_assign();
1211 
1212  // Doc what has actually been specified on the command line
1213  CommandLineArgs::doc_specified_flags();
1214 
1215  // Create problem in initial configuration
1217  ProjectableTaylorHoodElement<MyTaylorHoodElement> > problem;
1218 
1219  //Initially ensure that the nodal positions are consistent with
1220  //their weak imposition
1222 
1223  // Initialise timestepper
1224  double dt=0.05;
1225  problem.initialise_dt(dt);
1226 
1227  // Perform impulsive start
1228  problem.assign_initial_values_impulsive();
1229 
1230  // Output initial conditions
1231  problem.doc_solution();
1232 
1233  // Solve problem a few times on given mesh
1234  unsigned nstep=3;
1235  for (unsigned i=0;i<nstep;i++)
1236  {
1237  // Solve the problem
1238  problem.unsteady_newton_solve(dt);
1239  problem.doc_solution();
1240  }
1241 
1242  // Now do a couple of adaptations
1243  unsigned ncycle=200;
1244  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1245  {
1246  ncycle=1;
1247  oomph_info << "Only doing one cycle during validation\n";
1248  }
1249 
1250  for (unsigned j=0;j<ncycle;j++)
1251  {
1252  // Adapt the problem
1253  problem.adapt();
1254 
1255  //Solve problem a few times
1256  for (unsigned i=0;i<nstep;i++)
1257  {
1258  // Solve the problem
1259  problem.unsteady_newton_solve(dt);
1260  problem.doc_solution();
1261  }
1262  }
1263 
1264 } //end of main
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
~GeneralEllipse()
Empty Destructor.
GeneralEllipse(const double &centre_x, const double &centre_y, const double &a, const double &b)
Simple Constructor that transfers appropriate geometric parameters into internal data.
void position(const Vector< double > &xi, Vector< double > &r) const
Return the position of the ellipse boundary as a function of the angle xi[0].
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
void pin_rigid_body()
Pin the degrees of freedom associated with the solid bodies.
void doc_solution(const bool &project=false)
Doc the solution.
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
void create_drag_elements()
Create elements that calculate the drag and torque on the boundaries.
DocInfo Doc_info
Internal DocInfo object.
void output_exact_solution(std::ofstream &output_file)
Output the exact solution.
void solve_for_consistent_nodal_positions()
Function that solves a simplified problem to ensure that the positions of the boundary nodes are init...
Mesh * Rigid_body_mesh_pt
Mesh of the generalised elements for the rigid bodies.
void actions_before_implicit_timestep()
Reset the boundary conditions when timestepping.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion for the pseudo-solid fluid mesh by Lagrange m...
ofstream Cog_exact_file
File to document the exact motion of the centre of gravity.
void complete_problem_setup()
Set boundary condition, assign auxiliary node update fct. Complete the build of all elements,...
void unpin_rigid_body()
Unpin the degrees of freedom associated with the solid bodies.
void delete_lagrange_multiplier_elements()
Delete elements that impose the prescribed boundary displacement and wipe the associated mesh.
UnstructuredImmersedEllipseProblem()
Constructor.
ofstream Norm_file
File to document the norm of the solution (for validation purposes)
void delete_drag_elements()
Delete elements that calculate the drag and torque on the boundaries.
Vector< Mesh * > Drag_mesh_pt
Mesh of drag elements.
void actions_before_adapt()
Wipe the meshes of Lagrange multiplier and drag elements.
Vector< GeomObject * > Rigid_body_pt
Storage for the geom object.
void actions_after_adapt()
Rebuild the meshes of Lagrange multiplier and drag elements.
void set_boundary_velocity()
Set the boundary velocity.
ofstream Cog_file
File to document the motion of the centre of gravity.
SolidMesh * Lagrange_multiplier_mesh_pt
Pointers to mesh of Lagrange multiplier elements.
TriangleMeshPolygon * Outer_boundary_polygon_pt
Triangle mesh polygon for outer boundary.
void actions_before_newton_convergence_check()
Re-apply the no slip condition (imposed indirectly via dependent velocities)
int main(int argc, char **argv)
Driver code for immersed ellipse problem.
Exact solution for the rotation of an ellipse in unbounded shear flow as computed by Jeffery (1922)
double velocity(const double &t)
Angular velocity as function of time t.
double angle(const double &t)
Angular position as a function of time t.
double acceleration(const double &t)
Angular acceleration as a function of time t (should always be zero)
Namespace for Problem Parameters.
double Lambda_sq
Pseudo-solid (mesh) "density" Set to zero because we don't want inertia in the node update!
double A
Initial axis of the elliptical solid in x-direction.
double Density_ratio
Density ratio (Solid density / Fluid density)
double St
Strouhal number.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
double Nu
Pseudo-solid (mesh) Poisson ratio.
double B
Initial axis of the elliptical solid in y-direction (N.B. 2B = 1 is the reference length scale)
double Re
Reynolds number.