43 template<
class ELEMENT>
44 Vector<std::pair<unsigned, int>> ExtrudedCubeMeshFromQuadMesh<
45 ELEMENT>::get_element_boundary_information(QuadMeshBase* quad_mesh_pt,
46 FiniteElement* quad_el_pt)
49 Vector<std::pair<unsigned, int>> boundary_info;
52 unsigned n_boundary = quad_mesh_pt->nboundary();
55 for (
unsigned b = 0; b < n_boundary; b++)
58 unsigned n_boundary_element = quad_mesh_pt->nboundary_element(b);
61 for (
unsigned e = 0; e < n_boundary_element; e++)
64 if (quad_el_pt == quad_mesh_pt->boundary_element_pt(b, e))
67 std::pair<unsigned, int> boundary_and_face_index;
70 boundary_and_face_index.first = b;
73 boundary_and_face_index.second =
74 quad_mesh_pt->face_index_at_boundary(b, e);
77 boundary_info.push_back(boundary_and_face_index);
98 template<
class ELEMENT>
100 QuadMeshBase* quad_mesh_pt, TimeStepper* time_stepper_pt)
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);
122 double t_mesh_setup_start = TimingHelpers::timer();
125 double t_mesh_construct_start = TimingHelpers::timer();
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();
142 Vector<Vector<std::pair<unsigned, int>>> boundary_information(
143 n_element_in_quad_mesh);
146 for (
unsigned e = 0; e < n_element_in_quad_mesh; e++)
149 FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(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;
175 if ((n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode_1d()) ||
176 (n_node_1d * n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode()))
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++)
199 if (quad_mesh_pt->node_pt(i)->is_hanging())
202 std::ostringstream warning_message_stream;
205 warning_message_stream <<
"Extrusion machinery is not currently able "
206 <<
"to deal with hanging nodes!" << std::endl;
209 throw OomphLibError(warning_message_stream.str(),
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++)
251 FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(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;
269 Vector<std::pair<unsigned, int>> el_boundary_information =
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;
304 FiniteElement* copy_cube_el_pt =
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);
357 Vector<std::pair<Edge, double>> edges_and_z_value;
361 std::vector<bool> has_edge_been_done(n_edge,
false);
366 int last_edge = int(QuadTreeNames::N + n_edge);
369 for (
int i_edge = QuadTreeNames::N; i_edge < last_edge; i_edge++)
380 case QuadTreeNames::N:
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));
389 case QuadTreeNames::E:
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);
398 case QuadTreeNames::S:
400 node1_pt = quad_el_pt->node_pt(0);
403 node2_pt = quad_el_pt->node_pt(n_node_1d - 1);
407 case QuadTreeNames::W:
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"
423 <<
"be either N/E/S/W, i.e. " << QuadTreeNames::N <<
","
424 << QuadTreeNames::E <<
"," << QuadTreeNames::S <<
" or "
425 << QuadTreeNames::W << std::endl;
428 throw OomphLibError(error_message_stream.str(),
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())
451 Vector<Node*> edge_nodes_pt = it->second;
455 if (edge_nodes_pt.size() != n_node_1d)
458 throw OomphLibError(
"Not enough nodes in edge_nodes_pt!",
459 OOMPH_CURRENT_FUNCTION,
460 OOMPH_EXCEPTION_LOCATION);
467 case QuadTreeNames::N:
469 has_edge_been_done[0] =
true;
473 case QuadTreeNames::E:
475 has_edge_been_done[1] =
true;
479 case QuadTreeNames::S:
481 has_edge_been_done[2] =
true;
485 case QuadTreeNames::W:
487 has_edge_been_done[3] =
true;
493 throw OomphLibError(
"Incorrect edge type!",
494 OOMPH_CURRENT_FUNCTION,
495 OOMPH_EXCEPTION_LOCATION);
506 for (
unsigned i_node = 0; i_node < n_node_1d; i_node++)
512 case QuadTreeNames::N:
514 quad_local_node_number =
515 ((n_node_1d * n_node_1d) - 1) - i_node;
519 case QuadTreeNames::E:
521 quad_local_node_number =
522 (n_node_1d - 1) + (i_node * n_node_1d);
526 case QuadTreeNames::S:
528 quad_local_node_number = i_node;
532 case QuadTreeNames::W:
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);
695 if (quad_node_pt->is_on_boundary())
698 std::set<unsigned>* boundaries_pt;
701 quad_node_pt->get_boundaries_pt(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);
814 if (quad_node_pt->is_on_boundary())
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;
832 quad_node_pt->get_boundaries_pt(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])
886 Vector<Node*> edge_nodes_pt(n_node_1d, 0);
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!"
1013 throw OomphLibError(error_message_stream.str(),
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 <<
"!";
1044 throw OomphLibError(error_message_stream.str(),
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++)
1064 unsigned n_boundary_node = quad_mesh_pt->nboundary_node(b);
1067 for (
unsigned n = 0; n < n_boundary_node; n++)
1070 if (all_quad_boundary_nodes_pt.find(quad_mesh_pt->boundary_node_pt(
1071 b, n)) == all_quad_boundary_nodes_pt.end())
1074 n_quad_mesh_boundary_node++;
1077 all_quad_boundary_nodes_pt.insert(
1078 quad_mesh_pt->boundary_node_pt(b, n));
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;
1144 throw OomphLibError(error_message_stream.str(),
1145 OOMPH_CURRENT_FUNCTION,
1146 OOMPH_EXCEPTION_LOCATION);
1153 oomph_info <<
"Time taken to extrude mesh [sec]: "
1154 << TimingHelpers::timer() - t_mesh_construct_start
1159 double t_mesh_macro_start = TimingHelpers::timer();
1170 std::map<Domain*, ExtrudedDomain*> domain_to_extruded_domain_map;
1173 for (
unsigned e = 0; e < n_element_in_quad_mesh; e++)
1178 QElementBase* quad_el_pt =
1179 dynamic_cast<QElementBase*
>(quad_mesh_pt->finite_element_pt(e));
1182 if (quad_el_pt->macro_elem_pt() != 0)
1185 Domain* domain_pt = quad_el_pt->macro_elem_pt()->domain_pt();
1188 unsigned n_macro_element = domain_pt->nmacro_element();
1191 ExtrudedDomain* extruded_domain_pt = 0;
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 +
1225 quad_el_pt->macro_elem_pt()->macro_element_number());
1228 ExtrudedMacroElement* extruded_macro_element_pt =
1229 extruded_domain_pt->macro_element_pt(macro_id);
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++)
1248 el_pt->s_macro_ll(i) = quad_el_pt->s_macro_ll(i);
1252 el_pt->s_macro_ur(i) = quad_el_pt->s_macro_ur(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]: "
1271 << TimingHelpers::timer() - t_mesh_macro_start << std::endl;
1280 if (check_for_repeated_nodes())
1283 OomphLibWarning(
"The mesh contains 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!"
1302 OomphLibWarning(error_message_stream.str(),
1303 OOMPH_CURRENT_FUNCTION,
1304 OOMPH_EXCEPTION_LOCATION);
1310 setup_boundary_element_info();
1316 oomph_info <<
"Total time for extruded mesh creation/setup [sec]: "
1317 << TimingHelpers::timer() - t_mesh_setup_start << std::endl;
void build_mesh(QuadMeshBase *quad_mesh_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
class oomph::MeshExtrusionHelpers::ExtrusionHelper Mesh_extrusion_helper
////////////////////////////////////////////////////////////////////// //////////////////////////////...