adaptive_drop_in_channel.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 
27 //Generic routines
28 #include "generic.h"
29 
30 
31 // The equations
32 #include "navier_stokes.h"
33 #include "solid.h"
34 #include "constitutive.h"
35 #include "fluid_interface.h"
36 
37 // The mesh
38 #include "meshes/triangle_mesh.h"
39 
40 using namespace std;
41 using namespace oomph;
42 
43 
44 namespace oomph
45 {
46 
47 
48 //==============================================================
49 /// Overload CrouzeixRaviart element to modify output
50 //==============================================================
52  public virtual PseudoSolidNodeUpdateElement<TCrouzeixRaviartElement<2>,
53  TPVDBubbleEnrichedElement<2,3> >
54  {
55 
56  private:
57 
58  /// Storage for elemental error estimate -- used for post-processing
59  double Error;
60 
61  public:
62 
63  /// Constructor initialise error
64  MyCrouzeixRaviartElement() : PseudoSolidNodeUpdateElement<TCrouzeixRaviartElement<2>,
65  TPVDBubbleEnrichedElement<2,3> >()
66  {
67  Error=0.0;
68  }
69 
70  /// Set error value for post-processing
71  void set_error(const double& error){Error=error;}
72 
73  /// Return variable identifier
74  std::string variable_identifier()
75  {
76  std::string txt="VARIABLES=";
77  txt+="\"x\",";
78  txt+="\"y\",";
79  txt+="\"u\",";
80  txt+="\"v\",";
81  txt+="\"p\",";
82  txt+="\"du/dt\",";
83  txt+="\"dv/dt\",";
84  txt+="\"u_m\",";
85  txt+="\"v_m\",";
86  txt+="\"x_h1\",";
87  txt+="\"y_h1\",";
88  txt+="\"x_h2\",";
89  txt+="\"y_h2\",";
90  txt+="\"u_h1\",";
91  txt+="\"v_h1\",";
92  txt+="\"u_h2\",";
93  txt+="\"v_h2\",";
94  txt+="\"error\",";
95  txt+="\"size\",";
96  txt+="\n";
97  return txt;
98  }
99 
100 
101  /// Overload output function
102  void output(std::ostream &outfile,
103  const unsigned &nplot)
104  {
105 
106  // Assign dimension
107  unsigned el_dim=2;
108 
109  // Vector of local coordinates
110  Vector<double> s(el_dim);
111 
112  // Acceleration
113  Vector<double> dudt(el_dim);
114 
115  // Mesh elocity
116  Vector<double> mesh_veloc(el_dim,0.0);
117 
118  // Tecplot header info
119  outfile << tecplot_zone_string(nplot);
120 
121  // Find out how many nodes there are
122  unsigned n_node = nnode();
123 
124  //Set up memory for the shape functions
125  Shape psif(n_node);
126  DShape dpsifdx(n_node,el_dim);
127 
128  // Loop over plot points
129  unsigned num_plot_points=nplot_points(nplot);
130  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
131  {
132 
133  // Get local coordinates of plot point
134  get_s_plot(iplot,nplot,s);
135 
136  //Call the derivatives of the shape and test functions
137  dshape_eulerian(s,psif,dpsifdx);
138 
139  //Allocate storage
140  Vector<double> mesh_veloc(el_dim);
141  Vector<double> dudt(el_dim);
142  Vector<double> dudt_ALE(el_dim);
143  DenseMatrix<double> interpolated_dudx(el_dim,el_dim);
144 
145  //Initialise everything to zero
146  for(unsigned i=0;i<el_dim;i++)
147  {
148  mesh_veloc[i]=0.0;
149  dudt[i]=0.0;
150  dudt_ALE[i]=0.0;
151  for(unsigned j=0;j<el_dim;j++)
152  {
153  interpolated_dudx(i,j) = 0.0;
154  }
155  }
156 
157  //Calculate velocities and derivatives
158 
159  //Loop over directions
160  for(unsigned i=0;i<el_dim;i++)
161  {
162  //Get the index at which velocity i is stored
163  unsigned u_nodal_index = u_index_nst(i);
164  // Loop over nodes
165  for(unsigned l=0;l<n_node;l++)
166  {
167  dudt[i]+=du_dt_nst(l,u_nodal_index)*psif[l];
168  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
169 
170  //Loop over derivative directions for velocity gradients
171  for(unsigned j=0;j<el_dim;j++)
172  {
173  interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*
174  dpsifdx(l,j);
175  }
176  }
177  }
178 
179 
180  // Get dudt in ALE form (incl mesh veloc)
181  for(unsigned i=0;i<el_dim;i++)
182  {
183  dudt_ALE[i]=dudt[i];
184  for (unsigned k=0;k<el_dim;k++)
185  {
186  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
187  }
188  }
189 
190 
191  // Coordinates
192  for(unsigned i=0;i<el_dim;i++)
193  {
194  outfile << interpolated_x(s,i) << " ";
195  }
196 
197  // Velocities
198  for(unsigned i=0;i<el_dim;i++)
199  {
200  outfile << interpolated_u_nst(s,i) << " ";
201  }
202 
203  // Pressure
204  outfile << interpolated_p_nst(s) << " ";
205 
206  // Accelerations
207  for(unsigned i=0;i<el_dim;i++)
208  {
209  outfile << dudt_ALE[i] << " ";
210  }
211 
212  // Mesh velocity
213  for(unsigned i=0;i<el_dim;i++)
214  {
215  outfile << mesh_veloc[i] << " ";
216  }
217 
218  // History values of coordinates
219  unsigned n_prev=node_pt(0)->position_time_stepper_pt()->ntstorage();
220  for (unsigned t=1;t<n_prev;t++)
221  {
222  for(unsigned i=0;i<el_dim;i++)
223  {
224  outfile << interpolated_x(t,s,i) << " ";
225  }
226  }
227 
228  // History values of velocities
229  n_prev=node_pt(0)->time_stepper_pt()->ntstorage();
230  for (unsigned t=1;t<n_prev;t++)
231  {
232  for(unsigned i=0;i<el_dim;i++)
233  {
234  outfile << interpolated_u_nst(t,s,i) << " ";
235  }
236  }
237 
238  outfile << Error << " "
239  << size() << std::endl;
240  }
241 
242  // Write tecplot footer (e.g. FE connectivity lists)
243  write_tecplot_zone_footer(outfile,nplot);
244  }
245 
246 
247  /// Get square of L2 norm of velocity
249  {
250 
251  // Assign dimension
252  unsigned el_dim=2;
253  // Initalise
254  double sum=0.0;
255 
256  //Find out how many nodes there are
257  unsigned n_node = nnode();
258 
259  //Find the indices at which the local velocities are stored
260  unsigned u_nodal_index[el_dim];
261  for(unsigned i=0;i<el_dim;i++) {u_nodal_index[i] = u_index_nst(i);}
262 
263  //Set up memory for the velocity shape fcts
264  Shape psif(n_node);
265  DShape dpsidx(n_node,el_dim);
266 
267  //Number of integration points
268  unsigned n_intpt = integral_pt()->nweight();
269 
270  //Set the Vector to hold local coordinates
271  Vector<double> s(el_dim);
272 
273  //Loop over the integration points
274  for(unsigned ipt=0;ipt<n_intpt;ipt++)
275  {
276  //Assign values of s
277  for(unsigned i=0;i<el_dim;i++) s[i] = integral_pt()->knot(ipt,i);
278 
279  //Get the integral weight
280  double w = integral_pt()->weight(ipt);
281 
282  // Call the derivatives of the veloc shape functions
283  // (Derivs not needed but they are free)
284  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
285 
286  //Premultiply the weights and the Jacobian
287  double W = w*J;
288 
289  //Calculate velocities
290  Vector<double> interpolated_u(el_dim,0.0);
291 
292  // Loop over nodes
293  for(unsigned l=0;l<n_node;l++)
294  {
295  //Loop over directions
296  for(unsigned i=0;i<el_dim;i++)
297  {
298  //Get the nodal value
299  double u_value = raw_nodal_value(l,u_nodal_index[i]);
300  interpolated_u[i] += u_value*psif[l];
301  }
302  }
303 
304  //Assemble square of L2 norm
305  for(unsigned i=0;i<el_dim;i++)
306  {
307  sum+=interpolated_u[i]*interpolated_u[i]*W;
308  }
309  }
310 
311  return sum;
312 
313  }
314 
315  };
316 
317 
318 //=======================================================================
319 /// Face geometry for element is the same as that for the underlying
320 /// wrapped element
321 //=======================================================================
322  template<>
323  class FaceGeometry<MyCrouzeixRaviartElement>
324  : public virtual SolidTElement<1,3>
325  {
326  public:
327  FaceGeometry() : SolidTElement<1,3>() {}
328  };
329 
330 
331 //=======================================================================
332 /// Face geometry of Face geometry for element is the same
333 /// as that for the underlying
334 /// wrapped element
335 //=======================================================================
336  template<>
337  class FaceGeometry<FaceGeometry<MyCrouzeixRaviartElement> >
338  : public virtual SolidPointElement
339  {
340  public:
341  FaceGeometry() : SolidPointElement() {}
342  };
343 
344 
345 } //End of namespace extension
346 
347 
348 
349 /// ////////////////////////////////////////////////////////
350 /// ////////////////////////////////////////////////////////
351 /// ////////////////////////////////////////////////////////
352 
353 
354 //==start_of_namespace==============================
355 /// Namespace for Problem Parameter
356 //==================================================
357  namespace Problem_Parameter
358  {
359  /// Doc info
360  DocInfo Doc_info;
361 
362  /// Reynolds number
363  double Re=0.0;
364 
365  /// Capillary number
366  double Ca = 1.0;
367 
368  /// Pseudo-solid Poisson ratio
369  double Nu=0.3;
370 
371  /// Initial radius of drop
372  double Radius = 0.25;
373 
374  /// Viscosity ratio of the droplet to the surrounding fluid
375  double Viscosity_ratio = 10.0;
376 
377  /// Volume of the interface
378  /// (Negative because the constraint is computed from the bulk side)
379  double Volume = -MathematicalConstants::Pi*Radius*Radius;
380 
381  /// Scaling factor for inflow velocity (allows it to be switched off
382  /// to do hydrostatics)
383  double Inflow_veloc_magnitude = 0.0;
384 
385  /// Length of the channel
386  double Length = 3.0;
387 
388  /// Constitutive law used to determine the mesh deformation
389  ConstitutiveLaw *Constitutive_law_pt=0;
390 
391  /// Trace file
392  ofstream Trace_file;
393 
394  /// File to document the norm of the solution (for validation
395  /// purposes -- triangle doesn't give fully reproducible results so
396  /// mesh generation/adaptation may generate slightly different numbers
397  /// of elements on different machines!)
398  ofstream Norm_file;
399 
400  } // end_of_namespace
401 
402 
403 
404 //==start_of_problem_class============================================
405 /// Problem class to simulate viscous drop propagating along 2D channel
406 //====================================================================
407 template<class ELEMENT>
408 class DropInChannelProblem : public Problem
409 {
410 
411 public:
412 
413  /// Constructor
415 
416  /// Destructor
418  {
419  // Fluid timestepper
420  delete this->time_stepper_pt(0);
421 
422  // Kill data associated with outer boundary
423  unsigned n=Outer_boundary_polyline_pt->npolyline();
424  for (unsigned j=0;j<n;j++)
425  {
426  delete Outer_boundary_polyline_pt->polyline_pt(j);
427  }
428  delete Outer_boundary_polyline_pt;
429 
430  //Kill data associated with drops
431  unsigned n_drop = Drop_polygon_pt.size();
432  for(unsigned idrop=0;idrop<n_drop;idrop++)
433  {
434  unsigned n=Drop_polygon_pt[idrop]->npolyline();
435  for (unsigned j=0;j<n;j++)
436  {
437  delete Drop_polygon_pt[idrop]->polyline_pt(j);
438  }
439  delete Drop_polygon_pt[idrop];
440  }
441 
442  // Flush element of free surface elements
443  delete_free_surface_elements();
444  delete Free_surface_mesh_pt;
445  delete_volume_constraint_elements();
446  delete Volume_constraint_mesh_pt;
447 
448  // Delete error estimator
449  delete Fluid_mesh_pt->spatial_error_estimator_pt();
450 
451  // Delete fluid mesh
452  delete Fluid_mesh_pt;
453 
454  // Delete the global pressure drop data
455  delete Drop_pressure_data_pt;
456 
457  // Kill const eqn
459 
460  }
461 
462 
463  /// Actions before adapt: Wipe the mesh of free surface elements
465  {
466  // Kill the elements and wipe surface mesh
467  delete_free_surface_elements();
468  delete_volume_constraint_elements();
469 
470  // Rebuild the Problem's global mesh from its various sub-meshes
471  this->rebuild_global_mesh();
472 
473  }// end of actions_before_adapt
474 
475 
476  /// Actions after adapt: Rebuild the mesh of free surface elements
478  {
479  // Create the elements that impose the displacement constraint
480  create_free_surface_elements();
481  create_volume_constraint_elements();
482 
483 
484  // Rebuild the Problem's global mesh from its various sub-meshes
485  this->rebuild_global_mesh();
486 
487  // Setup the problem again -- remember that fluid mesh has been
488  // completely rebuilt and its element's don't have any
489  // pointers to Re etc. yet
490  complete_problem_setup();
491 
492  }// end of actions_after_adapt
493 
494 
495  /// Update the after solve (empty)
497 
498  /// Update the problem specs before solve
500  {
501  //Reset the Lagrangian coordinates of the nodes to be the current
502  //Eulerian coordinates -- this makes the current configuration
503  //stress free
504  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
505  }
506 
507  /// Set boundary conditions and complete the build of all elements
508  void complete_problem_setup();
509 
510  /// Doc the solution
511  void doc_solution(const std::string& comment="");
512 
513  /// Compute the error estimates and assign to elements for plotting
514  void compute_error_estimate(double& max_err,
515  double& min_err);
516 
517 
518  /// Change the boundary conditions to remove the volume constraint
520  {
521  // Ignore all volume constraint stuff from here onwards
522  Use_volume_constraint=false;
523 
524  //Unhijack the data in the internal element
525  Hijacked_element_pt->unhijack_all_data();
526 
527  //Delete the volume constraint elements
528  delete_volume_constraint_elements();
529 
530  // Internal pressure is gone -- null it out
531  Drop_pressure_data_pt=0;
532 
533  // Kill the mesh too
534  delete Volume_constraint_mesh_pt;
535  Volume_constraint_mesh_pt=0;
536 
537  //Remove the sub meshes
538  this->flush_sub_meshes();
539 
540  // Add Fluid_mesh_pt sub meshes
541  this->add_sub_mesh(Fluid_mesh_pt);
542 
543  // Add Free_surface sub meshes
544  this->add_sub_mesh(this->Free_surface_mesh_pt);
545 
546  //Rebuild the global mesh
547  this->rebuild_global_mesh();
548 
549  //Renumber the equations
550  std::cout << "Removed volume constraint to obtain "
551  << this->assign_eqn_numbers() << " new equation numbers\n";
552  }
553 
554 private:
555 
556 
557  /// Create free surface elements
558  void create_free_surface_elements();
559 
560  /// Delete free surface elements
562  {
563  // How many surface elements are in the surface mesh
564  unsigned n_element = Free_surface_mesh_pt->nelement();
565 
566  // Loop over the surface elements
567  for(unsigned e=0;e<n_element;e++)
568  {
569  // Kill surface element
570  delete Free_surface_mesh_pt->element_pt(e);
571  }
572 
573 
574  // Wipe the mesh
575  Free_surface_mesh_pt->flush_element_and_node_storage();
576 
577 
578  } // end of delete_free_surface_elements
579 
580 
581 /// Create elements that impose volume constraint on the drop
582  void create_volume_constraint_elements();
583 
584  /// Delete volume constraint elements
586  {
587  if (Volume_constraint_mesh_pt==0) return;
588 
589  // Backup
590  if (Vol_constraint_el_pt!=0)
591  {
592  Initial_value_for_drop_pressure=Vol_constraint_el_pt->p_traded();
593  }
594 
595  // How many surface elements are in the surface mesh
596  unsigned n_element = Volume_constraint_mesh_pt->nelement();
597 
598  // Loop over all
599  unsigned first_el_to_be_killed=0;
600  for(unsigned e=first_el_to_be_killed;e<n_element;e++)
601  {
602  delete Volume_constraint_mesh_pt->element_pt(e);
603  }
604 
605  // We've just killed the volume constraint element
606  Vol_constraint_el_pt=0;
607 
608  // Wipe the mesh
609  Volume_constraint_mesh_pt->flush_element_and_node_storage();
610 
611  } // end of delete_volume_constraint_elements
612 
613  /// Pointers to mesh of free surface elements
615 
616  /// Pointer to mesh containing elements that impose volume constraint
618 
619  /// Pointer to Fluid_mesh
620  RefineableSolidTriangleMesh<ELEMENT>* Fluid_mesh_pt;
621 
622  /// Vector storing pointer to the drop polygons
623  Vector<TriangleMeshPolygon*> Drop_polygon_pt;
624 
625  /// Triangle mesh polygon for outer boundary
626  TriangleMeshPolygon* Outer_boundary_polyline_pt;
627 
628  /// Pointer to a global drop pressure datum
630 
631  /// Pointer to element that imposes volume constraint for drop
632  VolumeConstraintElement* Vol_constraint_el_pt;
633 
634  /// Backed up drop pressure between adaptations
636 
637  /// Pointer to hijacked element
639 
640  /// Bool to indicate if volume constraint is applied (only for steady run)
642 
643  /// Enumeration of channel boundaries
644  enum
645  {
646  Inflow_boundary_id=0,
647  Upper_wall_boundary_id=1,
648  Outflow_boundary_id=2,
649  Bottom_wall_boundary_id=3,
650  First_drop_boundary_id=4,
651  Second_drop_boundary_id=5
652  };
653 
654 
655 }; // end_of_problem_class
656 
657 
658 //==start_constructor=====================================================
659 /// Constructor
660 //========================================================================
661 template<class ELEMENT>
663 {
664  // Output directory
665  Problem_Parameter::Doc_info.set_directory("RESLT");
666 
667  // Allocate the timestepper -- this constructs the Problem's
668  // time object with a sufficient amount of storage to store the
669  // previous timsteps.
670  this->add_time_stepper_pt(new BDF<2>);
671 
672 
673 
674  // Build the boundary segments for outer boundary, consisting of
675  //--------------------------------------------------------------
676  // four separate polylines
677  //------------------------
678  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
679 
680  // Each polyline only has two vertices -- provide storage for their
681  // coordinates
682  Vector<Vector<double> > vertex_coord(2);
683  for(unsigned i=0;i<2;i++)
684  {
685  vertex_coord[i].resize(2);
686  }
687 
688  // First polyline: Inflow
689  vertex_coord[0][0]=0.0;
690  vertex_coord[0][1]=0.0;
691  vertex_coord[1][0]=0.0;
692  vertex_coord[1][1]=1.0;
693 
694  // Build the 1st boundary polyline
695  boundary_polyline_pt[0] = new TriangleMeshPolyLine(vertex_coord,
696  Inflow_boundary_id);
697 
698  // Second boundary polyline: Upper wall
699  vertex_coord[0][0]=0.0;
700  vertex_coord[0][1]=1.0;
701  vertex_coord[1][0]=Problem_Parameter::Length;
702  vertex_coord[1][1]=1.0;
703 
704  // Build the 2nd boundary polyline
705  boundary_polyline_pt[1] = new TriangleMeshPolyLine(vertex_coord,
706  Upper_wall_boundary_id);
707 
708  // Third boundary polyline: Outflow
709  vertex_coord[0][0]=Problem_Parameter::Length;
710  vertex_coord[0][1]=1.0;
711  vertex_coord[1][0]=Problem_Parameter::Length;
712  vertex_coord[1][1]=0.0;
713 
714  // Build the 3rd boundary polyline
715  boundary_polyline_pt[2] = new TriangleMeshPolyLine(vertex_coord,
716  Outflow_boundary_id);
717 
718  // Fourth boundary polyline: Bottom wall
719  vertex_coord[0][0]=Problem_Parameter::Length;
720  vertex_coord[0][1]=0.0;
721  vertex_coord[1][0]=0.0;
722  vertex_coord[1][1]=0.0;
723 
724  // Build the 4th boundary polyline
725  boundary_polyline_pt[3] = new TriangleMeshPolyLine(vertex_coord,
726  Bottom_wall_boundary_id);
727 
728  // Create the triangle mesh polygon for outer boundary
729  Outer_boundary_polyline_pt = new TriangleMeshPolygon(boundary_polyline_pt);
730 
731 
732  // Now define initial shape of drop(s) with polygon
733  //---------------------------------------------------
734 
735  // We have one drop
736  Drop_polygon_pt.resize(1);
737 
738  // Place it smack in the middle of the channel
739  double x_center = 0.5*Problem_Parameter::Length;
740  double y_center = 0.5;
741  Ellipse * drop_pt = new Ellipse(Problem_Parameter::Radius,
743 
744  // Intrinsic coordinate along GeomObject defining the drop
745  Vector<double> zeta(1);
746 
747  // Position vector to GeomObject defining the drop
748  Vector<double> coord(2);
749 
750  // Number of points defining drop
751  unsigned npoints = 16;
752  double unit_zeta = MathematicalConstants::Pi/double(npoints-1);
753 
754  // This drop is bounded by two distinct boundaries, each
755  // represented by its own polyline
756  Vector<TriangleMeshCurveSection*> drop_polyline_pt(2);
757 
758  // Vertex coordinates
759  Vector<Vector<double> > drop_vertex(npoints);
760 
761  // Create points on boundary
762  for(unsigned ipoint=0; ipoint<npoints;ipoint++)
763  {
764  // Resize the vector
765  drop_vertex[ipoint].resize(2);
766 
767  // Get the coordinates
768  zeta[0]=unit_zeta*double(ipoint);
769  drop_pt->position(zeta,coord);
770 
771  // Shift
772  drop_vertex[ipoint][0]=coord[0]+x_center;
773  drop_vertex[ipoint][1]=coord[1]+y_center;
774  }
775 
776  // Build the 1st drop polyline
777  drop_polyline_pt[0] = new TriangleMeshPolyLine(drop_vertex,
778  First_drop_boundary_id);
779 
780  // Second boundary of drop
781  for(unsigned ipoint=0; ipoint<npoints;ipoint++)
782  {
783  // Resize the vector
784  drop_vertex[ipoint].resize(2);
785 
786  // Get the coordinates
787  zeta[0]=(unit_zeta*double(ipoint))+MathematicalConstants::Pi;
788  drop_pt->position(zeta,coord);
789 
790  // Shift
791  drop_vertex[ipoint][0]=coord[0]+x_center;
792  drop_vertex[ipoint][1]=coord[1]+y_center;
793  }
794 
795  // Build the 2nd drop polyline
796  drop_polyline_pt[1] = new TriangleMeshPolyLine(drop_vertex,
797  Second_drop_boundary_id);
798 
799  // Create closed polygon from two polylines
800  Drop_polygon_pt[0] = new TriangleMeshPolygon(
801  drop_polyline_pt);
802 
803  // Now build the mesh, based on the boundaries specified by
804  //---------------------------------------------------------
805  // polygons just created
806  //----------------------
807 
808  // Convert to "closed curve" objects
809  TriangleMeshClosedCurve* outer_closed_curve_pt=Outer_boundary_polyline_pt;
810  unsigned nb=Drop_polygon_pt.size();
811  Vector<TriangleMeshClosedCurve*> drop_closed_curve_pt(nb);
812  for (unsigned i=0;i<nb;i++)
813  {
814  drop_closed_curve_pt[i]=Drop_polygon_pt[i];
815  }
816 
817  // Target area for initial mesh
818  double uniform_element_area=0.2;
819 
820  // Define coordinates of a point inside the drop
821  Vector<double> drop_center(2);
822  drop_center[0]=x_center;
823  drop_center[1]=y_center;
824 
825  // Use the TriangleMeshParameters object for gathering all
826  // the necessary arguments for the TriangleMesh object
827  TriangleMeshParameters triangle_mesh_parameters(
828  outer_closed_curve_pt);
829 
830  // Define the holes on the boundary
831  triangle_mesh_parameters.internal_closed_curve_pt() =
832  drop_closed_curve_pt;
833 
834  // Define the maximum element area
835  triangle_mesh_parameters.element_area() =
836  uniform_element_area;
837 
838  // Define the region
839  triangle_mesh_parameters.add_region_coordinates(1, drop_center);
840 
841  // Create the mesh
842  Fluid_mesh_pt =
843  new RefineableSolidTriangleMesh<ELEMENT>(
844  triangle_mesh_parameters, this->time_stepper_pt());
845 
846  // Set error estimator for bulk mesh
847  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
848  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
849 
850  // Set targets for spatial adaptivity
851  Fluid_mesh_pt->max_permitted_error()=0.005;
852  Fluid_mesh_pt->min_permitted_error()=0.001;
853  Fluid_mesh_pt->max_element_size()=0.2;
854  Fluid_mesh_pt->min_element_size()=0.001;
855 
856  // Use coarser mesh during validation
857  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
858  {
859  Fluid_mesh_pt->min_element_size()=0.01;
860  }
861 
862  // Output boundary and mesh initial mesh for information
863  this->Fluid_mesh_pt->output_boundaries("boundaries.dat");
864  this->Fluid_mesh_pt->output("mesh.dat");
865 
866  // Provide initial value for drop pressure
867  Initial_value_for_drop_pressure=Problem_Parameter::Ca/
869 
870  // Set boundary condition and complete the build of all elements
871  complete_problem_setup();
872 
873  // Construct the mesh of elements that impose the volume constraint
874  Use_volume_constraint=true;
875  Volume_constraint_mesh_pt = new Mesh;
876  create_volume_constraint_elements();
877 
878  // Construct the mesh of free surface elements
879  Free_surface_mesh_pt=new Mesh;
880  create_free_surface_elements();
881 
882  // Combine meshes
883  //---------------
884 
885  // Add volume constraint sub mesh
886  this->add_sub_mesh(this->Volume_constraint_mesh_pt);
887 
888  // Add Fluid_mesh_pt sub meshes
889  this->add_sub_mesh(Fluid_mesh_pt);
890 
891  // Add Free_surface sub meshes
892  this->add_sub_mesh(this->Free_surface_mesh_pt);
893 
894  // Build global mesh
895  this->build_global_mesh();
896 
897  // Setup equation numbering scheme
898  cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
899 
900 } // end_of_constructor
901 
902 
903 //============start_of_create_free_surface_elements===============
904 /// Create elements that impose the kinematic and dynamic bcs
905 /// for the pseudo-solid fluid mesh
906 //=======================================================================
907 template<class ELEMENT>
909 {
910 
911  //Loop over the free surface boundaries
912  unsigned nb=Fluid_mesh_pt->nboundary();
913  for(unsigned b=First_drop_boundary_id;b<nb;b++)
914  {
915  // Note: region is important
916  // How many bulk fluid elements are adjacent to boundary b in region 0?
917  unsigned n_element = Fluid_mesh_pt->nboundary_element_in_region(b,0);
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 in region 0
924  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
925  Fluid_mesh_pt->boundary_element_in_region_pt(b,0,e));
926 
927  //Find the index of the face of element e along boundary b in region 0
928  int face_index = Fluid_mesh_pt->face_index_at_boundary_in_region(b,0,e);
929 
930  // Create new element
931  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
932  new ElasticLineFluidInterfaceElement<ELEMENT>(
933  bulk_elem_pt,face_index);
934 
935  // Add it to the mesh
936  Free_surface_mesh_pt->add_element_pt(el_pt);
937 
938  //Add the appropriate boundary number
939  el_pt->set_boundary_number_in_bulk_mesh(b);
940 
941  //Specify the capillary number
942  el_pt->ca_pt() = &Problem_Parameter::Ca;
943  }
944  }
945 }
946 // end of create_free_surface_elements
947 
948 
949 
950 
951 
952 //============start_of_create_volume_constraint_elements=================
953 /// Create elements that impose volume constraint on the drop
954 //=======================================================================
955 template<class ELEMENT>
957 {
958 
959  // Do we need it?
960  if (!Use_volume_constraint) return;
961 
962  // Store pointer to element whose pressure we're trading/hi-jacking:
963  // Element 0 in region 1
964  Hijacked_element_pt= dynamic_cast<ELEMENT*>(
965  Fluid_mesh_pt->region_element_pt(1,0));
966 
967  // Set the global pressure data by hijacking one of the pressure values
968  // from inside the droplet
969  unsigned index_of_traded_pressure=0;
970  Drop_pressure_data_pt=Hijacked_element_pt->
971  hijack_internal_value(0,index_of_traded_pressure);
972 
973  // Build volume constraint element -- pass traded pressure to it
974  Vol_constraint_el_pt=
975  new VolumeConstraintElement(&Problem_Parameter::Volume,
976  Drop_pressure_data_pt,
977  index_of_traded_pressure);
978 
979  //Provide a reasonable initial guess for drop pressure (hydrostatics):
980  Drop_pressure_data_pt->set_value(index_of_traded_pressure,
981  Initial_value_for_drop_pressure);
982 
983  // Add volume constraint element to the mesh
984  Volume_constraint_mesh_pt->add_element_pt(Vol_constraint_el_pt);
985 
986  //Loop over the free surface boundaries
987  unsigned nb=Fluid_mesh_pt->nboundary();
988  for(unsigned b=First_drop_boundary_id;b<nb;b++)
989  {
990  // Note region is important
991  // How many bulk fluid elements are adjacent to boundary b in region 0?
992  unsigned n_element = Fluid_mesh_pt->nboundary_element_in_region(b,0);
993 
994  // Loop over the bulk fluid elements adjacent to boundary b?
995  for(unsigned e=0;e<n_element;e++)
996  {
997  // Get pointer to the bulk fluid element that is
998  // adjacent to boundary b in region 0?
999  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1000  Fluid_mesh_pt->boundary_element_in_region_pt(b,0,e));
1001 
1002  //Find the index of the face of element e along boundary b in region 0
1003  int face_index = Fluid_mesh_pt->face_index_at_boundary_in_region(b,0,e);
1004 
1005  // Create new element
1006  ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
1007  new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
1008  bulk_elem_pt,face_index);
1009 
1010  //Set the "master" volume constraint element
1011  el_pt->set_volume_constraint_element(Vol_constraint_el_pt);
1012 
1013  // Add it to the mesh
1014  Volume_constraint_mesh_pt->add_element_pt(el_pt);
1015  }
1016  }
1017 }
1018 // end of create_volume_constraint_elements
1019 
1020 
1021 //==start_of_complete_problem_setup=======================================
1022 /// Set boundary conditions and complete the build of all elements
1023 //========================================================================
1024 template<class ELEMENT>
1026  {
1027  // Map to record if a given boundary is on a drop or not
1028  map<unsigned,bool> is_on_drop_bound;
1029 
1030  // Loop over the drops
1031  unsigned ndrop=Drop_polygon_pt.size();
1032  for(unsigned idrop=0;idrop<ndrop;idrop++)
1033  {
1034  // Get the vector all boundary IDs associated with the polylines that
1035  // make up the closed polygon
1036  Vector<unsigned> drop_bound_id=this->Drop_polygon_pt[idrop]->
1037  polygon_boundary_id();
1038 
1039  // Get the number of boundary
1040  unsigned nbound=drop_bound_id.size();
1041 
1042  // Fill in the map
1043  for(unsigned ibound=0;ibound<nbound;ibound++)
1044  {
1045  // This boundary...
1046  unsigned bound_id=drop_bound_id[ibound];
1047 
1048  // ...is on the drop
1049  is_on_drop_bound[bound_id]=true;
1050  }
1051  }
1052 
1053  // Re-set the boundary conditions for fluid problem: All nodes are
1054  // free by default -- just pin the ones that have Dirichlet conditions
1055  // here.
1056  unsigned nbound=Fluid_mesh_pt->nboundary();
1057  for(unsigned ibound=0;ibound<nbound;ibound++)
1058  {
1059  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
1060  for (unsigned inod=0;inod<num_nod;inod++)
1061  {
1062  // Get node
1063  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1064 
1065  //Pin both velocities on inflow (0) and side boundaries (1 and 3)
1066  if((ibound==0) || (ibound==1) || (ibound==3))
1067  {
1068  nod_pt->pin(0);
1069  nod_pt->pin(1);
1070  }
1071 
1072  //If it's the outflow pin only the vertical velocity
1073  if(ibound==2) {nod_pt->pin(1);}
1074 
1075  // Pin pseudo-solid positions apart from drop boundary which
1076  // we allow to move
1077  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
1078  if(is_on_drop_bound[ibound])
1079  {
1080  solid_node_pt->unpin_position(0);
1081  solid_node_pt->unpin_position(1);
1082  }
1083  else
1084  {
1085  solid_node_pt->pin_position(0);
1086  solid_node_pt->pin_position(1);
1087  }
1088  }
1089  } // end loop over boundaries
1090 
1091  // Complete the build of all elements so they are fully functional
1092  // Remember that adaptation for triangle meshes involves a complete
1093  // regneration of the mesh (rather than splitting as in tree-based
1094  // meshes where such parameters can be passed down from the father
1095  // element!)
1096  unsigned n_element = Fluid_mesh_pt->nelement();
1097  for(unsigned e=0;e<n_element;e++)
1098  {
1099  // Upcast from GeneralisedElement to the present element
1100  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
1101 
1102  // Set the Reynolds number
1103  el_pt->re_pt() = &Problem_Parameter::Re;
1104 
1105  // Set the Womersley number (same as Re since St=1)
1106  el_pt->re_st_pt() = &Problem_Parameter::Re;
1107 
1108  // Set the constitutive law for pseudo-elastic mesh deformation
1109  el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
1110  }
1111 
1112  //For the elements within the droplet (region 1),
1113  //set the viscosity ratio
1114  n_element = Fluid_mesh_pt->nregion_element(1);
1115  for(unsigned e=0;e<n_element;e++)
1116  {
1117  // Upcast from GeneralisedElement to the present element
1118  ELEMENT* el_pt =
1119  dynamic_cast<ELEMENT*>(Fluid_mesh_pt->region_element_pt(1,e));
1120 
1121  el_pt->viscosity_ratio_pt() = &Problem_Parameter::Viscosity_ratio;
1122  }
1123 
1124 
1125  // Re-apply boundary values on Dirichlet boundary conditions
1126  // (Boundary conditions are ignored when the solution is transferred
1127  // from the old to the new mesh by projection; this leads to a slight
1128  // change in the boundary values (which are, of course, never changed,
1129  // unlike the actual unknowns for which the projected values only
1130  // serve as an initial guess)
1131 
1132  // Set velocity and history values of velocity on walls
1133  nbound=this->Fluid_mesh_pt->nboundary();
1134  for(unsigned ibound=0;ibound<nbound;++ibound)
1135  {
1136  if ((ibound==Upper_wall_boundary_id)||
1137  (ibound==Bottom_wall_boundary_id)||
1138  (ibound==Outflow_boundary_id) ||
1139  (ibound==Inflow_boundary_id))
1140  {
1141  // Loop over nodes on this boundary
1142  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
1143  for (unsigned inod=0;inod<num_nod;inod++)
1144  {
1145  // Get node
1146  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1147 
1148  // Get number of previous (history) values
1149  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
1150 
1151  // Velocity is and was zero at all previous times
1152  for (unsigned t=0;t<=n_prev;t++)
1153  {
1154  if (ibound!=Inflow_boundary_id)
1155  {
1156  // Parallel outflow
1157  if (ibound!=Outflow_boundary_id)
1158  {
1159  nod_pt->set_value(t,0,0.0);
1160  }
1161  nod_pt->set_value(t,1,0.0);
1162  }
1163 
1164  // Nodes have always been there...
1165  nod_pt->x(t,0)=nod_pt->x(0,0);
1166  nod_pt->x(t,1)=nod_pt->x(0,1);
1167  }
1168  }
1169  }
1170  }
1171 
1172  // Re-assign prescribed inflow velocity at inlet
1173  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(Inflow_boundary_id);
1174  for (unsigned inod=0;inod<num_nod;inod++)
1175  {
1176  // Get node
1177  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(Inflow_boundary_id,
1178  inod);
1179  //Now set the boundary velocity
1180  double y = nod_pt->x(1);
1181  nod_pt->set_value(0,Problem_Parameter::Inflow_veloc_magnitude*y*(1-y));
1182  }
1183  }
1184 
1185 
1186 
1187 //==start_of_doc_solution=================================================
1188 /// Doc the solution
1189 //========================================================================
1190 template<class ELEMENT>
1191 void DropInChannelProblem<ELEMENT>::doc_solution(const std::string& comment)
1192 {
1193 
1194  oomph_info << "Docing step: " << Problem_Parameter::Doc_info.number()
1195  << std::endl;
1196 
1197  ofstream some_file;
1198  char filename[100];
1199  sprintf(filename,"%s/soln%i.dat",
1200  Problem_Parameter::Doc_info.directory().c_str(),
1201  Problem_Parameter::Doc_info.number());
1202 
1203  // Number of plot points
1204  unsigned npts;
1205  npts=5;
1206 
1207  // Compute errors and assign to each element for plotting
1208  double max_err;
1209  double min_err;
1210  compute_error_estimate(max_err,min_err);
1211 
1212  // Assemble square of L2 norm
1213  double square_of_l2_norm=0.0;
1214  unsigned nel=Fluid_mesh_pt->nelement();
1215  for (unsigned e=0;e<nel;e++)
1216  {
1217  square_of_l2_norm+=
1218  dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(e))->
1219  square_of_l2_norm();
1220  }
1221  Problem_Parameter::Norm_file << sqrt(square_of_l2_norm) << std::endl;
1222 
1223 
1224  some_file.open(filename);
1225  some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1226  ->variable_identifier();
1227  this->Fluid_mesh_pt->output(some_file,npts);
1228  some_file << "TEXT X = 25, Y = 78, CS=FRAME T = \"Global Step "
1229  << Problem_Parameter::Doc_info.number() << " "
1230  << comment << "\"\n";
1231  some_file.close();
1232 
1233  // Output boundaries
1234  sprintf(filename,"%s/boundaries%i.dat",
1235  Problem_Parameter::Doc_info.directory().c_str(),
1236  Problem_Parameter::Doc_info.number());
1237  some_file.open(filename);
1238  this->Fluid_mesh_pt->output_boundaries(some_file);
1239  some_file.close();
1240 
1241 
1242  // Get max/min area
1243  double max_area;
1244  double min_area;
1245  Fluid_mesh_pt->max_and_min_element_size(max_area, min_area);
1246 
1247  // Write trace file
1249  << this->time_pt()->time() << " "
1250  << Fluid_mesh_pt->nelement() << " "
1251  << max_area << " "
1252  << min_area << " "
1253  << max_err << " "
1254  << min_err << " "
1255  << sqrt(square_of_l2_norm) << " "
1256  << std::endl;
1257 
1258  // Increment the doc_info number
1259  Problem_Parameter::Doc_info.number()++;
1260 
1261 
1262 }
1263 
1264 //========================================================================
1265 /// Compute error estimates and assign to elements for plotting
1266 //========================================================================
1267 template<class ELEMENT>
1269  double& min_err)
1270 {
1271  // Get error estimator
1272  ErrorEstimator* err_est_pt=Fluid_mesh_pt->spatial_error_estimator_pt();
1273 
1274  // Get/output error estimates
1275  unsigned nel=Fluid_mesh_pt->nelement();
1276  Vector<double> elemental_error(nel);
1277 
1278  // We need a dynamic cast, get_element_errors from the Fluid_mesh_pt
1279  // Dynamic cast is used because get_element_errors require a Mesh* ans
1280  // not a SolidMesh*
1281  Mesh* fluid_mesh_pt=dynamic_cast<Mesh*>(Fluid_mesh_pt);
1282  err_est_pt->get_element_errors(fluid_mesh_pt,
1283  elemental_error);
1284 
1285  // Set errors for post-processing and find extrema
1286  max_err=0.0;
1287  min_err=DBL_MAX;
1288  for (unsigned e=0;e<nel;e++)
1289  {
1290  dynamic_cast<MyCrouzeixRaviartElement*>(Fluid_mesh_pt->element_pt(e))->
1291  set_error(elemental_error[e]);
1292 
1293  max_err=std::max(max_err,elemental_error[e]);
1294  min_err=std::min(min_err,elemental_error[e]);
1295  }
1296 
1297 }
1298 
1299 
1300 //============================================================
1301 /// Driver code for moving block problem
1302 //============================================================
1303 int main(int argc, char **argv)
1304 {
1305 
1306  // Store command line arguments
1307  CommandLineArgs::setup(argc,argv);
1308 
1309  // Define possible command line arguments and parse the ones that
1310  // were actually specified
1311 
1312  // Validation?
1313  CommandLineArgs::specify_command_line_flag("--validation");
1314 
1315  // Parse command line
1316  CommandLineArgs::parse_and_assign();
1317 
1318  // Doc what has actually been specified on the command line
1319  CommandLineArgs::doc_specified_flags();
1320 
1321  // Create generalised Hookean constitutive equations
1323  new GeneralisedHookean(&Problem_Parameter::Nu);
1324 
1325  // Open trace file
1326  Problem_Parameter::Trace_file.open("RESLT/trace.dat");
1327 
1328  // Open norm file
1329  Problem_Parameter::Norm_file.open("RESLT/norm.dat");
1330 
1331 
1332  // Create problem in initial configuration
1333  DropInChannelProblem<Hijacked<ProjectableCrouzeixRaviartElement<
1334  MyCrouzeixRaviartElement> > > problem;
1335 
1336  // Before starting the time-integration we want to "inflate" it to form
1337  // a proper circular drop. We do this by setting the inflow to zero
1338  // and doing a steady solve (with one adaptation)
1340 
1341  problem.steady_newton_solve(1);
1342 
1343  //Set the Capillary number to 10 and resolve
1344  Problem_Parameter::Ca = 10.0;
1345  problem.steady_newton_solve();
1346 
1347  // If all went well, this should show us a nice circular drop
1348  // in a stationary fluid
1349  problem.doc_solution();
1350 
1351  // Switch off volume constraint
1352  problem.remove_volume_constraint();
1353 
1354 
1355  // Initialise timestepper
1356  double dt=0.025;
1357  problem.initialise_dt(dt);
1358 
1359  // Perform impulsive start from current state
1360  problem.assign_initial_values_impulsive();
1361 
1362  // Output initial conditions
1363  problem.doc_solution();
1364 
1365  // Now switch on the inflow and re-assign the boundary conditions
1366  // (Call to complete_problem_setup() is a bit expensive given that we
1367  // we only want to set the inflow velocity but who cares -- it's just
1368  // a one off.
1370  problem.complete_problem_setup();
1371 
1372  // Solve problem on fixed mesh
1373  unsigned nstep=6;
1374  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1375  {
1376  nstep=2;
1377  oomph_info << "Remeshing after every second step during validation\n";
1378  }
1379  for (unsigned i=0;i<nstep;i++)
1380  {
1381  // Solve the problem
1382  problem.unsteady_newton_solve(dt);
1383  problem.doc_solution();
1384  } // done solution on fixed mesh
1385 
1386 
1387  // Now do a proper loop, doing nstep timesteps before adapting/remeshing
1388  // and repeating the lot ncycle times
1389  unsigned ncycle=1000;
1390  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1391  {
1392  ncycle=1;
1393  oomph_info << "Only doing one cycle during validation\n";
1394  }
1395 
1396 
1397  // Do the cycles
1398  for(unsigned j=0;j<ncycle;j++)
1399  {
1400  // Allow up to one level of refinement for next solve
1401  unsigned max_adapt=1;
1402 
1403  //Solve problem a few times
1404  for (unsigned i=0;i<nstep;i++)
1405  {
1406  // Solve the problem
1407  problem.unsteady_newton_solve(dt,max_adapt,false);
1408 
1409 
1410  // Build the label for doc and output solution
1411  std::stringstream label;
1412  label << "Adaptation " <<j << " Step "<< i;
1413  problem.doc_solution(label.str());
1414 
1415  // No more refinement for the next nstep steps
1416  max_adapt=0;
1417  }
1418  }
1419 
1420 
1421 } //End of main
int main(int argc, char **argv)
Driver code for moving block problem.
Problem class to simulate viscous drop propagating along 2D channel.
double Initial_value_for_drop_pressure
Backed up drop pressure between adaptations.
void actions_before_newton_solve()
Update the problem specs before solve.
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
void complete_problem_setup()
Set boundary conditions and complete the build of all elements.
void doc_solution(const std::string &comment="")
Doc the solution.
void delete_free_surface_elements()
Delete free surface elements.
void compute_error_estimate(double &max_err, double &min_err)
Compute the error estimates and assign to elements for plotting.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of free surface elements.
void create_volume_constraint_elements()
Create elements that impose volume constraint on the drop.
Mesh * Free_surface_mesh_pt
Pointers to mesh of free surface elements.
void create_free_surface_elements()
Create free surface elements.
bool Use_volume_constraint
Bool to indicate if volume constraint is applied (only for steady run)
void actions_after_newton_solve()
Update the after solve (empty)
void remove_volume_constraint()
Change the boundary conditions to remove the volume constraint.
Vector< TriangleMeshPolygon * > Drop_polygon_pt
Vector storing pointer to the drop polygons.
Mesh * Volume_constraint_mesh_pt
Pointer to mesh containing elements that impose volume constraint.
VolumeConstraintElement * Vol_constraint_el_pt
Pointer to element that imposes volume constraint for drop.
Data * Drop_pressure_data_pt
Pointer to a global drop pressure datum.
TriangleMeshPolygon * Outer_boundary_polyline_pt
Triangle mesh polygon for outer boundary.
void delete_volume_constraint_elements()
Delete volume constraint elements.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of free surface elements.
ELEMENT * Hijacked_element_pt
Pointer to hijacked element.
Overload CrouzeixRaviart element to modify output.
void set_error(const double &error)
Set error value for post-processing.
double square_of_l2_norm()
Get square of L2 norm of velocity.
MyCrouzeixRaviartElement()
Constructor initialise error.
double Error
Storage for elemental error estimate – used for post-processing.
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
std::string variable_identifier()
Return variable identifier.
////////////////////////////////////////////////////////////// //////////////////////////////////////...
DocInfo Doc_info
Doc info object.
ofstream Norm_file
File to document the norm of the solution (for validation purposes – triangle doesn't give fully repr...
ofstream Trace_file
Trace file.
double Inflow_veloc_magnitude
Scaling factor for inflow velocity (allows it to be switched off to do hydrostatics)
double Length
Length of the channel.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
double Radius
Initial radius of bubble.
double Volume
Volume of the bubble (negative because it's outside the fluid!)
double Nu
Pseudo-solid Poisson ratio.
double Re
Reynolds number.
double Ca
Capillary number.
double Viscosity_ratio
Viscosity ratio of the droplet to the surrounding fluid.