Generate reference mesh on single element.
Each element is subdivided on smaller elements, i.e. a reference mesh on single element is created. Nodes of such reference mesh are used as integration points at which field values are calculated and to each node a "moab" tag is attached to store those values.
301 {
303
304 auto get_nb_of_ref_levels_from_options = [this] {
306 int max_level = 0;
307 PetscBool flg = PETSC_TRUE;
309 "-my_max_post_proc_ref_level", &max_level, &flg);
310 return max_level;
311 } else {
313 }
314 return 0;
315 };
316
317 auto generate_for_hex = [&]() {
319
320 moab::Core core_ref;
321 moab::Interface &moab_ref = core_ref;
322
323 auto create_reference_element = [&moab_ref]() {
325 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
326 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
328 for (int nn = 0; nn < 8; nn++)
329 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
331 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
333 };
334
335 auto add_ho_nodes = [&]() {
338 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
340 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
341 CHKERR moab_ref.add_entities(meshset, hexes);
342 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
343 CHKERR moab_ref.delete_entities(&meshset, 1);
345 };
346
347 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
350 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
352 CHKERR moab_ref.get_connectivity(hexes, hexes_nodes,
false);
354 gauss_pts.resize(hexes_nodes.size(), 4, false);
355 size_t gg = 0;
356 for (auto node : hexes_nodes) {
357 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
358 little_map[node] = gg;
359 ++gg;
360 }
361 gauss_pts = trans(gauss_pts);
363 };
364
365 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
368 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
369 size_t hh = 0;
371 for (auto hex : hexes) {
373 int num_nodes;
374 CHKERR moab_ref.get_connectivity(hex, conn, num_nodes,
false);
375 if (ref_hexes.size2() != num_nodes) {
376 ref_hexes.resize(hexes.size(), num_nodes);
377 }
378 for (int nn = 0; nn != num_nodes; ++nn) {
379 ref_hexes(hh, nn) = little_map[conn[nn]];
380 }
381 ++hh;
382 }
384 };
385
386 auto set_shape_functions = [&]() {
390 const auto nb_gauss_pts = gauss_pts.size2();
391 shape_functions.resize(nb_gauss_pts, 8);
392 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
393 const double ksi = gauss_pts(0, gg);
394 const double zeta = gauss_pts(1, gg);
395 const double eta = gauss_pts(2, gg);
404 }
406 };
407
411
412 CHKERR create_reference_element();
415 std::map<EntityHandle, int> little_map;
416 CHKERR set_gauss_pts(little_map);
417 CHKERR set_ref_hexes(little_map);
418 CHKERR set_shape_functions();
419
421 };
422
423 auto generate_for_tet = [&]() {
425
426 const int max_level = get_nb_of_ref_levels_from_options();
427
428 moab::Core core_ref;
429 moab::Interface &moab_ref = core_ref;
430
431 auto create_reference_element = [&moab_ref]() {
433 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
435 for (int nn = 0; nn < 4; nn++) {
437 moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
438 }
440 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
442 };
443
446
447 auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
449
450
452 m_field_ref.
getInterface<BitRefManager>()->setBitRefLevelByDim(
454 for (int ll = 0; ll != max_level; ++ll) {
456 "Refine Level %d", ll);
459 ->getEntitiesByTypeAndRefLevel(
463 ->getEntitiesByTypeAndRefLevel(
465
466 MeshRefinement *m_ref;
468 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
471 }
473 };
474
475 auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
476 &m_field_ref]() {
478 for (int ll = 0; ll != max_level + 1; ++ll) {
482 ->getEntitiesByTypeAndRefLevel(
486 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
487 CHKERR moab_ref.add_entities(meshset, tets);
488 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
489 CHKERR moab_ref.delete_entities(&meshset, 1);
490 }
492 CHKERR moab_ref.get_connectivity(tets, elem_nodes,
false);
493
495 gauss_pts.resize(elem_nodes.size(), 4, false);
496 std::map<EntityHandle, int> little_map;
497 Range::iterator nit = elem_nodes.begin();
498 for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
499 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
500 little_map[*nit] = gg;
501 }
502 gauss_pts = trans(gauss_pts);
503
505 Range::iterator tit = tets.begin();
506 for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
508 int num_nodes;
509 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
510 if (tt == 0) {
511 ref_tets.resize(tets.size(), num_nodes);
512 }
513 for (int nn = 0; nn != num_nodes; ++nn) {
514 ref_tets(tt, nn) = little_map[conn[nn]];
515 }
516 }
517
519 shape_functions.resize(elem_nodes.size(), 4);
521 &gauss_pts(1, 0), &gauss_pts(2, 0),
522 elem_nodes.size());
523 }
525 };
526
530
531 CHKERR create_reference_element();
532 CHKERR refine_ref_tetrahedron();
533 CHKERR get_ref_gauss_pts_and_shape_functions();
534
536 };
537
538 CHKERR generate_for_hex();
539 CHKERR generate_for_tet();
540
542 }
#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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
double zeta
Viscous hardening.
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
std::vector< MatrixDouble > levelGaussPtsOnRefMeshTets
std::vector< ublas::matrix< int > > levelRefTets
std::vector< MatrixDouble > levelShapeFunctionsHexes
std::vector< MatrixDouble > levelGaussPtsOnRefMeshHexes
std::vector< MatrixDouble > levelShapeFunctionsTets
std::vector< ublas::matrix< int > > levelRefHexes