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-2022 Matthias Heil and Andrew Hazel
7//LIC//
8//LIC// This library is free software; you can redistribute it and/or
9//LIC// modify it under the terms of the GNU Lesser General Public
10//LIC// License as published by the Free Software Foundation; either
11//LIC// version 2.1 of the License, or (at your option) any later version.
12//LIC//
13//LIC// This library is distributed in the hope that it will be useful,
14//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16//LIC// Lesser General Public License for more details.
17//LIC//
18//LIC// You should have received a copy of the GNU Lesser General Public
19//LIC// License along with this library; if not, write to the Free Software
20//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21//LIC// 02110-1301 USA.
22//LIC//
23//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24//LIC//
25//LIC//====================================================================
26
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
40using namespace std;
41using namespace oomph;
42
43
44namespace 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//====================================================================
407template<class ELEMENT>
408class DropInChannelProblem : public Problem
409{
410
411public:
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 }
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
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
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
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
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
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
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
523
524 //Unhijack the data in the internal element
525 Hijacked_element_pt->unhijack_all_data();
526
527 //Delete the volume constraint elements
529
530 // Internal pressure is gone -- null it out
532
533 // Kill the mesh too
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
554private:
555
556
557 /// 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
583
584 /// Delete volume constraint elements
586 {
587 if (Volume_constraint_mesh_pt==0) return;
588
589 // Backup
591 {
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
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 {
652 };
653
654
655}; // end_of_problem_class
656
657
658//==start_constructor=====================================================
659/// Constructor
660//========================================================================
661template<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//=======================================================================
907template<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//=======================================================================
955template<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//========================================================================
1024template<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//========================================================================
1190template<class ELEMENT>
1191void 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(),
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(),
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//========================================================================
1267template<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//============================================================
1303int 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.