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.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
unsigned nnode() const
Return the number of nodes.
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
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.
unsigned nboundary() const
Return number of boundaries.
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
unsigned long nnode() const
Return number of nodes in the mesh.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
unsigned long nelement() const
Return number of elements in the mesh.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
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...
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_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
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...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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...