14 : hoNodes(PETSC_TRUE), defMaxLevel(0), optPrefix(
""), countEle(0),
15 countVertEle(0), nbVertices(0) {}
27 "Wrong parameter -max_post_proc_ref_level "
28 "should be positive number");
38 const int max_level = defMaxLevel;
41 moab::Interface &moab_ref = core_ref;
43 auto create_reference_element = [&moab_ref]() {
45 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
47 for (
int nn = 0; nn < 4; nn++) {
49 moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
52 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
59 auto refine_ref_tetrahedron = [
this, &m_field_ref, max_level]() {
66 for (
int ll = 0; ll != max_level; ++ll) {
71 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
75 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
80 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
87 auto get_ref_gauss_pts_and_shape_functions = [
this, max_level, &moab_ref,
90 for (
int ll = 0; ll != max_level + 1; ++ll) {
93 m_field_ref.
getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
97 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
98 CHKERR moab_ref.add_entities(meshset, tets);
99 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
100 CHKERR moab_ref.delete_entities(&meshset, 1);
103 CHKERR moab_ref.get_connectivity(tets, elem_nodes,
false);
105 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
106 gauss_pts.resize(elem_nodes.size(), 4,
false);
107 std::map<EntityHandle, int> little_map;
108 Range::iterator nit = elem_nodes.begin();
109 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
110 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
111 little_map[*nit] = gg;
113 gauss_pts = trans(gauss_pts);
115 auto &ref_tets = levelRef[ll];
116 Range::iterator tit = tets.begin();
117 for (
int tt = 0; tit != tets.end(); ++tit, ++tt) {
120 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
122 ref_tets.resize(tets.size(), num_nodes);
124 for (
int nn = 0; nn != num_nodes; ++nn) {
125 ref_tets(tt, nn) = little_map[conn[nn]];
129 auto &shape_functions = levelShapeFunctions[ll];
130 shape_functions.resize(elem_nodes.size(), 4);
132 &gauss_pts(1, 0), &gauss_pts(2, 0), elem_nodes.size());
137 levelRef.resize(max_level + 1);
138 levelGaussPtsOnRefMesh.resize(max_level + 1);
139 levelShapeFunctions.resize(max_level + 1);
141 CHKERR create_reference_element();
142 CHKERR refine_ref_tetrahedron();
143 CHKERR get_ref_gauss_pts_and_shape_functions();
155 <<
"Refinement for hexes is not implemented";
159 moab::Interface &moab_ref = core_ref;
161 auto create_reference_element = [&moab_ref]() {
163 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
164 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
166 for (
int nn = 0; nn < 8; nn++)
167 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
169 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
173 auto add_ho_nodes = [&]() {
176 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
178 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
179 CHKERR moab_ref.add_entities(meshset, hexes);
180 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
181 CHKERR moab_ref.delete_entities(&meshset, 1);
185 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
188 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
190 CHKERR moab_ref.get_connectivity(hexes, hexes_nodes,
false);
191 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
192 gauss_pts.resize(hexes_nodes.size(), 4,
false);
194 for (
auto node : hexes_nodes) {
195 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
196 little_map[node] = gg;
199 gauss_pts = trans(gauss_pts);
203 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
206 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
208 auto &ref_hexes = levelRef[0];
209 for (
auto hex : hexes) {
212 CHKERR moab_ref.get_connectivity(hex, conn, num_nodes,
false);
213 if (ref_hexes.size2() != num_nodes) {
214 ref_hexes.resize(hexes.size(), num_nodes);
216 for (
int nn = 0; nn != num_nodes; ++nn) {
217 ref_hexes(hh, nn) = little_map[conn[nn]];
224 auto set_shape_functions = [&]() {
226 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
227 auto &shape_functions = levelShapeFunctions[0];
228 const auto nb_gauss_pts = gauss_pts.size2();
229 shape_functions.resize(nb_gauss_pts, 8);
230 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
231 const double ksi = gauss_pts(0, gg);
232 const double zeta = gauss_pts(1, gg);
233 const double eta = gauss_pts(2, gg);
247 levelGaussPtsOnRefMesh.resize(1);
248 levelShapeFunctions.resize(1);
250 CHKERR create_reference_element();
253 std::map<EntityHandle, int> little_map;
254 CHKERR set_gauss_pts(little_map);
255 CHKERR set_ref_hexes(little_map);
256 CHKERR set_shape_functions();
265 const int max_level = defMaxLevel;
268 moab::Interface &moab_ref = core_ref;
270 auto create_reference_element = [&moab_ref]() {
272 constexpr double base_coords[] = {
285 for (
int nn = 0; nn != 3; ++nn)
286 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
288 CHKERR moab_ref.create_element(MBTRI, nodes, 3, tri);
291 CHKERR moab_ref.get_adjacencies(&tri, 1, 1,
true, edges,
292 moab::Interface::UNION);
297 CHKERR create_reference_element();
302 auto refine_ref_triangles = [
this, &m_field_ref, max_level]() {
310 for (
int ll = 0; ll != max_level; ++ll) {
315 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
319 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
323 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
348 auto get_ref_gauss_pts_and_shape_functions = [
this, max_level, &moab_ref,
351 for (
int ll = 0; ll != max_level + 1; ++ll) {
360 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
361 CHKERR moab_ref.add_entities(meshset, tris);
362 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
363 CHKERR moab_ref.delete_entities(&meshset, 1);
367 CHKERR moab_ref.get_connectivity(tris, elem_nodes,
false);
369 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
370 gauss_pts.resize(elem_nodes.size(), 3,
false);
371 std::map<EntityHandle, int> little_map;
372 Range::iterator nit = elem_nodes.begin();
373 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
374 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
375 little_map[*nit] = gg;
377 gauss_pts = trans(gauss_pts);
379 auto &ref_tris = levelRef[ll];
380 Range::iterator tit = tris.begin();
381 for (
int tt = 0; tit != tris.end(); ++tit, ++tt) {
384 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
386 ref_tris.resize(tris.size(), num_nodes);
388 for (
int nn = 0; nn != num_nodes; ++nn) {
389 ref_tris(tt, nn) = little_map[conn[nn]];
393 auto &shape_functions = levelShapeFunctions[ll];
394 shape_functions.resize(elem_nodes.size(), 3);
396 &gauss_pts(1, 0), elem_nodes.size());
401 levelRef.resize(max_level + 1);
402 levelGaussPtsOnRefMesh.resize(max_level + 1);
403 levelShapeFunctions.resize(max_level + 1);
405 CHKERR refine_ref_triangles();
406 CHKERR get_ref_gauss_pts_and_shape_functions();
418 <<
"Refinement for quad is not implemented";
422 moab::Interface &moab_ref = core_ref;
424 auto create_reference_element = [&moab_ref]() {
426 constexpr double base_coords[] = {
442 for (
int nn = 0; nn < 4; nn++)
443 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
445 CHKERR moab_ref.create_element(MBQUAD, nodes, 4, quad);
449 auto add_ho_nodes = [&]() {
452 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
454 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
455 CHKERR moab_ref.add_entities(meshset, quads);
456 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
457 CHKERR moab_ref.delete_entities(&meshset, 1);
461 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
464 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
466 CHKERR moab_ref.get_connectivity(quads, quads_nodes,
false);
467 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
468 gauss_pts.resize(quads_nodes.size(), 4,
false);
470 for (
auto node : quads_nodes) {
471 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
472 little_map[node] = gg;
475 gauss_pts = trans(gauss_pts);
479 auto set_ref_quads = [&](std::map<EntityHandle, int> &little_map) {
482 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
484 auto &ref_quads = levelRef[0];
485 for (
auto quad : quads) {
488 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes,
false);
489 if (ref_quads.size2() != num_nodes) {
490 ref_quads.resize(quads.size(), num_nodes);
492 for (
int nn = 0; nn != num_nodes; ++nn) {
493 ref_quads(hh, nn) = little_map[conn[nn]];
500 auto set_shape_functions = [&]() {
502 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
503 auto &shape_functions = levelShapeFunctions[0];
504 const auto nb_gauss_pts = gauss_pts.size2();
505 shape_functions.resize(nb_gauss_pts, 4);
506 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
507 const double ksi = gauss_pts(0, gg);
508 const double zeta = gauss_pts(1, gg);
518 levelGaussPtsOnRefMesh.resize(1);
519 levelShapeFunctions.resize(1);
521 CHKERR create_reference_element();
524 std::map<EntityHandle, int> little_map;
525 CHKERR set_gauss_pts(little_map);
526 CHKERR set_ref_quads(little_map);
527 CHKERR set_shape_functions();
539 <<
"Refinement for edges is not implemented";
542 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
549 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
550 gauss_pts.resize(2, nb_nodes,
false);
554 for (; nn != 2; ++nn) {
555 gauss_pts(0, nn) =
static_cast<double>(nn);
560 gauss_pts(0, nn) = 0.5;
567 auto set_ref_edges = [&](std::map<EntityHandle, int> &little_map) {
571 int nb_edges = level + 1;
577 auto &ref_edges = levelRef[level];
578 ref_edges.resize(nb_edges, nb_nodes,
false);
580 for (
int ee = 0; ee != nb_edges; ++ee) {
582 for (; nn != 2; ++nn) {
583 ref_edges(ee, nn) = nb_nodes * ee + nn;
586 ref_edges(ee, nn) = nb_nodes * ee + 2;
593 auto set_shape_functions = [&]() {
595 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
596 auto &shape_functions = levelShapeFunctions[0];
597 const auto nb_gauss_pts = gauss_pts.size2();
598 shape_functions.resize(nb_gauss_pts, 2);
599 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
600 const double ksi = gauss_pts(0, gg);
608 levelGaussPtsOnRefMesh.resize(1);
609 levelShapeFunctions.resize(1);
611 std::map<EntityHandle, int> little_map;
612 CHKERR set_gauss_pts(little_map);
613 CHKERR set_ref_edges(little_map);
614 CHKERR set_shape_functions();
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Deprecated interface functions.
Mesh refinement interface.
PostProcGenerateRefMeshBase()
MoFEMErrorCode getOptions()
Element for postprocessing. Uses MoAB to generate post-processing mesh.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.