30 #ifndef OOMPH_QUARTER_TUBE_MESH_TEMPLATE_CC 31 #define OOMPH_QUARTER_TUBE_MESH_TEMPLATE_CC 49 template <
class ELEMENT>
52 const double& fract_mid,
54 const unsigned& nlayer,
56 Wall_pt(wall_pt), Xi_lo(xi_lo), Fract_mid(fract_mid), Xi_hi(xi_hi)
60 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3);
72 unsigned nelem=3*nlayer;
76 ELEMENT* dummy_el_pt=
new ELEMENT;
79 unsigned n_p = dummy_el_pt->nnode_1d();
85 unsigned nnodes_total=
86 (n_p*n_p+(n_p-1)*n_p+(n_p-1)*(n_p-1))*(1+nlayer*(n_p-1));
97 for (
unsigned ielem=0;ielem<nelem;ielem++)
104 for (
unsigned i2=0;i2<n_p;i2++)
108 for (
unsigned i1=0;i1<n_p;i1++)
112 for (
unsigned i0=0;i0<n_p;i0++)
116 unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
120 construct_node(jnod_local,time_stepper_pt);
123 s[0]=-1.0+2.0*double(i0)/double(n_p-1);
124 s[1]=-1.0+2.0*double(i1)/double(n_p-1);
125 s[2]=-1.0+2.0*double(i2)/double(n_p-1);
126 Domain_pt->macro_element_pt(ielem)->macro_map(s,r);
128 node_pt->
x(0) = r[0];
129 node_pt->
x(1) = r[1];
130 node_pt->
x(2) = r[2];
139 unsigned node_count=0;
142 double node_kill_tol=1.0e-12;
148 for (
unsigned ielem=0;ielem<nelem;ielem++)
152 unsigned ilayer=unsigned(ielem/3);
164 for (
unsigned i2=0;i2<n_p;i2++)
168 for (
unsigned i1=0;i1<n_p;i1++)
172 for (
unsigned i0=0;i0<n_p;i0++)
176 unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
183 if ((i2==0)&&(ilayer>0))
187 unsigned ielem_neigh=ielem-3;
190 unsigned i0_neigh=i0;
191 unsigned i1_neigh=i1;
192 unsigned i2_neigh=n_p-1;
193 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
196 for (
unsigned i=0;
i<3;
i++)
198 double error=std::fabs(
202 if (error>node_kill_tol)
205 <<
i <<
" " << error << std::endl;
228 if ((i2==0)&&(ilayer==0))
235 if ((i2==n_p-1)&&(ilayer==nlayer-1))
273 for (
unsigned i2=0;i2<n_p;i2++)
277 for (
unsigned i1=0;i1<n_p;i1++)
281 for (
unsigned i0=0;i0<n_p;i0++)
285 unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
292 if ((i2==0)&&(ilayer>0))
296 unsigned ielem_neigh=ielem-3;
299 unsigned i0_neigh=i0;
300 unsigned i1_neigh=i1;
301 unsigned i2_neigh=n_p-1;
302 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
306 for (
unsigned i=0;
i<3;
i++)
308 double error=std::fabs(
312 if (error>node_kill_tol)
315 <<
i <<
" " << error << std::endl;
338 unsigned ielem_neigh=ielem-1;
341 unsigned i0_neigh=n_p-1;
342 unsigned i1_neigh=i1;
343 unsigned i2_neigh=i2;
344 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
348 for (
unsigned i=0;
i<3;
i++)
350 double error=std::fabs(
354 if (error>node_kill_tol)
357 <<
i <<
" " << error << std::endl;
381 if ((i2==0)&&(ilayer==0))
388 if ((i2==n_p-1)&&(ilayer==nlayer-1))
410 (double(ilayer)+double(i2)/double(n_p-1))*
415 double(i1)/double(n_p-1)*0.5*(
Xi_hi[1]-
Xi_lo[1]);
417 Node_pt[node_count]->set_coordinates_on_boundary(3,zeta);
437 for (
unsigned i2=0;i2<n_p;i2++)
441 for (
unsigned i1=0;i1<n_p;i1++)
445 for (
unsigned i0=0;i0<n_p;i0++)
449 unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
456 if ((i2==0)&&(ilayer>0))
460 unsigned ielem_neigh=ielem-3;
463 unsigned i0_neigh=i0;
464 unsigned i1_neigh=i1;
465 unsigned i2_neigh=n_p-1;
466 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
469 for (
unsigned i=0;
i<3;
i++)
471 double error=std::fabs(
475 if (error>node_kill_tol)
478 <<
i <<
" " << error << std::endl;
502 unsigned ielem_neigh=ielem-1;
505 unsigned i0_neigh=i1;
506 unsigned i1_neigh=n_p-1;
507 unsigned i2_neigh=i2;
508 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
511 for (
unsigned i=0;
i<3;
i++)
513 double error=std::fabs(
517 if (error>node_kill_tol)
520 <<
i <<
" " << error << std::endl;
537 if ((i1==0)&&(i0!=n_p-1))
541 unsigned ielem_neigh=ielem-2;
544 unsigned i0_neigh=i0;
545 unsigned i1_neigh=n_p-1;
546 unsigned i2_neigh=i2;
547 unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
550 for (
unsigned i=0;
i<3;
i++)
552 double error=std::fabs(
556 if (error>node_kill_tol)
559 <<
i <<
" " << error << std::endl;
583 if ((i2==0)&&(ilayer==0))
590 if ((i2==n_p-1)&&(ilayer==nlayer-1))
613 (double(ilayer)+double(i2)/double(n_p-1))*
618 double(i0)/double(n_p-1)*0.5*(
Xi_hi[1]-
Xi_lo[1]);
620 Node_pt[node_count]->set_coordinates_on_boundary(3,zeta);
641 std::ostringstream error_stream;
642 error_stream <<
"Error in killing nodes\n" 643 <<
"The most probable cause is that the domain is not\n" 644 <<
"compatible with the mesh.\n" 645 <<
"For the QuarterTubeMesh, the domain must be\n" 646 <<
"topologically consistent with a quarter tube with a\n" 647 <<
"non-curved centreline.\n";
649 OOMPH_CURRENT_FUNCTION,
650 OOMPH_EXCEPTION_LOCATION);
689 template <
class ELEMENT>
699 if (algebraic_element_pt==0)
701 std::ostringstream error_message;
703 "Element in AlgebraicRefineableQuarterTubeMesh must be\n ";
704 error_message <<
"derived from AlgebraicElementBase\n";
705 error_message<<
"but it is of type: " 708 "AlgebraicRefineableQuarterTubeMesh::";
709 function_name +=
"setup_algebraic_node_update()";
711 OOMPH_CURRENT_FUNCTION,
712 OOMPH_EXCEPTION_LOCATION);
721 unsigned nnodes_2d = nnodes_1d*nnodes_1d;
725 unsigned tl_node = nnodes_2d-nnodes_1d;
726 unsigned br_node = nnodes_1d-1;
742 Lambda_x = Centre_box_size*x_c_element/x_wall;
743 Lambda_y = Centre_box_size*y_c_element/y_wall;
749 for (
unsigned e=0;
e<nelements;
e++)
755 unsigned region_id =
e%3;
759 unsigned nstart=nnodes_2d;
766 for (
unsigned n=nstart; n<nnodes_3d; n++)
804 geom_object_pt[0] = obj_br_pt;
807 geom_object_pt[1] = obj_tl_pt;
813 ref_value[0]=x/x_c_element;
816 ref_value[1]=y/y_c_element;
827 ref_value[3]=s_br[0];
828 ref_value[4]=s_br[1];
832 ref_value[5]=s_tl[0];
833 ref_value[6]=s_tl[1];
837 add_node_update_info(
846 else if (region_id==1)
849 double fract = 1.0/double(nnodes_1d-1);
852 double rho_0 = fract * double(n%nnodes_1d);
855 double rho_1 = fract * double((n/nnodes_1d)%nnodes_1d);
873 geom_object_pt[0] = obj_br_pt;
876 geom_object_pt[1] = obj_tl_pt;
879 geom_object_pt[2] = obj_wall_pt;
895 ref_value[3]=s_br[0];
896 ref_value[4]=s_br[1];
900 ref_value[5]=s_tl[0];
901 ref_value[6]=s_tl[1];
905 ref_value[7]=s_wall[0];
906 ref_value[8]=s_wall[1];
910 add_node_update_info(
919 else if (region_id==2)
922 double fract = 1.0/double(nnodes_1d-1);
925 double rho_0 = fract * double(n%nnodes_1d);
928 double rho_1 = fract * double((n/nnodes_1d)%nnodes_1d);
946 geom_object_pt[0] = obj_br_pt;
949 geom_object_pt[1] = obj_tl_pt;
952 geom_object_pt[2] = obj_wall_pt;
968 ref_value[3]=s_br[0];
969 ref_value[4]=s_br[1];
973 ref_value[5]=s_tl[0];
974 ref_value[6]=s_tl[1];
978 ref_value[7]=s_wall[0];
979 ref_value[8]=s_wall[1];
983 add_node_update_info(
998 template <
class ELEMENT>
1020 "Trying to update the nodal position at a time level\n";
1021 error_message +=
"beyond the number of previous values in the nodes'\n";
1022 error_message +=
"position timestepper. This seems highly suspect!\n";
1023 error_message +=
"If you're sure the code behaves correctly\n";
1024 error_message +=
"in your application, remove this warning \n";
1025 error_message +=
"or recompile with PARNOID switched off.\n";
1028 "AlgebraicRefineableQuarterTubeMesh::";
1029 function_name +=
"node_update_central_region()",
1031 OOMPH_CURRENT_FUNCTION,
1032 OOMPH_EXCEPTION_LOCATION);
1041 vector_geom_object_pt(Central_region));
1044 double rho_x=ref_value[0];
1047 double rho_y=ref_value[1];
1056 s_br[0]=ref_value[3];
1057 s_br[1]=ref_value[4];
1060 unsigned n_dim = obj_br_pt->
ndim();
1071 s_tl[0]=ref_value[5];
1072 s_tl[1]=ref_value[6];
1079 node_pt->
x(t,0)=r_br[0]*Lambda_x*rho_x;
1080 node_pt->
x(t,1)=r_tl[1]*Lambda_y*rho_y;
1081 node_pt->
x(t,2)=(r_br[2]+r_tl[2])*0.5;
1089 template <
class ELEMENT>
1111 "Trying to update the nodal position at a time level";
1112 error_message +=
"beyond the number of previous values in the nodes'";
1113 error_message +=
"position timestepper. This seems highly suspect!";
1114 error_message +=
"If you're sure the code behaves correctly";
1115 error_message +=
"in your application, remove this warning ";
1116 error_message +=
"or recompile with PARNOID switched off.";
1119 "AlgebraicRefineableQuarterTubeSectorMesh::";
1120 function_name +=
"node_update_lower_right_region()",
1122 OOMPH_CURRENT_FUNCTION,
1123 OOMPH_EXCEPTION_LOCATION);
1136 double rho_0=ref_value[0];
1142 double rho_1=ref_value[1];
1151 s_br[0]=ref_value[3];
1152 s_br[1]=ref_value[4];
1155 unsigned n_dim=obj_br_pt->
ndim();
1168 s_tl[0]=ref_value[5];
1169 s_tl[1]=ref_value[6];
1181 s_wall[0]=ref_value[7];
1182 s_wall[1]=ref_value[8];
1185 obj_wall_pt->
position(t,s_wall,r_wall);
1189 r_corner[0]=Lambda_x*r_br[0];
1190 r_corner[1]=Lambda_y*r_tl[1];
1191 r_corner[2]=(r_br[2]+r_tl[2])*0.5;
1195 r_left[0]=r_corner[0];
1196 r_left[1]=rho_1*r_corner[1];
1197 r_left[2]=r_corner[2];
1200 node_pt->
x(t,0)=r_left[0]+rho_0*(r_wall[0]-r_left[0]);
1201 node_pt->
x(t,1)=r_left[1]+rho_0*(r_wall[1]-r_left[1]);
1202 node_pt->
x(t,2)=r_left[2]+rho_0*(r_wall[2]-r_left[2]);
1209 template <
class ELEMENT>
1231 "Trying to update the nodal position at a time level";
1232 error_message +=
"beyond the number of previous values in the nodes'";
1233 error_message +=
"position timestepper. This seems highly suspect!";
1234 error_message +=
"If you're sure the code behaves correctly";
1235 error_message +=
"in your application, remove this warning ";
1236 error_message +=
"or recompile with PARNOID switched off.";
1239 "AlgebraicRefineableQuarterTubeMesh::";
1240 function_name +=
"node_update_upper_left_region()";
1243 OOMPH_CURRENT_FUNCTION,
1244 OOMPH_EXCEPTION_LOCATION);
1257 double rho_0=ref_value[0];
1260 double rho_1=ref_value[1];
1271 unsigned n_dim=obj_br_pt->
ndim();
1275 s_br[0]=ref_value[3];
1276 s_br[1]=ref_value[4];
1289 s_tl[0]=ref_value[5];
1290 s_tl[1]=ref_value[6];
1302 s_wall[0]=ref_value[7];
1303 s_wall[1]=ref_value[8];
1306 obj_wall_pt->
position(t,s_wall,r_wall);
1310 r_corner[0]=Lambda_x*r_br[0];
1311 r_corner[1]=Lambda_y*r_tl[1];
1312 r_corner[2]=(r_br[2]+r_tl[2])*0.5;
1316 r_top[0]=rho_0*r_corner[0];
1317 r_top[1]=r_corner[1];
1318 r_top[2]=r_corner[2];
1321 node_pt->
x(t,0)=r_top[0]+rho_1*(r_wall[0]-r_top[0]);
1322 node_pt->
x(t,1)=r_top[1]+rho_1*(r_wall[1]-r_top[1]);
1323 node_pt->
x(t,2)=r_top[2]+rho_1*(r_wall[2]-r_top[2]);
1332 template <
class ELEMENT>
1357 double z = ref_value[2];
1373 geom_object_pt[0] = obj_br_pt;
1376 ref_value[3]=s_br[0];
1377 ref_value[4]=s_br[1];
1392 geom_object_pt[1] = obj_tl_pt;
1395 ref_value[5]=s_tl[0];
1396 ref_value[6]=s_tl[1];
1398 if (region_id!=Central_region)
1404 if (region_id==Lower_right_region)
1407 double rho_1=ref_value[1];
1410 xi[1]=xi_lo+rho_1*fract_mid*(xi_hi-xi_lo);
1415 double rho_0=ref_value[0];
1418 xi[1]=xi_hi-rho_0*(1.0-fract_mid)*(xi_hi-xi_lo);
1429 geom_object_pt[2] = obj_wall_pt;
1432 ref_value[7]=s_wall[0];
1433 ref_value[8]=s_wall[1];
void setup_boundary_element_info()
void node_update_lower_right_region(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for a node that is located in the lower-right region.
GeomObject * geom_object_pt(const unsigned &i)
Return pointer to i-th geometric object involved in default (usually first) update function...
Vector< Node * > Node_pt
Vector of pointers to nodes.
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement *> &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this...
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
void kill_node_update_info(const int &id=0)
Erase algebraic node update information for id-th node update function. Id defaults to 0...
void add_node_update_info(const int &id, AlgebraicMesh *mesh_pt, const Vector< GeomObject *> &geom_object_pt, const Vector< double > &ref_value, const bool &called_from_constructor=false)
Add algebraic update information for node: What's the ID of the mesh update function (typically used ...
A general Finite Element class.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
QuarterTubeMesh(GeomObject *wall_pt, const Vector< double > &xi_lo, const double &fract_mid, const Vector< double > &xi_hi, const unsigned &nlayer, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the wall, start and end coordinates on t...
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
void update_node_update_in_region(AlgebraicNode *&node_pt, int ®ion_id)
Update algebraic node update function for nodes in the region defined by region_id.
Vector< double > Xi_hi
Upper limits for the coordinates along the wall.
unsigned long nelement() const
Return number of elements in the mesh.
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
Vector< double > Xi_lo
Lower limits for the coordinates along the wall.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
3D quarter tube mesh class. The domain is specified by the GeomObject that identifies boundary 3...
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
void node_update_central_region(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for a node that is located in the central region.
GeomObject * Wall_pt
Pointer to the geometric object that represents the curved wall.
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Quarter tube as domain. Domain is bounded by curved boundary which is represented by a GeomObject...
double s_squashed(const double &s)
Function that squashes the outer two macro elements towards the wall by mapping the input value of th...
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
unsigned ndim() const
Access function to # of Eulerian coordinates.
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overloa...
void node_update_upper_left_region(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for a node that is located in the upper-left region.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
QuarterTubeDomain * Domain_pt
Pointer to domain.
virtual void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
A geometric object may be composed of may sub-objects (e.g. a finite-element representation of a boun...
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
unsigned nnode() const
Return the number of nodes.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
void setup_algebraic_node_update()
Setup algebraic update operation for all nodes.
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.