35 #include "navier_stokes.h" 38 #include "meshes/tube_mesh.h" 42 using namespace oomph;
54 GeomObject(3,3), Radius(radius) { }
60 void position (
const Vector<double>& xi, Vector<double>& r)
const 62 r[0] = xi[2]*Radius*cos(xi[1]);
64 r[2] = -xi[2]*Radius*sin(xi[1]);
71 const Vector<double>& xi, Vector<double>& r)
const 97 template<
class ELEMENT>
105 const double& max_error_target);
111 void actions_before_newton_solve();
117 RefineableNavierStokesEquations<3>::
118 pin_redundant_nodal_pressures(mesh_pt()->element_pt());
128 return dynamic_cast<RefineableTubeMesh<ELEMENT>*
>(Problem::mesh_pt());
147 template<
class ELEMENT>
149 const double& min_error_target,
150 const double& max_error_target)
162 const double pi = MathematicalConstants::Pi;
166 Vector<double> centreline_limits(2);
167 centreline_limits[0] = 0.0;
168 centreline_limits[1] = pi;
173 Vector<double> theta_positions(4);
174 theta_positions[0] = -0.75*pi;
175 theta_positions[1] = -0.25*pi;
176 theta_positions[2] = 0.25*pi;
177 theta_positions[3] = 0.75*pi;
181 Vector<double> radial_frac(4,0.5);
187 Problem::mesh_pt()=
new RefineableTubeMesh<ELEMENT>(
Volume_pt,
194 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
195 mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
198 mesh_pt()->max_permitted_error()=max_error_target;
199 mesh_pt()->min_permitted_error()=min_error_target;
206 ELEMENT::Gamma[0] = 0.0;
207 ELEMENT::Gamma[1] = 0.0;
208 ELEMENT::Gamma[2] = 0.0;
211 unsigned num_bound =
mesh_pt()->nboundary();
212 for(
unsigned ibound=0;ibound<num_bound;ibound++)
214 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
215 for (
unsigned inod=0;inod<num_nod;inod++)
220 if((ibound==0) || (ibound==1))
222 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
223 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
224 mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
231 unsigned n_element =
mesh_pt()->nelement();
232 for(
unsigned i=0;i<n_element;i++)
235 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(
mesh_pt()->element_pt(i));
242 RefineableNavierStokesEquations<3>::
243 pin_redundant_nodal_pressures(
mesh_pt()->element_pt());
246 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
254 template<
class ELEMENT>
261 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
262 for (
unsigned inod=0;inod<num_nod;inod++)
265 double x=
mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
266 double z=
mesh_pt()->boundary_node_pt(ibound,inod)->x(2);
268 double r=sqrt(x*x+z*z);
271 mesh_pt()->boundary_node_pt(ibound,inod)->
272 set_value(1,(1.0-pow(r,10.0)));
281 template<
class ELEMENT>
295 sprintf(filename,
"%s/soln_Re%g.dat",
Doc_info.directory().c_str(),
297 some_file.open(filename);
298 mesh_pt()->output(some_file,npts);
314 int main(
int argc,
char* argv[])
319 unsigned max_adapt=1;
320 double max_error_target=0.02 , min_error_target=0.002;
330 doc_info.set_directory(
"RESLT_TH");
337 problem(doc_info,min_error_target,max_error_target);
339 cout <<
" Doing Taylor-Hood elements " << std::endl;
341 for(
unsigned i=0;i<2;i++)
344 problem.newton_solve(max_adapt);
358 doc_info.set_directory(
"RESLT_CR");
365 problem(doc_info,min_error_target,max_error_target);
367 cout <<
" Doing Crouzeix-Raviart elements " << std::endl;
369 for(
unsigned i=0;i<2;i++)
372 problem.newton_solve(max_adapt);
int main(int argc, char *argv[])
Driver for 3D entry flow into a straight tube.
~SteadyTubeProblem()
Destructor (empty)
MyCylinder(const double &radius)
Constructor that takes the radius of the tube as its argument.
double Radius
Storage for the radius of the tube.
virtual ~MyCylinder()
Destructor.
Entry flow problem in tapered tube domain.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
SteadyTubeProblem(DocInfo &doc_info, const double &min_error_target, const double &max_error_target)
Constructor: Pass DocInfo object and target errors.
double Re
Reynolds number.
void doc_solution()
Doc the solution.
GeomObject * Volume_pt
Pointer to GeomObject that specifies the domain volume.
RefineableTubeMesh< ELEMENT > * mesh_pt()
Overload generic access function by one that returns a pointer to the specific mesh.
Namespace for physical parameters.
void actions_after_adapt()
After adaptation: Pin redudant pressure dofs.
void actions_before_newton_solve()
Update the problem specs before solve.
DocInfo Doc_info
Doc info object.
void position(const Vector< double > &xi, Vector< double > &r) const
Lagrangian coordinate xi.