14 : hoNodes(PETSC_TRUE), defMaxLevel(0), countEle(0), countVertEle(0),
24 "-max_post_ho_nodes", &
hoNodes, PETSC_NULL);
28 "Wrong parameter -max_post_proc_ref_level "
29 "should be positive number");
37 const int max_level = defMaxLevel;
40 moab::Interface &moab_ref = core_ref;
42 auto create_reference_element = [&moab_ref]() {
44 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
46 for (
int nn = 0; nn < 4; nn++) {
48 moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
51 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
58 auto refine_ref_tetrahedron = [
this, &m_field_ref, max_level]() {
63 m_field_ref.
getInterface<BitRefManager>()->setBitRefLevelByDim(
65 for (
int ll = 0; ll != max_level; ++ll) {
70 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
74 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
77 MeshRefinement *m_ref;
79 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
86 auto get_ref_gauss_pts_and_shape_functions = [
this, max_level, &moab_ref,
89 for (
int ll = 0; ll != max_level + 1; ++ll) {
92 m_field_ref.
getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
96 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
97 CHKERR moab_ref.add_entities(meshset, tets);
98 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
99 CHKERR moab_ref.delete_entities(&meshset, 1);
102 CHKERR moab_ref.get_connectivity(tets, elem_nodes,
false);
104 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
105 gauss_pts.resize(elem_nodes.size(), 4,
false);
106 std::map<EntityHandle, int> little_map;
107 Range::iterator nit = elem_nodes.begin();
108 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
109 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
110 little_map[*nit] = gg;
112 gauss_pts = trans(gauss_pts);
114 auto &ref_tets = levelRef[ll];
115 Range::iterator tit = tets.begin();
116 for (
int tt = 0; tit != tets.end(); ++tit, ++tt) {
119 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
121 ref_tets.resize(tets.size(), num_nodes);
123 for (
int nn = 0; nn != num_nodes; ++nn) {
124 ref_tets(tt, nn) = little_map[conn[nn]];
128 auto &shape_functions = levelShapeFunctions[ll];
129 shape_functions.resize(elem_nodes.size(), 4);
131 &gauss_pts(1, 0), &gauss_pts(2, 0), elem_nodes.size());
136 levelRef.resize(max_level + 1);
137 levelGaussPtsOnRefMesh.resize(max_level + 1);
138 levelShapeFunctions.resize(max_level + 1);
140 CHKERR create_reference_element();
141 CHKERR refine_ref_tetrahedron();
142 CHKERR get_ref_gauss_pts_and_shape_functions();
153 <<
"Refinement for hexes is not implemented";
157 moab::Interface &moab_ref = core_ref;
159 auto create_reference_element = [&moab_ref]() {
161 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
162 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
164 for (
int nn = 0; nn < 8; nn++)
165 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
167 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
171 auto add_ho_nodes = [&]() {
174 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
176 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
177 CHKERR moab_ref.add_entities(meshset, hexes);
178 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
179 CHKERR moab_ref.delete_entities(&meshset, 1);
183 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
186 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
188 CHKERR moab_ref.get_connectivity(hexes, hexes_nodes,
false);
189 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
190 gauss_pts.resize(hexes_nodes.size(), 4,
false);
192 for (
auto node : hexes_nodes) {
193 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
194 little_map[node] = gg;
197 gauss_pts = trans(gauss_pts);
201 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
204 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
206 auto &ref_hexes = levelRef[0];
207 for (
auto hex : hexes) {
210 CHKERR moab_ref.get_connectivity(hex, conn, num_nodes,
false);
211 if (ref_hexes.size2() != num_nodes) {
212 ref_hexes.resize(hexes.size(), num_nodes);
214 for (
int nn = 0; nn != num_nodes; ++nn) {
215 ref_hexes(hh, nn) = little_map[conn[nn]];
222 auto set_shape_functions = [&]() {
224 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
225 auto &shape_functions = levelShapeFunctions[0];
226 const auto nb_gauss_pts = gauss_pts.size2();
227 shape_functions.resize(nb_gauss_pts, 8);
228 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
229 const double ksi = gauss_pts(0, gg);
230 const double zeta = gauss_pts(1, gg);
231 const double eta = gauss_pts(2, gg);
245 levelGaussPtsOnRefMesh.resize(1);
246 levelShapeFunctions.resize(1);
248 CHKERR create_reference_element();
251 std::map<EntityHandle, int> little_map;
252 CHKERR set_gauss_pts(little_map);
253 CHKERR set_ref_hexes(little_map);
254 CHKERR set_shape_functions();
262 const int max_level = defMaxLevel;
265 moab::Interface &moab_ref = core_ref;
267 auto create_reference_element = [&moab_ref]() {
269 constexpr double base_coords[] = {
282 for (
int nn = 0; nn != 3; ++nn)
283 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
285 CHKERR moab_ref.create_element(MBTRI, nodes, 3, tri);
288 CHKERR moab_ref.get_adjacencies(&tri, 1, 1,
true, edges,
289 moab::Interface::UNION);
294 CHKERR create_reference_element();
299 auto refine_ref_triangles = [
this, &m_field_ref, max_level]() {
307 for (
int ll = 0; ll != max_level; ++ll) {
312 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
316 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
320 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
345 auto get_ref_gauss_pts_and_shape_functions = [
this, max_level, &moab_ref,
348 for (
int ll = 0; ll != max_level + 1; ++ll) {
357 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
358 CHKERR moab_ref.add_entities(meshset, tris);
359 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
360 CHKERR moab_ref.delete_entities(&meshset, 1);
364 CHKERR moab_ref.get_connectivity(tris, elem_nodes,
false);
366 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
367 gauss_pts.resize(elem_nodes.size(), 3,
false);
368 std::map<EntityHandle, int> little_map;
369 Range::iterator nit = elem_nodes.begin();
370 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
371 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
372 little_map[*nit] = gg;
374 gauss_pts = trans(gauss_pts);
376 auto &ref_tris = levelRef[ll];
377 Range::iterator tit = tris.begin();
378 for (
int tt = 0; tit != tris.end(); ++tit, ++tt) {
381 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
383 ref_tris.resize(tris.size(), num_nodes);
385 for (
int nn = 0; nn != num_nodes; ++nn) {
386 ref_tris(tt, nn) = little_map[conn[nn]];
390 auto &shape_functions = levelShapeFunctions[ll];
391 shape_functions.resize(elem_nodes.size(), 3);
393 &gauss_pts(1, 0), elem_nodes.size());
398 levelRef.resize(max_level + 1);
399 levelGaussPtsOnRefMesh.resize(max_level + 1);
400 levelShapeFunctions.resize(max_level + 1);
402 CHKERR refine_ref_triangles();
403 CHKERR get_ref_gauss_pts_and_shape_functions();
414 <<
"Refinement for quad is not implemented";
418 moab::Interface &moab_ref = core_ref;
420 auto create_reference_element = [&moab_ref]() {
422 constexpr double base_coords[] = {
438 for (
int nn = 0; nn < 4; nn++)
439 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
441 CHKERR moab_ref.create_element(MBQUAD, nodes, 4, quad);
445 auto add_ho_nodes = [&]() {
448 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
450 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
451 CHKERR moab_ref.add_entities(meshset, quads);
452 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
453 CHKERR moab_ref.delete_entities(&meshset, 1);
457 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
460 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
462 CHKERR moab_ref.get_connectivity(quads, quads_nodes,
false);
463 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
464 gauss_pts.resize(quads_nodes.size(), 4,
false);
466 for (
auto node : quads_nodes) {
467 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
468 little_map[node] = gg;
471 gauss_pts = trans(gauss_pts);
475 auto set_ref_quads = [&](std::map<EntityHandle, int> &little_map) {
478 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
480 auto &ref_quads = levelRef[0];
481 for (
auto quad : quads) {
484 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes,
false);
485 if (ref_quads.size2() != num_nodes) {
486 ref_quads.resize(quads.size(), num_nodes);
488 for (
int nn = 0; nn != num_nodes; ++nn) {
489 ref_quads(hh, nn) = little_map[conn[nn]];
496 auto set_shape_functions = [&]() {
498 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
499 auto &shape_functions = levelShapeFunctions[0];
500 const auto nb_gauss_pts = gauss_pts.size2();
501 shape_functions.resize(nb_gauss_pts, 4);
502 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
503 const double ksi = gauss_pts(0, gg);
504 const double zeta = gauss_pts(1, gg);
514 levelGaussPtsOnRefMesh.resize(1);
515 levelShapeFunctions.resize(1);
517 CHKERR create_reference_element();
520 std::map<EntityHandle, int> little_map;
521 CHKERR set_gauss_pts(little_map);
522 CHKERR set_ref_quads(little_map);
523 CHKERR set_shape_functions();
534 <<
"Refinement for edges is not implemented";
537 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
544 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
545 gauss_pts.resize(2, nb_nodes,
false);
549 for (; nn != 2; ++nn) {
550 gauss_pts(0, nn) =
static_cast<double>(nn);
555 gauss_pts(0, nn) = 0.5;
562 auto set_ref_edges = [&](std::map<EntityHandle, int> &little_map) {
566 int nb_edges = level + 1;
572 auto &ref_edges = levelRef[level];
573 ref_edges.resize(nb_edges, nb_nodes,
false);
575 for (
int ee = 0; ee != nb_edges; ++ee) {
577 for (; nn != 2; ++nn) {
578 ref_edges(ee, nn) = nb_nodes * ee + nn;
581 ref_edges(ee, nn) = nb_nodes * ee + 2;
588 auto set_shape_functions = [&]() {
590 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
591 auto &shape_functions = levelShapeFunctions[0];
592 const auto nb_gauss_pts = gauss_pts.size2();
593 shape_functions.resize(nb_gauss_pts, 2);
594 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
595 const double ksi = gauss_pts(0, gg);
603 levelGaussPtsOnRefMesh.resize(1);
604 levelShapeFunctions.resize(1);
606 std::map<EntityHandle, int> little_map;
607 CHKERR set_gauss_pts(little_map);
608 CHKERR set_ref_edges(little_map);
609 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)
double zeta
Viscous hardening.
Deprecated interface functions.
Mesh refinement interface.
virtual MoFEMErrorCode getOptions(std::string prefix)
PostProcGenerateRefMeshBase()
Element for postprocessing. Uses MoAB to generate post-processing mesh.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.