30 #include "navier_stokes.h"
33 #include "meshes/quarter_tube_mesh.h"
37 using namespace oomph;
52 template<
class ELEMENT>
60 const double& max_error_target);
76 void actions_before_newton_solve();
82 RefineableNavierStokesEquations<3>::
83 pin_redundant_nodal_pressures(mesh_pt()->element_pt());
91 RefineableQuarterTubeMesh<ELEMENT>*
mesh_pt()
93 return dynamic_cast<RefineableQuarterTubeMesh<ELEMENT>*
>(Problem::mesh_pt());
112 template<
class ELEMENT>
114 const double& min_error_target,
115 const double& max_error_target)
124 GeomObject* Wall_pt=
new EllipticalTube(radius,radius);
127 Vector<double> xi_lo(2);
133 Vector<double> xi_hi(2);
137 xi_hi[1]=2.0*atan(1.0);
147 new RefineableQuarterTubeMesh<ELEMENT>(Wall_pt,xi_lo,frac_mid,xi_hi,nlayer);
151 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
152 mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
155 mesh_pt()->max_permitted_error()=max_error_target;
156 mesh_pt()->min_permitted_error()=min_error_target;
163 sprintf(filename,
"boundaries.dat");
164 some_file.open(filename);
165 mesh_pt()->output_boundaries(some_file);
172 unsigned num_bound =
mesh_pt()->nboundary();
173 for(
unsigned ibound=0;ibound<num_bound;ibound++)
175 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
176 for (
unsigned inod=0;inod<num_nod;inod++)
180 if(ibound!=1)
mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
184 if(ibound!=2)
mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
189 if((ibound==0) || (ibound==3))
191 mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
198 unsigned n_element =
mesh_pt()->nelement();
199 for(
unsigned i=0;i<n_element;i++)
202 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(
mesh_pt()->element_pt(i));
209 RefineableNavierStokesEquations<3>::
210 pin_redundant_nodal_pressures(
mesh_pt()->element_pt());
214 unsigned n_nod=
mesh_pt()->nnode();
215 for (
unsigned j=0;j<n_nod;j++)
217 Node* node_pt=
mesh_pt()->node_pt(j);
219 double x=node_pt->x(0);
220 double y=node_pt->x(1);
221 double r=sqrt(x*x+y*y );
224 node_pt->set_value(0,0.0);
225 node_pt->set_value(1,0.0);
226 node_pt->set_value(2,(1.0-r*r));
234 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
242 template<
class ELEMENT>
251 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
252 for (
unsigned inod=0;inod<num_nod;inod++)
255 double x=mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
256 double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
257 double r=sqrt(x*x+y*y);
260 mesh_pt()->boundary_node_pt(ibound,inod)->
261 set_value(2,(1.0-pow(r,Alpha)));
270 template<
class ELEMENT>
282 sprintf(filename,
"%s/soln%i.dat",Doc_info.directory().c_str(),
284 some_file.open(filename);
285 mesh_pt()->output(some_file,npts);
303 int main(
int argc,
char* argv[])
307 CommandLineArgs::setup(argc,argv);
314 double max_error_target,min_error_target;
318 if (CommandLineArgs::Argc==1)
324 max_error_target=0.005;
325 min_error_target=0.0005;
336 max_error_target=0.02;
337 min_error_target=0.002;
349 doc_info.set_directory(
"RESLT_TH");
356 problem(doc_info,min_error_target,max_error_target);
358 cout <<
" Doing Taylor-Hood elements " << std::endl;
365 problem.newton_solve(max_adapt);
373 doc_info.set_directory(
"RESLT_CR");
380 problem(doc_info,min_error_target,max_error_target);
382 cout <<
" Doing Crouzeix-Raviart elements " << std::endl;
385 problem.newton_solve(max_adapt);
Entry flow problem in quarter tube domain.
void actions_after_newton_solve()
Doc the solution after solve.
EntryFlowProblem(DocInfo &doc_info, const double &min_error_target, const double &max_error_target)
Constructor: Pass DocInfo object and target errors.
~EntryFlowProblem()
Destructor (empty)
RefineableQuarterTubeMesh< ELEMENT > * mesh_pt()
Overload generic access function by one that returns a pointer to the specific mesh.
int Alpha
Exponent for bluntness of velocity profile.
void doc_solution()
Doc the solution.
DocInfo Doc_info
Doc info object.
void actions_after_adapt()
After adaptation: Pin redudant pressure dofs.
void actions_before_newton_solve()
Update the problem specs before solve.
Namespace for physical parameters.
double Re
Reynolds number.
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...