32 #include "navier_stokes.h"
34 #include "constitutive.h"
35 #include "fluid_interface.h"
38 #include "meshes/triangle_mesh.h"
41 using namespace oomph;
52 public virtual PseudoSolidNodeUpdateElement<TTaylorHoodElement<2>,
75 std::string txt=
"VARIABLES=";
106 const unsigned &nplot)
113 Vector<double> s(el_dim);
116 Vector<double> dudt(el_dim);
119 Vector<double> mesh_veloc(el_dim,0.0);
122 outfile << tecplot_zone_string(nplot);
125 unsigned n_node = nnode();
129 DShape dpsifdx(n_node,el_dim);
132 unsigned num_plot_points=nplot_points(nplot);
133 for (
unsigned iplot=0;iplot<num_plot_points;iplot++)
137 get_s_plot(iplot,nplot,s);
140 dshape_eulerian(s,psif,dpsifdx);
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);
149 for(
unsigned i=0;i<el_dim;i++)
154 for(
unsigned j=0;j<el_dim;j++)
156 interpolated_dudx(i,j) = 0.0;
163 for(
unsigned i=0;i<el_dim;i++)
166 unsigned u_nodal_index = u_index_nst(i);
168 for(
unsigned l=0;l<n_node;l++)
170 dudt[i]+=du_dt_nst(l,u_nodal_index)*psif[l];
171 mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
174 for(
unsigned j=0;j<el_dim;j++)
176 interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*
184 for(
unsigned i=0;i<el_dim;i++)
187 for (
unsigned k=0;k<el_dim;k++)
189 dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
195 for(
unsigned i=0;i<el_dim;i++)
197 outfile << interpolated_x(s,i) <<
" ";
201 for(
unsigned i=0;i<el_dim;i++)
203 outfile << interpolated_u_nst(s,i) <<
" ";
207 outfile << interpolated_p_nst(s) <<
" ";
210 for(
unsigned i=0;i<el_dim;i++)
212 outfile << dudt_ALE[i] <<
" ";
216 for(
unsigned i=0;i<el_dim;i++)
218 outfile << mesh_veloc[i] <<
" ";
222 unsigned n_prev=node_pt(0)->position_time_stepper_pt()->ntstorage();
223 for (
unsigned t=1;t<n_prev;t++)
225 for(
unsigned i=0;i<el_dim;i++)
227 outfile << interpolated_x(t,s,i) <<
" ";
232 n_prev=node_pt(0)->time_stepper_pt()->ntstorage();
233 for (
unsigned t=1;t<n_prev;t++)
235 for(
unsigned i=0;i<el_dim;i++)
237 outfile << interpolated_u_nst(t,s,i) <<
" ";
242 for(
unsigned i=0;i<el_dim;i++)
244 for(
unsigned j=0;j<el_dim;j++)
246 outfile << interpolated_dudx(i,j) <<
" ";
250 outfile << Error <<
" "
251 << size() << std::endl;
255 write_tecplot_zone_footer(outfile,nplot);
269 unsigned n_node = nnode();
272 unsigned u_nodal_index[el_dim];
273 for(
unsigned i=0;i<el_dim;i++) {u_nodal_index[i] = u_index_nst(i);}
277 DShape dpsidx(n_node,el_dim);
280 unsigned n_intpt = integral_pt()->nweight();
283 Vector<double> s(el_dim);
286 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
289 for(
unsigned i=0;i<el_dim;i++) s[i] = integral_pt()->knot(ipt,i);
292 double w = integral_pt()->weight(ipt);
296 double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
302 Vector<double> interpolated_u(el_dim,0.0);
305 for(
unsigned l=0;l<n_node;l++)
308 for(
unsigned i=0;i<el_dim;i++)
311 double u_value = raw_nodal_value(l,u_nodal_index[i]);
312 interpolated_u[i] += u_value*psif[l];
317 for(
unsigned i=0;i<el_dim;i++)
319 sum+=interpolated_u[i]*interpolated_u[i]*W;
336 :
public virtual SolidTElement<1,3>
350 :
public virtual SolidPointElement
423 template<
class ELEMENT>
436 delete this->time_stepper_pt(0);
439 unsigned n=Outer_boundary_polyline_pt->npolyline();
440 for (
unsigned j=0;j<n;j++)
442 delete Outer_boundary_polyline_pt->polyline_pt(j);
444 delete Outer_boundary_polyline_pt;
447 unsigned n_bubble = Bubble_polygon_pt.size();
448 for(
unsigned ibubble=0;ibubble<n_bubble;ibubble++)
450 unsigned n=Bubble_polygon_pt[ibubble]->npolyline();
451 for (
unsigned j=0;j<n;j++)
453 delete Bubble_polygon_pt[ibubble]->polyline_pt(j);
455 delete Bubble_polygon_pt[ibubble];
459 delete_free_surface_elements();
460 delete Free_surface_mesh_pt;
461 delete_volume_constraint_elements();
462 delete Volume_constraint_mesh_pt;
465 delete Fluid_mesh_pt->spatial_error_estimator_pt();
468 delete Fluid_mesh_pt;
471 delete Bubble_pressure_data_pt;
483 delete_free_surface_elements();
484 delete_volume_constraint_elements();
487 this->rebuild_global_mesh();
496 create_free_surface_elements();
497 create_volume_constraint_elements();
500 this->rebuild_global_mesh();
505 complete_problem_setup();
519 Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
523 void complete_problem_setup();
526 void doc_solution(
const std::string& comment=
"");
529 void compute_error_estimate(
double& max_err,
536 void create_free_surface_elements();
542 unsigned n_element = Free_surface_mesh_pt->nelement();
545 for(
unsigned e=0;e<n_element;e++)
548 delete Free_surface_mesh_pt->element_pt(e);
552 Free_surface_mesh_pt->flush_element_and_node_storage();
558 void create_volume_constraint_elements();
564 unsigned n_element = Volume_constraint_mesh_pt->nelement();
568 unsigned first_el_to_be_killed=1;
569 for(
unsigned e=first_el_to_be_killed;e<n_element;e++)
571 delete Volume_constraint_mesh_pt->element_pt(e);
575 Volume_constraint_mesh_pt->flush_element_and_node_storage();
603 Inflow_boundary_id=0,
604 Upper_wall_boundary_id=1,
605 Outflow_boundary_id=2,
606 Bottom_wall_boundary_id=3,
607 First_bubble_boundary_id=4,
608 Second_bubble_boundary_id=5
618 template<
class ELEMENT>
627 this->add_time_stepper_pt(
new BDF<2>);
638 Bubble_pressure_data_pt =
new Data(1);
639 unsigned index_of_traded_pressure=0;
640 this->add_global_data(Bubble_pressure_data_pt);
643 Vol_constraint_el_pt=
645 Bubble_pressure_data_pt,
646 index_of_traded_pressure);
649 Bubble_pressure_data_pt->set_value(index_of_traded_pressure,
655 Vol_constraint_el_pt=
659 unsigned index=Vol_constraint_el_pt->index_of_traded_pressure();
662 Bubble_pressure_data_pt=Vol_constraint_el_pt->p_traded_data_pt();
665 Vol_constraint_el_pt->p_traded_data_pt()->
675 Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
679 Vector<Vector<double> > vertex_coord(2);
680 for(
unsigned i=0;i<2;i++)
682 vertex_coord[i].resize(2);
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;
692 boundary_polyline_pt[0] =
new TriangleMeshPolyLine(vertex_coord,
696 vertex_coord[0][0]=0.0;
697 vertex_coord[0][1]=1.0;
699 vertex_coord[1][1]=1.0;
702 boundary_polyline_pt[1] =
new TriangleMeshPolyLine(vertex_coord,
703 Upper_wall_boundary_id);
707 vertex_coord[0][1]=1.0;
709 vertex_coord[1][1]=0.0;
712 boundary_polyline_pt[2] =
new TriangleMeshPolyLine(vertex_coord,
713 Outflow_boundary_id);
717 vertex_coord[0][1]=0.0;
718 vertex_coord[1][0]=0.0;
719 vertex_coord[1][1]=0.0;
722 boundary_polyline_pt[3] =
new TriangleMeshPolyLine(vertex_coord,
723 Bottom_wall_boundary_id);
726 Outer_boundary_polyline_pt =
new TriangleMeshPolygon(boundary_polyline_pt);
733 Bubble_polygon_pt.resize(1);
737 double y_center = 0.5;
742 Vector<double> zeta(1);
745 Vector<double> coord(2);
748 unsigned npoints = 16;
749 double unit_zeta = MathematicalConstants::Pi/double(npoints-1);
753 Vector<TriangleMeshCurveSection*> bubble_polyline_pt(2);
756 Vector<Vector<double> > bubble_vertex(npoints);
759 for(
unsigned ipoint=0; ipoint<npoints;ipoint++)
762 bubble_vertex[ipoint].resize(2);
765 zeta[0]=unit_zeta*double(ipoint);
766 bubble_pt->position(zeta,coord);
769 bubble_vertex[ipoint][0]=coord[0]+x_center;
770 bubble_vertex[ipoint][1]=coord[1]+y_center;
774 bubble_polyline_pt[0] =
new TriangleMeshPolyLine(bubble_vertex,
775 First_bubble_boundary_id);
778 for(
unsigned ipoint=0; ipoint<npoints;ipoint++)
781 bubble_vertex[ipoint].resize(2);
784 zeta[0]=(unit_zeta*double(ipoint))+MathematicalConstants::Pi;
785 bubble_pt->position(zeta,coord);
788 bubble_vertex[ipoint][0]=coord[0]+x_center;
789 bubble_vertex[ipoint][1]=coord[1]+y_center;
793 bubble_polyline_pt[1] =
new TriangleMeshPolyLine(bubble_vertex,
794 Second_bubble_boundary_id);
798 Vector<double> bubble_center(2);
799 bubble_center[0]=x_center;
800 bubble_center[1]=y_center;
804 Bubble_polygon_pt[0] =
new TriangleMeshPolygon(
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++)
819 bubble_closed_curve_pt[i]=Bubble_polygon_pt[i];
823 double uniform_element_area=0.2;
827 TriangleMeshParameters triangle_mesh_parameters(
828 outer_closed_curve_pt);
831 triangle_mesh_parameters.internal_closed_curve_pt() =
832 bubble_closed_curve_pt;
835 triangle_mesh_parameters.element_area() =
836 uniform_element_area;
840 new RefineableSolidTriangleMesh<ELEMENT>(
841 triangle_mesh_parameters, this->time_stepper_pt());
844 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
845 Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
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;
854 if (CommandLineArgs::command_line_flag_has_been_set(
"--validation"))
856 Fluid_mesh_pt->min_element_size()=0.01;
860 this->Fluid_mesh_pt->output_boundaries(
"boundaries.dat");
861 this->Fluid_mesh_pt->output(
"mesh.dat");
864 complete_problem_setup();
867 Free_surface_mesh_pt=
new Mesh;
868 create_free_surface_elements();
871 Volume_constraint_mesh_pt =
new Mesh;
872 create_volume_constraint_elements();
878 this->add_sub_mesh(this->Volume_constraint_mesh_pt);
881 this->add_sub_mesh(Fluid_mesh_pt);
884 this->add_sub_mesh(this->Free_surface_mesh_pt);
887 this->build_global_mesh();
890 cout <<
"Number of equations: " << this->assign_eqn_numbers() << std::endl;
899 template<
class ELEMENT>
906 unsigned p_traded_index=Vol_constraint_el_pt->index_of_traded_pressure();
909 unsigned nb=Fluid_mesh_pt->nboundary();
910 for(
unsigned b=First_bubble_boundary_id;b<nb;b++)
913 unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
916 for(
unsigned e=0;e<n_element;e++)
920 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
921 Fluid_mesh_pt->boundary_element_pt(b,e));
924 int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
927 ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
928 new ElasticLineFluidInterfaceElement<ELEMENT>(
929 bulk_elem_pt,face_index);
932 Free_surface_mesh_pt->add_element_pt(el_pt);
935 el_pt->set_boundary_number_in_bulk_mesh(b);
943 el_pt->set_external_pressure_data(
944 Vol_constraint_el_pt->p_traded_data_pt(),p_traded_index);
957 template<
class ELEMENT>
962 Volume_constraint_mesh_pt->add_element_pt(Vol_constraint_el_pt);
965 unsigned nb=Fluid_mesh_pt->nboundary();
966 for(
unsigned b=First_bubble_boundary_id;b<nb;b++)
969 unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
972 for(
unsigned e=0;e<n_element;e++)
976 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
977 Fluid_mesh_pt->boundary_element_pt(b,e));
980 int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
983 ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
984 new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
985 bulk_elem_pt,face_index);
988 el_pt->set_volume_constraint_element(Vol_constraint_el_pt);
991 Volume_constraint_mesh_pt->add_element_pt(el_pt);
1002 template<
class ELEMENT>
1006 map<unsigned,bool> is_on_bubble_bound;
1009 unsigned nbubble=Bubble_polygon_pt.size();
1010 for(
unsigned ibubble=0;ibubble<nbubble;ibubble++)
1014 Vector<unsigned> bubble_bound_id=this->Bubble_polygon_pt[ibubble]->
1015 polygon_boundary_id();
1018 unsigned nbound=bubble_bound_id.size();
1021 for(
unsigned ibound=0;ibound<nbound;ibound++)
1024 unsigned bound_id=bubble_bound_id[ibound];
1027 is_on_bubble_bound[bound_id]=
true;
1034 unsigned nbound=Fluid_mesh_pt->nboundary();
1035 for(
unsigned ibound=0;ibound<nbound;ibound++)
1037 unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
1038 for (
unsigned inod=0;inod<num_nod;inod++)
1041 Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1044 if((ibound==0) || (ibound==1) || (ibound==3))
1051 if(ibound==2) {nod_pt->pin(1);}
1055 SolidNode* solid_node_pt =
dynamic_cast<SolidNode*
>(nod_pt);
1056 if(is_on_bubble_bound[ibound])
1058 solid_node_pt->unpin_position(0);
1059 solid_node_pt->unpin_position(1);
1063 solid_node_pt->pin_position(0);
1064 solid_node_pt->pin_position(1);
1074 unsigned n_element = Fluid_mesh_pt->nelement();
1075 for(
unsigned e=0;e<n_element;e++)
1078 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(Fluid_mesh_pt->element_pt(e));
1098 nbound=this->Fluid_mesh_pt->nboundary();
1099 for(
unsigned ibound=0;ibound<nbound;++ibound)
1101 if ((ibound==Upper_wall_boundary_id)||
1102 (ibound==Bottom_wall_boundary_id)||
1103 (ibound==Outflow_boundary_id)||
1104 (ibound==Inflow_boundary_id))
1107 unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
1108 for (
unsigned inod=0;inod<num_nod;inod++)
1111 Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1114 unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
1117 for (
unsigned t=0;t<=n_prev;t++)
1119 if (ibound!=Inflow_boundary_id)
1122 if (ibound!=Outflow_boundary_id)
1124 nod_pt->set_value(t,0,0.0);
1126 nod_pt->set_value(t,1,0.0);
1130 nod_pt->x(t,0)=nod_pt->x(0,0);
1131 nod_pt->x(t,1)=nod_pt->x(0,1);
1138 unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(Inflow_boundary_id);
1139 for (
unsigned inod=0;inod<num_nod;inod++)
1142 Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(Inflow_boundary_id,
1145 double y = nod_pt->x(1);
1155 template<
class ELEMENT>
1163 sprintf(filename,
"%s/soln%i.dat",
1174 compute_error_estimate(max_err,min_err);
1177 double square_of_l2_norm=0.0;
1178 unsigned nel=Fluid_mesh_pt->nelement();
1179 for (
unsigned e=0;e<nel;e++)
1182 dynamic_cast<ELEMENT*
>(this->Fluid_mesh_pt->element_pt(e))->
1183 square_of_l2_norm();
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 "
1194 << comment <<
"\"\n";
1200 sprintf(filename,
"%s/boundaries%i.dat",
1203 some_file.open(filename);
1204 this->Fluid_mesh_pt->output_boundaries(some_file);
1210 Fluid_mesh_pt->max_and_min_element_size(max_area, min_area);
1218 << this->time_pt()->time() <<
" "
1219 << Fluid_mesh_pt->nelement() <<
" "
1224 << sqrt(square_of_l2_norm) <<
" "
1236 template<
class ELEMENT>
1241 ErrorEstimator* err_est_pt=Fluid_mesh_pt->spatial_error_estimator_pt();
1244 unsigned nel=Fluid_mesh_pt->nelement();
1245 Vector<double> elemental_error(nel);
1250 Mesh* fluid_mesh_pt=
dynamic_cast<Mesh*
>(Fluid_mesh_pt);
1251 err_est_pt->get_element_errors(fluid_mesh_pt,
1257 for (
unsigned e=0;e<nel;e++)
1260 set_error(elemental_error[e]);
1262 max_err=std::max(max_err,elemental_error[e]);
1263 min_err=std::min(min_err,elemental_error[e]);
1276 CommandLineArgs::setup(argc,argv);
1282 CommandLineArgs::specify_command_line_flag(
"--validation");
1285 CommandLineArgs::parse_and_assign();
1288 CommandLineArgs::doc_specified_flags();
1318 problem.steady_newton_solve(1);
1326 problem.initialise_dt(dt);
1329 problem.assign_initial_values_impulsive();
1342 if (CommandLineArgs::command_line_flag_has_been_set(
"--validation"))
1345 oomph_info <<
"Remeshing after every second step during validation\n";
1347 for (
unsigned i=0;i<nstep;i++)
1350 problem.unsteady_newton_solve(dt);
1357 unsigned ncycle=1000;
1358 if (CommandLineArgs::command_line_flag_has_been_set(
"--validation"))
1361 oomph_info <<
"Only doing one cycle during validation\n";
1365 for(
unsigned j=0;j<ncycle;j++)
1368 unsigned max_adapt=1;
1371 for (
unsigned i=0;i<nstep;i++)
1374 problem.unsteady_newton_solve(dt,max_adapt,
false);
1378 std::stringstream label;
1379 label <<
"Adaptation " <<j <<
" Step "<< i;
int main(int argc, char **argv)
Driver code for moving bubble problem.
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
BubbleInChannelProblem()
Constructor.
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.
~BubbleInChannelProblem()
Destructor.
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.