v0.14.0
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 
11 namespace MoFEM {
12 
13 PostProcGenerateRefMeshBase::PostProcGenerateRefMeshBase()
14  : hoNodes(PETSC_TRUE), defMaxLevel(0), countEle(0), countVertEle(0),
15  nbVertices(0) {}
16 
17 MoFEMErrorCode 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_NULL, opt1.c_str(), &defMaxLevel, PETSC_NULL);
23 
24  std::string opt2 = prefix.size() ? "-" + prefix + "_max_post_ho_nodes"
25  : "-max_post_ho_nodes";
26  CHKERR PetscOptionsGetBool(PETSC_NULL, opt2.c_str(), &hoNodes, PETSC_NULL);
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 
36 MoFEMErrorCode 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 
149 MoFEMErrorCode 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 
261 MoFEMErrorCode PostProcGenerateRefMesh<MBTRI>::generateReferenceElementMesh() {
262  MoFEMFunctionBegin;
263 
264  const int max_level = defMaxLevel;
265 
266  moab::Core core_ref;
267  moab::Interface &moab_ref = core_ref;
268 
269  auto create_reference_element = [&moab_ref]() {
270  MoFEMFunctionBegin;
271  constexpr double base_coords[] = {
272 
273  0, 0,
274  0,
275 
276  1, 0,
277  0,
278 
279  0, 1,
280  0
281 
282  };
283  EntityHandle nodes[3];
284  for (int nn = 0; nn != 3; ++nn)
285  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
286  EntityHandle tri;
287  CHKERR moab_ref.create_element(MBTRI, nodes, 3, tri);
288 
289  Range edges;
290  CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges,
291  moab::Interface::UNION);
292 
293  MoFEMFunctionReturn(0);
294  };
295 
296  CHKERR create_reference_element();
297 
298  MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
299  MoFEM::Interface &m_field_ref = m_core_ref;
300 
301  auto refine_ref_triangles = [this, &m_field_ref, max_level]() {
302  MoFEMFunctionBegin;
303  // seed ref mofem database by setting bit ref level to reference
304  // tetrahedron
305  CHKERR
306  m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
307  0, 2, BitRefLevel().set(0));
308 
309  for (int ll = 0; ll != max_level; ++ll) {
310  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "PostProc", "Refine Level %d",
311  ll);
312  Range edges;
313  CHKERR m_field_ref.getInterface<BitRefManager>()
314  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
315  BitRefLevel().set(), MBEDGE, edges);
316  Range tris;
317  CHKERR m_field_ref.getInterface<BitRefManager>()
318  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
319  BitRefLevel(ll).set(), MBTRI, tris);
320  // refine mesh
321  auto m_ref = m_field_ref.getInterface<MeshRefinement>();
322  CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
323  BitRefLevel().set(ll + 1));
324  CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
325  }
326  MoFEMFunctionReturn(0);
327  };
328 
329  // auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
330  // MoFEMFunctionBegin;
331  // Range faces;
332  // CHKERR moab_ref.get_entities_by_type(0, MBTRI, faces, true);
333  // Range faces_nodes;
334  // CHKERR moab_ref.get_connectivity(faces, faces_nodes, false);
335  // auto &gauss_pts = levelGaussPtsOnRefMesh[0];
336  // gauss_pts.resize(faces_nodes.size(), 4, false);
337  // size_t gg = 0;
338  // for (auto node : faces_nodes) {
339  // CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
340  // little_map[node] = gg;
341  // ++gg;
342  // }
343  // gauss_pts = trans(gauss_pts);
344  // MoFEMFunctionReturn(0);
345  // };
346 
347  auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
348  &m_field_ref]() {
349  MoFEMFunctionBegin;
350  for (int ll = 0; ll != max_level + 1; ++ll) {
351 
352  Range tris;
353  CHKERR
354  m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
355  BitRefLevel().set(ll), BitRefLevel().set(ll), MBTRI, tris);
356 
357  if (hoNodes) {
358  EntityHandle meshset;
359  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
360  CHKERR moab_ref.add_entities(meshset, tris);
361  CHKERR moab_ref.convert_entities(meshset, true, false, false);
362  CHKERR moab_ref.delete_entities(&meshset, 1);
363  }
364 
365  Range elem_nodes;
366  CHKERR moab_ref.get_connectivity(tris, elem_nodes, false);
367 
368  auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
369  gauss_pts.resize(elem_nodes.size(), 3, false);
370  std::map<EntityHandle, int> little_map;
371  Range::iterator nit = elem_nodes.begin();
372  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
373  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
374  little_map[*nit] = gg;
375  }
376  gauss_pts = trans(gauss_pts);
377 
378  auto &ref_tris = levelRef[ll];
379  Range::iterator tit = tris.begin();
380  for (int tt = 0; tit != tris.end(); ++tit, ++tt) {
381  const EntityHandle *conn;
382  int num_nodes;
383  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
384  if (tt == 0) {
385  ref_tris.resize(tris.size(), num_nodes);
386  }
387  for (int nn = 0; nn != num_nodes; ++nn) {
388  ref_tris(tt, nn) = little_map[conn[nn]];
389  }
390  }
391 
392  auto &shape_functions = levelShapeFunctions[ll];
393  shape_functions.resize(elem_nodes.size(), 3);
394  CHKERR ShapeMBTRI(&*shape_functions.data().begin(), &gauss_pts(0, 0),
395  &gauss_pts(1, 0), elem_nodes.size());
396  }
397  MoFEMFunctionReturn(0);
398  };
399 
400  levelRef.resize(max_level + 1);
401  levelGaussPtsOnRefMesh.resize(max_level + 1);
402  levelShapeFunctions.resize(max_level + 1);
403 
404  CHKERR refine_ref_triangles();
405  CHKERR get_ref_gauss_pts_and_shape_functions();
406 
407  MoFEMFunctionReturn(0);
408 }
409 
410 MoFEMErrorCode PostProcGenerateRefMesh<MBQUAD>::generateReferenceElementMesh() {
411  MoFEMFunctionBegin;
412 
413 #ifndef NDEBUG
414  if (defMaxLevel > 0)
415  MOFEM_LOG("WORLD", Sev::warning)
416  << "Refinement for quad is not implemented";
417 #endif
418 
419  moab::Core core_ref;
420  moab::Interface &moab_ref = core_ref;
421 
422  auto create_reference_element = [&moab_ref]() {
423  MoFEMFunctionBegin;
424  constexpr double base_coords[] = {
425 
426  0, 0,
427  0,
428 
429  1, 0,
430  0,
431 
432  1, 1,
433  0,
434 
435  0, 1,
436  0
437 
438  };
439  EntityHandle nodes[4];
440  for (int nn = 0; nn < 4; nn++)
441  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
442  EntityHandle quad;
443  CHKERR moab_ref.create_element(MBQUAD, nodes, 4, quad);
444  MoFEMFunctionReturn(0);
445  };
446 
447  auto add_ho_nodes = [&]() {
448  MoFEMFunctionBegin;
449  Range quads;
450  CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
451  EntityHandle meshset;
452  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
453  CHKERR moab_ref.add_entities(meshset, quads);
454  CHKERR moab_ref.convert_entities(meshset, true, true, true);
455  CHKERR moab_ref.delete_entities(&meshset, 1);
456  MoFEMFunctionReturn(0);
457  };
458 
459  auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
460  MoFEMFunctionBegin;
461  Range quads;
462  CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
463  Range quads_nodes;
464  CHKERR moab_ref.get_connectivity(quads, quads_nodes, false);
465  auto &gauss_pts = levelGaussPtsOnRefMesh[0];
466  gauss_pts.resize(quads_nodes.size(), 4, false);
467  size_t gg = 0;
468  for (auto node : quads_nodes) {
469  CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
470  little_map[node] = gg;
471  ++gg;
472  }
473  gauss_pts = trans(gauss_pts);
474  MoFEMFunctionReturn(0);
475  };
476 
477  auto set_ref_quads = [&](std::map<EntityHandle, int> &little_map) {
478  MoFEMFunctionBegin;
479  Range quads;
480  CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads, true);
481  size_t hh = 0;
482  auto &ref_quads = levelRef[0];
483  for (auto quad : quads) {
484  const EntityHandle *conn;
485  int num_nodes;
486  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
487  if (ref_quads.size2() != num_nodes) {
488  ref_quads.resize(quads.size(), num_nodes);
489  }
490  for (int nn = 0; nn != num_nodes; ++nn) {
491  ref_quads(hh, nn) = little_map[conn[nn]];
492  }
493  ++hh;
494  }
495  MoFEMFunctionReturn(0);
496  };
497 
498  auto set_shape_functions = [&]() {
499  MoFEMFunctionBegin;
500  auto &gauss_pts = levelGaussPtsOnRefMesh[0];
501  auto &shape_functions = levelShapeFunctions[0];
502  const auto nb_gauss_pts = gauss_pts.size2();
503  shape_functions.resize(nb_gauss_pts, 4);
504  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
505  const double ksi = gauss_pts(0, gg);
506  const double zeta = gauss_pts(1, gg);
507  shape_functions(gg, 0) = N_MBQUAD0(ksi, zeta);
508  shape_functions(gg, 1) = N_MBQUAD1(ksi, zeta);
509  shape_functions(gg, 2) = N_MBQUAD2(ksi, zeta);
510  shape_functions(gg, 3) = N_MBQUAD3(ksi, zeta);
511  }
512  MoFEMFunctionReturn(0);
513  };
514 
515  levelRef.resize(1);
516  levelGaussPtsOnRefMesh.resize(1);
517  levelShapeFunctions.resize(1);
518 
519  CHKERR create_reference_element();
520  if (hoNodes)
521  CHKERR add_ho_nodes();
522  std::map<EntityHandle, int> little_map;
523  CHKERR set_gauss_pts(little_map);
524  CHKERR set_ref_quads(little_map);
525  CHKERR set_shape_functions();
526 
527  MoFEMFunctionReturn(0);
528 }
529 
530 MoFEMErrorCode PostProcGenerateRefMesh<MBEDGE>::generateReferenceElementMesh() {
531  MoFEMFunctionBegin;
532 
533 #ifndef NDEBUG
534  if (defMaxLevel > 0)
535  MOFEM_LOG("WORLD", Sev::warning)
536  << "Refinement for edges is not implemented";
537 #endif
538 
539  auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
540  MoFEMFunctionBegin;
541 
542  int nb_nodes = 2;
543  if (hoNodes)
544  nb_nodes = 3;
545 
546  auto &gauss_pts = levelGaussPtsOnRefMesh[0];
547  gauss_pts.resize(2, nb_nodes, false);
548  gauss_pts.clear();
549 
550  int nn = 0;
551  for (; nn != 2; ++nn) {
552  gauss_pts(0, nn) = static_cast<double>(nn);
553  little_map[nn] = nn;
554  }
555 
556  if (nn < nb_nodes) {
557  gauss_pts(0, nn) = 0.5;
558  little_map[nn] = 2;
559  }
560 
561  MoFEMFunctionReturn(0);
562  };
563 
564  auto set_ref_edges = [&](std::map<EntityHandle, int> &little_map) {
565  MoFEMFunctionBegin;
566 
567  int level = 0;
568  int nb_edges = level + 1;
569 
570  int nb_nodes = 2;
571  if (hoNodes)
572  nb_nodes = 3;
573 
574  auto &ref_edges = levelRef[level];
575  ref_edges.resize(nb_edges, nb_nodes, false);
576 
577  for (int ee = 0; ee != nb_edges; ++ee) {
578  int nn = 0;
579  for (; nn != 2; ++nn) {
580  ref_edges(ee, nn) = nb_nodes * ee + nn;
581  }
582  if (nn < nb_nodes) {
583  ref_edges(ee, nn) = nb_nodes * ee + 2;
584  }
585  }
586 
587  MoFEMFunctionReturn(0);
588  };
589 
590  auto set_shape_functions = [&]() {
591  MoFEMFunctionBegin;
592  auto &gauss_pts = levelGaussPtsOnRefMesh[0];
593  auto &shape_functions = levelShapeFunctions[0];
594  const auto nb_gauss_pts = gauss_pts.size2();
595  shape_functions.resize(nb_gauss_pts, 2);
596  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
597  const double ksi = gauss_pts(0, gg);
598  shape_functions(gg, 0) = N_MBEDGE0(ksi);
599  shape_functions(gg, 1) = N_MBEDGE1(ksi);
600  }
601  MoFEMFunctionReturn(0);
602  };
603 
604  levelRef.resize(1);
605  levelGaussPtsOnRefMesh.resize(1);
606  levelShapeFunctions.resize(1);
607 
608  std::map<EntityHandle, int> little_map;
609  CHKERR set_gauss_pts(little_map);
610  CHKERR set_ref_edges(little_map);
611  CHKERR set_shape_functions();
612 
613  MoFEMFunctionReturn(0);
614 }
615 
616 } // namespace MoFEM
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
EntityHandle
N_MBHEX5
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:76
MoFEM::PostProcGenerateRefMeshBase::PostProcGenerateRefMeshBase
PostProcGenerateRefMeshBase()
Definition: PostProcBrokenMeshInMoabBase.cpp:13
MoFEM::PostProcGenerateRefMesh
Element for postprocessing. Uses MoAB to generate post-processing mesh.
Definition: PostProcBrokenMeshInMoabBase.hpp:62
N_MBHEX0
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM.hpp
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
N_MBHEX7
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:78
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
N_MBHEX3
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
MoFEM::PostProcGenerateRefMeshBase::getOptions
virtual MoFEMErrorCode getOptions(std::string prefix)
Definition: PostProcBrokenMeshInMoabBase.cpp:17
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::CoreTmp
Definition: Core.hpp:36
MoFEM::PostProcGenerateRefMeshBase::hoNodes
PetscBool hoNodes
Definition: PostProcBrokenMeshInMoabBase.hpp:49
eta
double eta
Definition: free_surface.cpp:170
N_MBHEX1
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:72
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
MoFEM::PostProcGenerateRefMeshBase::defMaxLevel
int defMaxLevel
Definition: PostProcBrokenMeshInMoabBase.hpp:50
Range
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
N_MBHEX2
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
N_MBHEX4
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:75
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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
N_MBHEX6
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
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
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow