43 template<
class ELEMENT>
44 Vector<std::pair<unsigned, int>> ExtrudedCubeMeshFromQuadMesh<
45 ELEMENT>::get_element_boundary_information(
QuadMeshBase* quad_mesh_pt,
52 unsigned n_boundary = quad_mesh_pt->
nboundary();
55 for (
unsigned b = 0; b < n_boundary; b++)
61 for (
unsigned e = 0;
e < n_boundary_element;
e++)
67 std::pair<unsigned, int> boundary_and_face_index;
70 boundary_and_face_index.first = b;
73 boundary_and_face_index.second =
77 boundary_info.push_back(boundary_and_face_index);
98 template<
class ELEMENT>
103 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
109 std::ostringstream error_message;
112 error_message <<
"Extrude2DQuadMeshTo3DCubeMesh needs at least two "
113 <<
"elements in each,\ncoordinate direction. You have "
114 <<
"specified Nz=" << Nz << std::endl;
118 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
128 unsigned n_boundary_in_quad_mesh = quad_mesh_pt->
nboundary();
131 unsigned n_boundary_in_extruded_mesh = n_boundary_in_quad_mesh + 2;
135 set_nboundary(n_boundary_in_extruded_mesh);
138 unsigned n_element_in_quad_mesh = quad_mesh_pt->
nelement();
143 n_element_in_quad_mesh);
146 for (
unsigned e = 0;
e < n_element_in_quad_mesh;
e++)
152 boundary_information[
e] =
153 get_element_boundary_information(quad_mesh_pt, quad_el_pt);
157 Element_pt.resize(n_element_in_quad_mesh * Nz);
160 ELEMENT* temp_el_pt =
new ELEMENT;
163 unsigned n_node_1d = temp_el_pt->nnode_1d();
166 N_node_1d = n_node_1d;
179 std::ostringstream error_message;
183 <<
"Extrude2DQuadMeshTo3DCubeMesh needs the number of "
184 <<
"nodes (in the 2D element) in one direction\n to match "
185 <<
"the number of nodes in one direction in the 3D element!";
189 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
193 unsigned long n_node_in_quad_mesh = quad_mesh_pt->
nnode();
196 for (
unsigned i = 0;
i < n_node_in_quad_mesh;
i++)
202 std::ostringstream warning_message_stream;
205 warning_message_stream <<
"Extrusion machinery is not currently able "
206 <<
"to deal with hanging nodes!" << std::endl;
210 OOMPH_CURRENT_FUNCTION,
211 OOMPH_EXCEPTION_LOCATION);
216 unsigned n_node_in_z_direction = (1 + (n_node_1d - 1) * Nz);
219 unsigned long n_node_in_cube_mesh =
220 n_node_in_quad_mesh * n_node_in_z_direction;
223 Node_pt.resize(n_node_in_cube_mesh);
226 unsigned long node_count = 0;
229 for (
unsigned i = 0;
i < Nz;
i++)
235 std::map<std::pair<Edge, double>,
Vector<Node*>> edge_to_nodes_map;
244 std::map<std::pair<Node*, double>,
Node*>
245 quad_corner_node_to_cube_node_map;
248 for (
unsigned e = 0;
e < n_element_in_quad_mesh;
e++)
254 unsigned element_index =
i * n_element_in_quad_mesh +
e;
257 ELEMENT* cube_el_pt =
new ELEMENT;
260 Element_pt[element_index] = cube_el_pt;
263 unsigned quad_local_node_number = 0;
266 unsigned cube_local_node_number = 0;
270 boundary_information[
e];
273 for (
unsigned n_2 = 0; n_2 < n_node_1d; n_2++)
277 double z_value = z_spacing_function(
i, n_2);
297 if ((
i > 0) && (n_2 == 0))
300 unsigned copy_element_index = (
i - 1) * n_element_in_quad_mesh +
e;
305 finite_element_pt(copy_element_index);
309 unsigned copy_cube_local_node_number = 0;
312 for (
unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
315 for (
unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
318 cube_local_node_number = n_0 + (n_node_1d * n_1);
321 quad_local_node_number = n_0 + (n_node_1d * n_1);
325 copy_cube_local_node_number =
326 cube_local_node_number +
327 n_node_1d * n_node_1d * (n_node_1d - 1);
330 cube_el_pt->node_pt(cube_local_node_number) =
331 copy_cube_el_pt->
node_pt(copy_cube_local_node_number);
351 std::vector<bool> has_node_been_setup(n_node_1d * n_node_1d,
false);
361 std::vector<bool> has_edge_been_done(n_edge,
false);
382 node1_pt = quad_el_pt->
node_pt((n_node_1d * n_node_1d) - 1);
385 node2_pt = quad_el_pt->
node_pt(n_node_1d * (n_node_1d - 1));
391 node1_pt = quad_el_pt->
node_pt(n_node_1d - 1);
394 node2_pt = quad_el_pt->
node_pt((n_node_1d * n_node_1d) - 1);
400 node1_pt = quad_el_pt->
node_pt(0);
403 node2_pt = quad_el_pt->
node_pt(n_node_1d - 1);
409 node1_pt = quad_el_pt->
node_pt(n_node_1d * (n_node_1d - 1));
412 node2_pt = quad_el_pt->
node_pt(0);
418 std::ostringstream error_message_stream;
422 <<
"Input is " << i_edge <<
" but it can only\n"
429 OOMPH_CURRENT_FUNCTION,
430 OOMPH_EXCEPTION_LOCATION);
435 Edge edge(node1_pt, node2_pt);
438 std::pair<Edge, double> edge_and_z_value_pair(edge, z_value);
441 edges_and_z_value.push_back(edge_and_z_value_pair);
444 std::map<std::pair<Edge, double>,
Vector<Node*>>::iterator it =
445 edge_to_nodes_map.find(edge_and_z_value_pair);
448 if (it != edge_to_nodes_map.end())
455 if (edge_nodes_pt.size() != n_node_1d)
459 OOMPH_CURRENT_FUNCTION,
460 OOMPH_EXCEPTION_LOCATION);
469 has_edge_been_done[0] =
true;
475 has_edge_been_done[1] =
true;
481 has_edge_been_done[2] =
true;
487 has_edge_been_done[3] =
true;
494 OOMPH_CURRENT_FUNCTION,
495 OOMPH_EXCEPTION_LOCATION);
506 for (
unsigned i_node = 0; i_node < n_node_1d; i_node++)
514 quad_local_node_number =
515 ((n_node_1d * n_node_1d) - 1) - i_node;
521 quad_local_node_number =
522 (n_node_1d - 1) + (i_node * n_node_1d);
528 quad_local_node_number = i_node;
534 quad_local_node_number =
535 (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
542 cube_local_node_number =
543 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
546 cube_el_pt->node_pt(cube_local_node_number) =
547 edge_nodes_pt[(n_node_1d - 1) - i_node];
550 has_node_been_setup[quad_local_node_number] =
true;
556 unsigned n_corner = 4;
559 for (
unsigned i_corner = 0; i_corner < n_corner; i_corner++)
567 quad_local_node_number = 0;
573 quad_local_node_number = n_node_1d - 1;
579 quad_local_node_number = n_node_1d * (n_node_1d - 1);
585 quad_local_node_number = (n_node_1d * n_node_1d) - 1;
592 if (!has_node_been_setup[quad_local_node_number])
595 Node* quad_corner_node_pt =
596 quad_el_pt->
node_pt(quad_local_node_number);
599 std::pair<Node*, double> quad_corner_node_and_z_value(
600 quad_corner_node_pt, z_value);
603 std::map<std::pair<Node*, double>,
Node*>::iterator it =
604 quad_corner_node_to_cube_node_map.find(
605 quad_corner_node_and_z_value);
608 if (it != quad_corner_node_to_cube_node_map.end())
611 cube_local_node_number =
612 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
615 cube_el_pt->node_pt(cube_local_node_number) = it->second;
618 has_node_been_setup[quad_local_node_number] =
true;
628 for (
unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
631 for (
unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
634 quad_local_node_number = n_0 + (n_node_1d * n_1);
637 cube_local_node_number =
638 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
641 if (has_node_been_setup[quad_local_node_number])
648 Node* cube_node_pt = 0;
659 if (((
i == 0) && (n_2 == 0)) ||
660 ((
i == Nz - 1) && (n_2 == n_node_1d - 1)))
663 if ((
i == 0) && (n_2 == 0))
666 cube_node_pt = cube_el_pt->construct_boundary_node(
667 cube_local_node_number, time_stepper_pt);
670 has_node_been_setup[quad_local_node_number] =
true;
673 add_boundary_node(n_boundary_in_quad_mesh, cube_node_pt);
676 else if ((
i == Nz - 1) && (n_2 == n_node_1d - 1))
679 cube_node_pt = cube_el_pt->construct_boundary_node(
680 cube_local_node_number, time_stepper_pt);
683 has_node_been_setup[quad_local_node_number] =
true;
686 add_boundary_node(n_boundary_in_quad_mesh + 1,
692 quad_el_pt->
node_pt(quad_local_node_number);
698 std::set<unsigned>* boundaries_pt;
704 std::set<unsigned>::const_iterator const_it;
707 for (const_it = boundaries_pt->begin();
708 const_it != boundaries_pt->end();
713 add_boundary_node(*const_it, cube_node_pt);
720 if (el_boundary_information.size() > 0)
723 unsigned n_boundaries = el_boundary_information.size();
726 for (
unsigned i_bound = 0; i_bound < n_boundaries; i_bound++)
729 if (!has_node_been_setup[quad_local_node_number])
732 cube_node_pt = cube_el_pt->construct_boundary_node(
733 cube_local_node_number, time_stepper_pt);
736 has_node_been_setup[quad_local_node_number] =
true;
741 std::pair<unsigned, int> boundary_and_face_index =
742 el_boundary_information[i_bound];
746 switch (boundary_and_face_index.second)
755 add_boundary_node(boundary_and_face_index.first,
765 if (n_0 == n_node_1d - 1)
768 add_boundary_node(boundary_and_face_index.first,
781 add_boundary_node(boundary_and_face_index.first,
791 if (n_1 == n_node_1d - 1)
794 add_boundary_node(boundary_and_face_index.first,
811 quad_el_pt->
node_pt(quad_local_node_number);
818 if (cube_node_pt == 0)
821 cube_node_pt = cube_el_pt->construct_boundary_node(
822 cube_local_node_number, time_stepper_pt);
825 has_node_been_setup[quad_local_node_number] =
true;
829 std::set<unsigned>* boundaries_pt;
835 std::set<unsigned>::const_iterator const_it;
838 for (const_it = boundaries_pt->begin();
839 const_it != boundaries_pt->end();
844 add_boundary_node(*const_it, cube_node_pt);
849 if (!has_node_been_setup[quad_local_node_number])
852 cube_node_pt = cube_el_pt->construct_node(
853 cube_local_node_number, time_stepper_pt);
856 has_node_been_setup[quad_local_node_number] =
true;
860 cube_node_pt->
x(0) = quad_node_pt->
x(0);
863 cube_node_pt->
x(1) = quad_node_pt->
x(1);
866 cube_node_pt->
x(2) = z_value;
869 Node_pt[node_count] = cube_node_pt;
880 for (
unsigned i_edge = 0; i_edge < n_edge; i_edge++)
883 if (!has_edge_been_done[i_edge])
889 for (
unsigned i_node = 0; i_node < n_node_1d; i_node++)
897 quad_local_node_number =
898 ((n_node_1d * n_node_1d) - 1) - i_node;
904 quad_local_node_number =
905 (n_node_1d - 1) + (i_node * n_node_1d);
911 quad_local_node_number = i_node;
917 quad_local_node_number =
918 (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
925 cube_local_node_number =
926 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
929 edge_nodes_pt[i_node] =
930 cube_el_pt->node_pt(cube_local_node_number);
934 edge_to_nodes_map[edges_and_z_value[i_edge]] = edge_nodes_pt;
943 for (
unsigned i_corner = 0; i_corner < n_corner; i_corner++)
951 quad_local_node_number = 0;
957 quad_local_node_number = n_node_1d - 1;
963 quad_local_node_number = n_node_1d * (n_node_1d - 1);
969 quad_local_node_number = (n_node_1d * n_node_1d) - 1;
976 cube_local_node_number =
977 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
980 Node* quad_corner_node_pt =
981 quad_el_pt->
node_pt(quad_local_node_number);
984 Node* cube_corner_node_pt =
985 cube_el_pt->node_pt(cube_local_node_number);
988 std::pair<Node*, double> quad_corner_node_and_z_value(
989 quad_corner_node_pt, z_value);
992 quad_corner_node_to_cube_node_map[quad_corner_node_and_z_value] =
998 for (
unsigned i_node = 0; i_node < (n_node_1d * n_node_1d);
1002 if (!has_node_been_setup[i_node])
1005 std::ostringstream error_message_stream;
1008 error_message_stream <<
"There are nodes in element " <<
e
1009 <<
" which have not been constructed!"
1014 OOMPH_CURRENT_FUNCTION,
1015 OOMPH_EXCEPTION_LOCATION);
1024 if (node_count != (1 + (n_node_1d - 1) * (
i + 1)) * n_node_in_quad_mesh)
1027 std::ostringstream error_message_stream;
1030 int expected_n_node =
1031 (1 + (n_node_1d - 1) * (
i + 1)) * n_node_in_quad_mesh;
1034 int node_diff = expected_n_node - node_count;
1037 error_message_stream
1038 <<
"There are meant to be " << expected_n_node
1039 <<
" nodes in the extruded mesh by the\nend of slice " <<
i
1040 <<
" but you have " << node_count <<
". "
1041 <<
"That's a difference of " << node_diff <<
"!";
1045 OOMPH_CURRENT_FUNCTION,
1046 OOMPH_EXCEPTION_LOCATION);
1055 std::set<const Node*> all_quad_boundary_nodes_pt;
1058 int n_quad_mesh_boundary_node = 0;
1061 for (
unsigned b = 0; b < n_boundary_in_quad_mesh; b++)
1067 for (
unsigned n = 0; n < n_boundary_node; n++)
1071 b, n)) == all_quad_boundary_nodes_pt.end())
1074 n_quad_mesh_boundary_node++;
1077 all_quad_boundary_nodes_pt.insert(
1089 int n_expected_boundary_node = 2 * quad_mesh_pt->
nnode();
1093 n_expected_boundary_node +=
1094 (((n_node_1d - 1) * Nz) - 1) * n_quad_mesh_boundary_node;
1097 int n_extruded_mesh_boundary_node = 0;
1100 std::set<const Node*> all_extruded_boundary_nodes_pt;
1103 for (
unsigned b = 0; b < n_boundary_in_extruded_mesh; b++)
1106 unsigned n_boundary_node = nboundary_node(b);
1109 for (
unsigned n = 0; n < n_boundary_node; n++)
1112 if (all_extruded_boundary_nodes_pt.find(this->boundary_node_pt(b, n)) ==
1113 all_extruded_boundary_nodes_pt.end())
1116 n_extruded_mesh_boundary_node++;
1119 all_extruded_boundary_nodes_pt.insert(this->boundary_node_pt(b, n));
1128 if (n_extruded_mesh_boundary_node != n_expected_boundary_node)
1131 std::ostringstream error_message_stream;
1134 int node_diff = n_expected_boundary_node - n_extruded_mesh_boundary_node;
1137 error_message_stream <<
"There should be " << n_expected_boundary_node
1138 <<
" boundary nodes in the extruded mesh but there"
1139 <<
"\nare only " << n_extruded_mesh_boundary_node
1140 <<
" boundary nodes. That's a difference of "
1141 << node_diff <<
"!" << std::endl;
1145 OOMPH_CURRENT_FUNCTION,
1146 OOMPH_EXCEPTION_LOCATION);
1153 oomph_info <<
"Time taken to extrude mesh [sec]: "
1170 std::map<Domain*, ExtrudedDomain*> domain_to_extruded_domain_map;
1173 for (
unsigned e = 0;
e < n_element_in_quad_mesh;
e++)
1194 std::map<Domain*, ExtrudedDomain*>::iterator it;
1197 it = domain_to_extruded_domain_map.find(domain_pt);
1200 if (it != domain_to_extruded_domain_map.end())
1203 extruded_domain_pt = it->second;
1209 extruded_domain_pt =
new ExtrudedDomain(domain_pt, Nz, Zmin, Zmax);
1212 domain_to_extruded_domain_map[domain_pt] = extruded_domain_pt;
1215 Extruded_domain_pt.push_back(extruded_domain_pt);
1219 for (
unsigned i = 0;
i < Nz;
i++)
1224 (
i * n_macro_element +
1232 unsigned element_index =
i * n_element_in_quad_mesh +
e;
1235 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(Element_pt[element_index]);
1238 el_pt->set_macro_elem_pt(extruded_macro_element_pt);
1244 for (
unsigned i = 0;
i < n_dim - 1;
i++)
1257 el_pt->s_macro_ll(n_dim - 1) = -1.0;
1261 el_pt->s_macro_ur(n_dim - 1) = 1.0;
1270 oomph_info <<
"Time taken to set up extruded macro-elements [sec]: "
1280 if (check_for_repeated_nodes())
1284 OOMPH_CURRENT_FUNCTION,
1285 OOMPH_EXCEPTION_LOCATION);
1289 if (node_count != n_node_in_cube_mesh)
1292 std::ostringstream error_message_stream;
1295 error_message_stream <<
"An incorrect number of nodes have been created! "
1296 <<
"There are\nmeant to be " << n_node_in_cube_mesh
1297 <<
" nodes in the extruded mesh but\nonly "
1298 << node_count <<
" have actually been constructed!"
1303 OOMPH_CURRENT_FUNCTION,
1304 OOMPH_EXCEPTION_LOCATION);
1310 setup_boundary_element_info();
1316 oomph_info <<
"Total time for extruded mesh creation/setup [sec]: "
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
unsigned nmacro_element()
Number of macro elements in domain.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void build_mesh(QuadMeshBase *quad_mesh_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Base class for ExtrudedDomains with curvilinear and/or time-dependent boundaries. ExtrudedDomain boun...
ExtrudedMacroElement * macro_element_pt(const unsigned &i)
Access to i-th extruded macro element.
DRAIG: FILL IN COMPLETE DESCRIPTION ONCE FINISHED...
A general Finite Element class.
unsigned nnode() const
Return the number of nodes.
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Domain *& domain_pt()
Access function to the Domain_pt.
unsigned & macro_element_number()
Access function to the Macro_element_number.
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
unsigned nboundary() const
Return number of boundaries.
unsigned long nnode() const
Return number of nodes in the mesh.
unsigned long nelement() const
Return number of elements in the mesh.
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
bool is_hanging() const
Test whether the node is geometrically hanging.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
class oomph::MeshExtrusionHelpers::ExtrusionHelper Mesh_extrusion_helper
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...