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