v0.16.0
Loading...
Searching...
No Matches
PostProcBrokenMeshInMoabBase.cpp
Go to the documentation of this file.
1/**
2 * @file PostProc.cpp
3 * @brief Post processing elements and operators
4 *
5 * @copyright Copyright (c) 2022
6 *
7 */
8
9#include <MoFEM.hpp>
10
11namespace MoFEM {
12
13PostProcGenerateRefMeshBase::PostProcGenerateRefMeshBase()
14 : hoNodes(PETSC_TRUE), defMaxLevel(0), countEle(0), countVertEle(0),
15 nbVertices(0) {}
16
17MoFEMErrorCode PostProcGenerateRefMeshBase::getOptions(std::string prefix) {
18 MoFEMFunctionBegin;
19
20 std::string opt1 = prefix.size() ? "-" + prefix + "_max_post_proc_ref_level"
21 : "-max_post_proc_ref_level";
22 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, opt1.c_str(), &defMaxLevel, PETSC_NULLPTR);
23
24 std::string opt2 = prefix.size() ? "-" + prefix + "_max_post_ho_nodes"
25 : "-max_post_ho_nodes";
26 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, opt2.c_str(), &hoNodes, PETSC_NULLPTR);
27
28 if (defMaxLevel < 0)
29 SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
30 "Wrong parameter -max_post_proc_ref_level "
31 "should be positive number");
32
33 MoFEMFunctionReturn(0);
34};
35
36MoFEMErrorCode PostProcGenerateRefMesh<MBTET>::generateReferenceElementMesh() {
37 MoFEMFunctionBegin;
38
39 const int max_level = defMaxLevel;
40
41 moab::Core core_ref;
42 moab::Interface &moab_ref = core_ref;
43
44 auto create_reference_element = [&moab_ref]() {
45 MoFEMFunctionBegin;
46 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
47 EntityHandle nodes[4];
48 for (int nn = 0; nn < 4; nn++) {
49 CHKERR
50 moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
51 }
52 EntityHandle tet;
53 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
54 MoFEMFunctionReturn(0);
55 };
56
57 MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
58 MoFEM::Interface &m_field_ref = m_core_ref;
59
60 auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
61 MoFEMFunctionBegin;
62 // seed ref mofem database by setting bit ref level to reference
63 // tetrahedron
64 CHKERR
65 m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
66 0, 3, BitRefLevel().set(0));
67 for (int ll = 0; ll != max_level; ++ll) {
68 MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "PostProc", "Refine Level %d",
69 ll);
70 Range edges;
71 CHKERR m_field_ref.getInterface<BitRefManager>()
72 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
73 BitRefLevel().set(), MBEDGE, edges);
74 Range tets;
75 CHKERR m_field_ref.getInterface<BitRefManager>()
76 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
77 BitRefLevel(ll).set(), MBTET, tets);
78 // refine mesh
79 MeshRefinement *m_ref;
80 CHKERR m_field_ref.getInterface(m_ref);
81 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
82 BitRefLevel().set(ll + 1));
83 CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
84 }
85 MoFEMFunctionReturn(0);
86 };
87
88 auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
89 &m_field_ref]() {
90 MoFEMFunctionBegin;
91 for (int ll = 0; ll != max_level + 1; ++ll) {
92 Range tets;
93 CHKERR
94 m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
95 BitRefLevel().set(ll), BitRefLevel().set(ll), MBTET, tets);
96 if (hoNodes) {
97 EntityHandle meshset;
98 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
99 CHKERR moab_ref.add_entities(meshset, tets);
100 CHKERR moab_ref.convert_entities(meshset, true, false, false);
101 CHKERR moab_ref.delete_entities(&meshset, 1);
102 }
103 Range elem_nodes;
104 CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
105
106 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
107 gauss_pts.resize(elem_nodes.size(), 4, false);
108 std::map<EntityHandle, int> little_map;
109 Range::iterator nit = elem_nodes.begin();
110 for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
111 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
112 little_map[*nit] = gg;
113 }
114 gauss_pts = trans(gauss_pts);
115
116 auto &ref_tets = levelRef[ll];
117 Range::iterator tit = tets.begin();
118 for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
119 const EntityHandle *conn;
120 int num_nodes;
121 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
122 if (tt == 0) {
123 ref_tets.resize(tets.size(), num_nodes);
124 }
125 for (int nn = 0; nn != num_nodes; ++nn) {
126 ref_tets(tt, nn) = little_map[conn[nn]];
127 }
128 }
129
130 auto &shape_functions = levelShapeFunctions[ll];
131 shape_functions.resize(elem_nodes.size(), 4);
132 CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
133 &gauss_pts(1, 0), &gauss_pts(2, 0), elem_nodes.size());
134 }
135 MoFEMFunctionReturn(0);
136 };
137
138 levelRef.resize(max_level + 1);
139 levelGaussPtsOnRefMesh.resize(max_level + 1);
140 levelShapeFunctions.resize(max_level + 1);
141
142 CHKERR create_reference_element();
143 CHKERR refine_ref_tetrahedron();
144 CHKERR get_ref_gauss_pts_and_shape_functions();
145
146 MoFEMFunctionReturn(0);
147}
148
149MoFEMErrorCode PostProcGenerateRefMesh<MBHEX>::generateReferenceElementMesh() {
150 MoFEMFunctionBegin;
151
152#ifndef NDEBUG
153 if (defMaxLevel > 0)
154 MOFEM_LOG("WORLD", Sev::warning)
155 << "Refinement for hexes is not implemented";
156#endif
157
158 moab::Core core_ref;
159 moab::Interface &moab_ref = core_ref;
160
161 auto create_reference_element = [&moab_ref]() {
162 MoFEMFunctionBegin;
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};
165 EntityHandle nodes[8];
166 for (int nn = 0; nn < 8; nn++)
167 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
168 EntityHandle hex;
169 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
170 MoFEMFunctionReturn(0);
171 };
172
173 auto add_ho_nodes = [&]() {
174 MoFEMFunctionBegin;
175 Range hexes;
176 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
177 EntityHandle meshset;
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);
182 MoFEMFunctionReturn(0);
183 };
184
185 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
186 MoFEMFunctionBegin;
187 Range hexes;
188 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
189 Range hexes_nodes;
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);
193 size_t gg = 0;
194 for (auto node : hexes_nodes) {
195 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
196 little_map[node] = gg;
197 ++gg;
198 }
199 gauss_pts = trans(gauss_pts);
200 MoFEMFunctionReturn(0);
201 };
202
203 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
204 MoFEMFunctionBegin;
205 Range hexes;
206 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
207 size_t hh = 0;
208 auto &ref_hexes = levelRef[0];
209 for (auto hex : hexes) {
210 const EntityHandle *conn;
211 int num_nodes;
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);
215 }
216 for (int nn = 0; nn != num_nodes; ++nn) {
217 ref_hexes(hh, nn) = little_map[conn[nn]];
218 }
219 ++hh;
220 }
221 MoFEMFunctionReturn(0);
222 };
223
224 auto set_shape_functions = [&]() {
225 MoFEMFunctionBegin;
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);
234 shape_functions(gg, 0) = N_MBHEX0(ksi, zeta, eta);
235 shape_functions(gg, 1) = N_MBHEX1(ksi, zeta, eta);
236 shape_functions(gg, 2) = N_MBHEX2(ksi, zeta, eta);
237 shape_functions(gg, 3) = N_MBHEX3(ksi, zeta, eta);
238 shape_functions(gg, 4) = N_MBHEX4(ksi, zeta, eta);
239 shape_functions(gg, 5) = N_MBHEX5(ksi, zeta, eta);
240 shape_functions(gg, 6) = N_MBHEX6(ksi, zeta, eta);
241 shape_functions(gg, 7) = N_MBHEX7(ksi, zeta, eta);
242 }
243 MoFEMFunctionReturn(0);
244 };
245
246 levelRef.resize(1);
247 levelGaussPtsOnRefMesh.resize(1);
248 levelShapeFunctions.resize(1);
249
250 CHKERR create_reference_element();
251 if (hoNodes)
252 CHKERR add_ho_nodes();
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();
257
258 MoFEMFunctionReturn(0);
259}
260
261MoFEMErrorCode
262PostProcGenerateRefMesh<MBPRISM>::generateReferenceElementMesh() {
263 MoFEMFunctionBegin;
264
265#ifndef NDEBUG
266 if (defMaxLevel > 0)
267 MOFEM_LOG("WORLD", Sev::warning)
268 << "Refinement for prisms is not implemented";
269#endif
270
271 moab::Core core_ref;
272 moab::Interface &moab_ref = core_ref;
273
274 auto create_reference_element = [&moab_ref]() {
275 MoFEMFunctionBegin;
276 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0,
277 0, 0, 1, 1, 0, 1, 0, 1, 1};
278 EntityHandle nodes[6];
279 for (int nn = 0; nn != 6; ++nn)
280 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
281 EntityHandle prism;
282 CHKERR moab_ref.create_element(MBPRISM, nodes, 6, prism);
283 MoFEMFunctionReturn(0);
284 };
285
286 auto add_ho_nodes = [&]() {
287 MoFEMFunctionBegin;
288 Range prisms;
289 CHKERR moab_ref.get_entities_by_type(0, MBPRISM, prisms, true);
290 EntityHandle meshset;
291 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
292 CHKERR moab_ref.add_entities(meshset, prisms);
293 CHKERR moab_ref.convert_entities(meshset, true, true, true);
294 CHKERR moab_ref.delete_entities(&meshset, 1);
295 MoFEMFunctionReturn(0);
296 };
297
298 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
299 MoFEMFunctionBegin;
300 Range prisms;
301 CHKERR moab_ref.get_entities_by_type(0, MBPRISM, prisms, true);
302 Range prism_nodes;
303 CHKERR moab_ref.get_connectivity(prisms, prism_nodes, false);
304 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
305 gauss_pts.resize(prism_nodes.size(), 4, false);
306 size_t gg = 0;
307 for (auto node : prism_nodes) {
308 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
309 little_map[node] = gg;
310 ++gg;
311 }
312 gauss_pts = trans(gauss_pts);
313 MoFEMFunctionReturn(0);
314 };
315
316 auto set_ref_prisms = [&](std::map<EntityHandle, int> &little_map) {
317 MoFEMFunctionBegin;
318 Range prisms;
319 CHKERR moab_ref.get_entities_by_type(0, MBPRISM, prisms, true);
320 size_t pp = 0;
321 auto &ref_prisms = levelRef[0];
322 for (auto prism : prisms) {
323 const EntityHandle *conn;
324 int num_nodes;
325 CHKERR moab_ref.get_connectivity(prism, conn, num_nodes, false);
326 if (ref_prisms.size2() != num_nodes)
327 ref_prisms.resize(prisms.size(), num_nodes);
328 for (int nn = 0; nn != num_nodes; ++nn)
329 ref_prisms(pp, nn) = little_map[conn[nn]];
330 ++pp;
331 }
332 MoFEMFunctionReturn(0);
333 };
334
335 auto set_shape_functions = [&]() {
336 MoFEMFunctionBegin;
337 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
338 auto &shape_functions = levelShapeFunctions[0];
339 const auto nb_gauss_pts = gauss_pts.size2();
340 shape_functions.resize(nb_gauss_pts, 6);
341 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
342 const double ksi = gauss_pts(0, gg);
343 const double eta = gauss_pts(1, gg);
344 const double zeta = gauss_pts(2, gg);
345 const double n0 = 1 - ksi - eta;
346 shape_functions(gg, 0) = n0 * (1 - zeta);
347 shape_functions(gg, 1) = ksi * (1 - zeta);
348 shape_functions(gg, 2) = eta * (1 - zeta);
349 shape_functions(gg, 3) = n0 * zeta;
350 shape_functions(gg, 4) = ksi * zeta;
351 shape_functions(gg, 5) = eta * zeta;
352 }
353 MoFEMFunctionReturn(0);
354 };
355
356 levelRef.resize(1);
357 levelGaussPtsOnRefMesh.resize(1);
358 levelShapeFunctions.resize(1);
359
360 CHKERR create_reference_element();
361 if (hoNodes)
362 CHKERR add_ho_nodes();
363 std::map<EntityHandle, int> little_map;
364 CHKERR set_gauss_pts(little_map);
365 CHKERR set_ref_prisms(little_map);
366 CHKERR set_shape_functions();
367
368 MoFEMFunctionReturn(0);
369}
370
371MoFEMErrorCode PostProcGenerateRefMesh<MBTRI>::generateReferenceElementMesh() {
372 MoFEMFunctionBegin;
373
374 const int max_level = defMaxLevel;
375
376 moab::Core core_ref;
377 moab::Interface &moab_ref = core_ref;
378
379 auto create_reference_element = [&moab_ref]() {
380 MoFEMFunctionBegin;
381 constexpr double base_coords[] = {
382
383 0, 0,
384 0,
385
386 1, 0,
387 0,
388
389 0, 1,
390 0
391
392 };
393 EntityHandle nodes[3];
394 for (int nn = 0; nn != 3; ++nn)
395 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
396 EntityHandle tri;
397 CHKERR moab_ref.create_element(MBTRI, nodes, 3, tri);
398
399 Range edges;
400 CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges,
401 moab::Interface::UNION);
402
403 MoFEMFunctionReturn(0);
404 };
405
406 CHKERR create_reference_element();
407
408 MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
409 MoFEM::Interface &m_field_ref = m_core_ref;
410
411 auto refine_ref_triangles = [this, &m_field_ref, max_level]() {
412 MoFEMFunctionBegin;
413 // seed ref mofem database by setting bit ref level to reference
414 // tetrahedron
415 CHKERR
416 m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
417 0, 2, BitRefLevel().set(0));
418
419 for (int ll = 0; ll != max_level; ++ll) {
420 MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "PostProc", "Refine Level %d",
421 ll);
422 Range edges;
423 CHKERR m_field_ref.getInterface<BitRefManager>()
424 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
425 BitRefLevel().set(), MBEDGE, edges);
426 Range tris;
427 CHKERR m_field_ref.getInterface<BitRefManager>()
428 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
429 BitRefLevel(ll).set(), MBTRI, tris);
430 // refine mesh
431 auto m_ref = m_field_ref.getInterface<MeshRefinement>();
432 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
433 BitRefLevel().set(ll + 1));
434 CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
435 }
436 MoFEMFunctionReturn(0);
437 };
438
439 // auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
440 // MoFEMFunctionBegin;
441 // Range faces;
442 // CHKERR moab_ref.get_entities_by_type(0, MBTRI, faces, true);
443 // Range faces_nodes;
444 // CHKERR moab_ref.get_connectivity(faces, faces_nodes, false);
445 // auto &gauss_pts = levelGaussPtsOnRefMesh[0];
446 // gauss_pts.resize(faces_nodes.size(), 4, false);
447 // size_t gg = 0;
448 // for (auto node : faces_nodes) {
449 // CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
450 // little_map[node] = gg;
451 // ++gg;
452 // }
453 // gauss_pts = trans(gauss_pts);
454 // MoFEMFunctionReturn(0);
455 // };
456
457 auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
458 &m_field_ref]() {
459 MoFEMFunctionBegin;
460 for (int ll = 0; ll != max_level + 1; ++ll) {
461
462 Range tris;
463 CHKERR
464 m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
465 BitRefLevel().set(ll), BitRefLevel().set(ll), MBTRI, tris);
466
467 if (hoNodes) {
468 EntityHandle meshset;
469 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
470 CHKERR moab_ref.add_entities(meshset, tris);
471 CHKERR moab_ref.convert_entities(meshset, true, false, false);
472 CHKERR moab_ref.delete_entities(&meshset, 1);
473 }
474
475 Range elem_nodes;
476 CHKERR moab_ref.get_connectivity(tris, elem_nodes, false);
477
478 auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
479 gauss_pts.resize(elem_nodes.size(), 3, false);
480 std::map<EntityHandle, int> little_map;
481 Range::iterator nit = elem_nodes.begin();
482 for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
483 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
484 little_map[*nit] = gg;
485 }
486 gauss_pts = trans(gauss_pts);
487
488 auto &ref_tris = levelRef[ll];
489 Range::iterator tit = tris.begin();
490 for (int tt = 0; tit != tris.end(); ++tit, ++tt) {
491 const EntityHandle *conn;
492 int num_nodes;
493 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
494 if (tt == 0) {
495 ref_tris.resize(tris.size(), num_nodes);
496 }
497 for (int nn = 0; nn != num_nodes; ++nn) {
498 ref_tris(tt, nn) = little_map[conn[nn]];
499 }
500 }
501
502 auto &shape_functions = levelShapeFunctions[ll];
503 shape_functions.resize(elem_nodes.size(), 3);
504 CHKERR ShapeMBTRI(&*shape_functions.data().begin(), &gauss_pts(0, 0),
505 &gauss_pts(1, 0), elem_nodes.size());
506 }
507 MoFEMFunctionReturn(0);
508 };
509
510 levelRef.resize(max_level + 1);
511 levelGaussPtsOnRefMesh.resize(max_level + 1);
512 levelShapeFunctions.resize(max_level + 1);
513
514 CHKERR refine_ref_triangles();
515 CHKERR get_ref_gauss_pts_and_shape_functions();
516
517 MoFEMFunctionReturn(0);
518}
519
520MoFEMErrorCode PostProcGenerateRefMesh<MBQUAD>::generateReferenceElementMesh() {
521 MoFEMFunctionBegin;
522
523#ifndef NDEBUG
524 if (defMaxLevel > 0)
525 MOFEM_LOG("WORLD", Sev::warning)
526 << "Refinement for quad is not implemented";
527#endif
528
529 moab::Core core_ref;
530 moab::Interface &moab_ref = core_ref;
531
532 auto create_reference_element = [&moab_ref]() {
533 MoFEMFunctionBegin;
534 constexpr double base_coords[] = {
535
536 0, 0,
537 0,
538
539 1, 0,
540 0,
541
542 1, 1,
543 0,
544
545 0, 1,
546 0
547
548 };
549 EntityHandle nodes[4];
550 for (int nn = 0; nn < 4; nn++)
551 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
552 EntityHandle quad;
553 CHKERR moab_ref.create_element(MBQUAD, nodes, 4, quad);
554 MoFEMFunctionReturn(0);
555 };
556
557 auto add_ho_nodes = [&]() {
558 MoFEMFunctionBegin;
559 Range quads;
560 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
561 EntityHandle meshset;
562 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
563 CHKERR moab_ref.add_entities(meshset, quads);
564 CHKERR moab_ref.convert_entities(meshset, true, true, true);
565 CHKERR moab_ref.delete_entities(&meshset, 1);
566 MoFEMFunctionReturn(0);
567 };
568
569 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
570 MoFEMFunctionBegin;
571 Range quads;
572 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
573 Range quads_nodes;
574 CHKERR moab_ref.get_connectivity(quads, quads_nodes, false);
575 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
576 gauss_pts.resize(quads_nodes.size(), 4, false);
577 size_t gg = 0;
578 for (auto node : quads_nodes) {
579 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
580 little_map[node] = gg;
581 ++gg;
582 }
583 gauss_pts = trans(gauss_pts);
584 MoFEMFunctionReturn(0);
585 };
586
587 auto set_ref_quads = [&](std::map<EntityHandle, int> &little_map) {
588 MoFEMFunctionBegin;
589 Range quads;
590 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
591 size_t hh = 0;
592 auto &ref_quads = levelRef[0];
593 for (auto quad : quads) {
594 const EntityHandle *conn;
595 int num_nodes;
596 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
597 if (ref_quads.size2() != num_nodes) {
598 ref_quads.resize(quads.size(), num_nodes);
599 }
600 for (int nn = 0; nn != num_nodes; ++nn) {
601 ref_quads(hh, nn) = little_map[conn[nn]];
602 }
603 ++hh;
604 }
605 MoFEMFunctionReturn(0);
606 };
607
608 auto set_shape_functions = [&]() {
609 MoFEMFunctionBegin;
610 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
611 auto &shape_functions = levelShapeFunctions[0];
612 const auto nb_gauss_pts = gauss_pts.size2();
613 shape_functions.resize(nb_gauss_pts, 4);
614 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
615 const double ksi = gauss_pts(0, gg);
616 const double zeta = gauss_pts(1, gg);
617 shape_functions(gg, 0) = N_MBQUAD0(ksi, zeta);
618 shape_functions(gg, 1) = N_MBQUAD1(ksi, zeta);
619 shape_functions(gg, 2) = N_MBQUAD2(ksi, zeta);
620 shape_functions(gg, 3) = N_MBQUAD3(ksi, zeta);
621 }
622 MoFEMFunctionReturn(0);
623 };
624
625 levelRef.resize(1);
626 levelGaussPtsOnRefMesh.resize(1);
627 levelShapeFunctions.resize(1);
628
629 CHKERR create_reference_element();
630 if (hoNodes)
631 CHKERR add_ho_nodes();
632 std::map<EntityHandle, int> little_map;
633 CHKERR set_gauss_pts(little_map);
634 CHKERR set_ref_quads(little_map);
635 CHKERR set_shape_functions();
636
637 MoFEMFunctionReturn(0);
638}
639
640MoFEMErrorCode PostProcGenerateRefMesh<MBEDGE>::generateReferenceElementMesh() {
641 MoFEMFunctionBegin;
642
643#ifndef NDEBUG
644 if (defMaxLevel > 0)
645 MOFEM_LOG("WORLD", Sev::warning)
646 << "Refinement for edges is not implemented";
647#endif
648
649 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
650 MoFEMFunctionBegin;
651
652 int nb_nodes = 2;
653 if (hoNodes)
654 nb_nodes = 3;
655
656 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
657 gauss_pts.resize(2, nb_nodes, false);
658 gauss_pts.clear();
659
660 int nn = 0;
661 for (; nn != 2; ++nn) {
662 gauss_pts(0, nn) = static_cast<double>(nn);
663 little_map[nn] = nn;
664 }
665
666 if (nn < nb_nodes) {
667 gauss_pts(0, nn) = 0.5;
668 little_map[nn] = 2;
669 }
670
671 MoFEMFunctionReturn(0);
672 };
673
674 auto set_ref_edges = [&](std::map<EntityHandle, int> &little_map) {
675 MoFEMFunctionBegin;
676
677 int level = 0;
678 int nb_edges = level + 1;
679
680 int nb_nodes = 2;
681 if (hoNodes)
682 nb_nodes = 3;
683
684 auto &ref_edges = levelRef[level];
685 ref_edges.resize(nb_edges, nb_nodes, false);
686
687 for (int ee = 0; ee != nb_edges; ++ee) {
688 int nn = 0;
689 for (; nn != 2; ++nn) {
690 ref_edges(ee, nn) = nb_nodes * ee + nn;
691 }
692 if (nn < nb_nodes) {
693 ref_edges(ee, nn) = nb_nodes * ee + 2;
694 }
695 }
696
697 MoFEMFunctionReturn(0);
698 };
699
700 auto set_shape_functions = [&]() {
701 MoFEMFunctionBegin;
702 auto &gauss_pts = levelGaussPtsOnRefMesh[0];
703 auto &shape_functions = levelShapeFunctions[0];
704 const auto nb_gauss_pts = gauss_pts.size2();
705 shape_functions.resize(nb_gauss_pts, 2);
706 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
707 const double ksi = gauss_pts(0, gg);
708 shape_functions(gg, 0) = N_MBEDGE0(ksi);
709 shape_functions(gg, 1) = N_MBEDGE1(ksi);
710 }
711 MoFEMFunctionReturn(0);
712 };
713
714 levelRef.resize(1);
715 levelGaussPtsOnRefMesh.resize(1);
716 levelShapeFunctions.resize(1);
717
718 std::map<EntityHandle, int> little_map;
719 CHKERR set_gauss_pts(little_map);
720 CHKERR set_ref_edges(little_map);
721 CHKERR set_shape_functions();
722
723 MoFEMFunctionReturn(0);
724}
725
726} // namespace MoFEM
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition LogManager.hpp:384
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition fem_tools.h:60
N_MBHEX7
#define N_MBHEX7(x, y, z)
Definition fem_tools.h:78
N_MBHEX3
#define N_MBHEX3(x, y, z)
Definition fem_tools.h:74
N_MBHEX5
#define N_MBHEX5(x, y, z)
Definition fem_tools.h:76
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition fem_tools.h:105
N_MBHEX4
#define N_MBHEX4(x, y, z)
Definition fem_tools.h:75
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition fem_tools.c:306
N_MBHEX0
#define N_MBHEX0(x, y, z)
Definition fem_tools.h:71
N_MBHEX6
#define N_MBHEX6(x, y, z)
Definition fem_tools.h:77
N_MBHEX2
#define N_MBHEX2(x, y, z)
Definition fem_tools.h:73
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition fem_tools.h:57
N_MBHEX1
#define N_MBHEX1(x, y, z)
Definition fem_tools.h:72
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition fem_tools.h:59
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition fem_tools.h:58
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition fem_tools.h:106
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
eta
double eta
Definition free_surface.cpp:192
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::MeshRefinement
Mesh refinement interface.
Definition MeshRefinement.hpp:26
MoFEM::PostProcGenerateRefMeshBase::getOptions
virtual MoFEMErrorCode getOptions(std::string prefix)
Definition PostProcBrokenMeshInMoabBase.cpp:17
MoFEM::PostProcGenerateRefMesh
Element for postprocessing. Uses MoAB to generate post-processing mesh.
Definition PostProcBrokenMeshInMoabBase.hpp:62
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
zeta
double zeta
Viscous hardening.
Definition plastic.cpp:131
Generated by Doxygen 1.9.8 and hosted at University of Glasgow