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