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