v0.15.0
Loading...
Searching...
No Matches
MyPostProc Struct Reference
Inheritance diagram for MyPostProc:
[legend]
Collaboration diagram for MyPostProc:
[legend]

Public Member Functions

MoFEMErrorCode generateReferenceElementMesh ()
 
MoFEMErrorCode setGaussPts (int order)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode postProcess ()
 

Protected Attributes

ublas::matrix< int > refEleMap
 
MatrixDouble shapeFunctions
 

Detailed Description

Examples
plot_base.cpp.

Definition at line 34 of file plot_base.cpp.

Member Function Documentation

◆ generateReferenceElementMesh()

MoFEMErrorCode MyPostProc::generateReferenceElementMesh ( )
Examples
plot_base.cpp.

Definition at line 427 of file plot_base.cpp.

427 {
429 moab::Core core_ref;
430 moab::Interface &moab_ref = core_ref;
431
432 char ref_mesh_file_name[255];
433
434 if (SPACE_DIM == 2)
435 strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
436 else if (SPACE_DIM == 3)
437 strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
438 else
439 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
440 "Dimension not implemented");
441
442 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-ref_file", ref_mesh_file_name,
443 255, PETSC_NULLPTR);
444 CHKERR moab_ref.load_file(ref_mesh_file_name, 0, "");
445
446 // Get elements
447 Range elems;
448 CHKERR moab_ref.get_entities_by_dimension(0, SPACE_DIM, elems);
449
450 // Add mid-nodes on edges
451 EntityHandle meshset;
452 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
453 CHKERR moab_ref.add_entities(meshset, elems);
454 CHKERR moab_ref.convert_entities(meshset, true, false, false);
455 CHKERR moab_ref.delete_entities(&meshset, 1);
456
457 // Get nodes on the mesh
458 Range elem_nodes;
459 CHKERR moab_ref.get_connectivity(elems, elem_nodes, false);
460
461 // Map node entity and Gauss pint number
462 std::map<EntityHandle, int> nodes_pts_map;
463
464 // Set gauss points coordinates from the reference mesh
465 gaussPts.resize(SPACE_DIM + 1, elem_nodes.size(), false);
466 gaussPts.clear();
467 Range::iterator nit = elem_nodes.begin();
468 for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
469 double coords[3];
470 CHKERR moab_ref.get_coords(&*nit, 1, coords);
471 for (auto d : {0, 1, 2})
472 gaussPts(d, gg) = coords[d];
473 nodes_pts_map[*nit] = gg;
474 }
475
476 if (SPACE_DIM == 2) {
477 // Set size of adjacency matrix (note ho order nodes 3 nodes and 3 nodes on
478 // edges)
479 refEleMap.resize(elems.size(), 3 + 3);
480 } else if (SPACE_DIM == 3) {
481 refEleMap.resize(elems.size(), 4 + 6);
482 }
483
484 // Set adjacency matrix
485 Range::iterator tit = elems.begin();
486 for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
487 const EntityHandle *conn;
488 int num_nodes;
489 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
490 for (int nn = 0; nn != num_nodes; ++nn) {
491 refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
492 }
493 }
494
496}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
constexpr int SPACE_DIM
Definition plot_base.cpp:26
ublas::matrix< int > refEleMap
Definition plot_base.cpp:44

◆ postProcess()

MoFEMErrorCode MyPostProc::postProcess ( )
Examples
plot_base.cpp.

Definition at line 651 of file plot_base.cpp.

651 {
653
654 auto resolve_shared_ents = [&]() {
656
657 ParallelComm *pcomm_post_proc_mesh =
658 ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
659 if (pcomm_post_proc_mesh == NULL) {
660 // wrapRefMeshComm =
661 // boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
662 pcomm_post_proc_mesh = new ParallelComm(
663 &(getPostProcMesh()),
664 PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
665 }
666
667 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
668
670 };
671
672 CHKERR resolve_shared_ents();
673
675}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

◆ preProcess()

MoFEMErrorCode MyPostProc::preProcess ( )
Examples
plot_base.cpp.

Definition at line 642 of file plot_base.cpp.

642 {
644 ParallelComm *pcomm_post_proc_mesh =
645 ParallelComm::get_pcomm(coreMeshPtr.get(), MYPCOMM_INDEX);
646 if (pcomm_post_proc_mesh != NULL)
647 delete pcomm_post_proc_mesh;
649};

◆ setGaussPts()

MoFEMErrorCode MyPostProc::setGaussPts ( int order)

pointers to memory allocated by MoAB for storing X, Y, and Z coordinates

Examples
plot_base.cpp.

Definition at line 498 of file plot_base.cpp.

498 {
500
501 const int num_nodes = gaussPts.size2();
502
503 // Calculate shape functions
504
505 switch (numeredEntFiniteElementPtr->getEntType()) {
506 case MBTRI:
507 shapeFunctions.resize(num_nodes, 3);
509 &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
510 break;
511 case MBQUAD: {
512 shapeFunctions.resize(num_nodes, 4);
513 for (int gg = 0; gg != num_nodes; gg++) {
514 double ksi = gaussPts(0, gg);
515 double eta = gaussPts(1, gg);
516 shapeFunctions(gg, 0) = N_MBQUAD0(ksi, eta);
517 shapeFunctions(gg, 1) = N_MBQUAD1(ksi, eta);
518 shapeFunctions(gg, 2) = N_MBQUAD2(ksi, eta);
519 shapeFunctions(gg, 3) = N_MBQUAD3(ksi, eta);
520 }
521 } break;
522 case MBTET: {
523 shapeFunctions.resize(num_nodes, 8);
525 &gaussPts(0, 0), &gaussPts(1, 0),
526 &gaussPts(2, 0), num_nodes);
527 } break;
528 case MBHEX: {
529 shapeFunctions.resize(num_nodes, 8);
530 for (int gg = 0; gg != num_nodes; gg++) {
531 double ksi = gaussPts(0, gg);
532 double eta = gaussPts(1, gg);
533 double zeta = gaussPts(2, gg);
534 shapeFunctions(gg, 0) = N_MBHEX0(ksi, eta, zeta);
535 shapeFunctions(gg, 1) = N_MBHEX1(ksi, eta, zeta);
536 shapeFunctions(gg, 2) = N_MBHEX2(ksi, eta, zeta);
537 shapeFunctions(gg, 3) = N_MBHEX3(ksi, eta, zeta);
538 shapeFunctions(gg, 4) = N_MBHEX4(ksi, eta, zeta);
539 shapeFunctions(gg, 5) = N_MBHEX5(ksi, eta, zeta);
540 shapeFunctions(gg, 6) = N_MBHEX6(ksi, eta, zeta);
541 shapeFunctions(gg, 7) = N_MBHEX7(ksi, eta, zeta);
542 }
543 } break;
544 default:
545 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
546 "Not implemented element type");
547 }
548
549 // Create physical nodes
550
551 // MoAB interface allowing for creating nodes and elements in the bulk
552 ReadUtilIface *iface;
553 CHKERR getPostProcMesh().query_interface(iface);
554
555 std::vector<double *> arrays; /// pointers to memory allocated by MoAB for
556 /// storing X, Y, and Z coordinates
557 EntityHandle startv; // Starting handle for vertex
558 // Allocate memory for num_nodes, and return starting handle, and access to
559 // memort.
560 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
561
562 mapGaussPts.resize(gaussPts.size2());
563 for (int gg = 0; gg != num_nodes; ++gg)
564 mapGaussPts[gg] = startv + gg;
565
566 Tag th;
567 int def_in_the_loop = -1;
568 CHKERR getPostProcMesh().tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
569 th, MB_TAG_CREAT | MB_TAG_SPARSE,
570 &def_in_the_loop);
571
572 // Create physical elements
573
574 const int num_el = refEleMap.size1();
575 const int num_nodes_on_ele = refEleMap.size2();
576
577 EntityHandle starte; // Starting handle to first created element
578 EntityHandle *conn; // Access to MOAB memory with connectivity of elements
579
580 // Create tris/tets in the bulk in MoAB database
581 if (SPACE_DIM == 2)
582 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
583 starte, conn);
584 else if (SPACE_DIM == 3)
585 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
586 starte, conn);
587 else
588 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
589 "Dimension not implemented");
590
591 // At this point elements (memory for elements) is allocated, at code bellow
592 // actual connectivity of elements is set.
593 for (unsigned int tt = 0; tt != refEleMap.size1(); ++tt) {
594 for (int nn = 0; nn != num_nodes_on_ele; ++nn)
595 conn[num_nodes_on_ele * tt + nn] = mapGaussPts[refEleMap(tt, nn)];
596 }
597
598 // Finalise elements creation. At that point MOAB updates adjacency tables,
599 // and elements are ready to use.
600 CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
601
602 auto physical_elements = Range(starte, starte + num_el - 1);
603 CHKERR getPostProcMesh().tag_clear_data(th, physical_elements, &(nInTheLoop));
604
605 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
606 int fe_num_nodes;
607 {
608 const EntityHandle *conn;
609 mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes, true);
610 coords.resize(3 * fe_num_nodes, false);
611 CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
612 }
613
614 // Set physical coordinates to physical nodes
615 FTensor::Index<'i', 3> i;
617 &*shapeFunctions.data().begin());
618
620 arrays[0], arrays[1], arrays[2]);
621 const double *t_coords_ele_x = &coords[0];
622 const double *t_coords_ele_y = &coords[1];
623 const double *t_coords_ele_z = &coords[2];
624 for (int gg = 0; gg != num_nodes; ++gg) {
626 t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
627 t_coords(i) = 0;
628 for (int nn = 0; nn != fe_num_nodes; ++nn) {
629 t_coords(i) += t_n * t_ele_coords(i);
630 for (auto ii : {0, 1, 2})
631 if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
632 t_coords(ii) = 0;
633 ++t_ele_coords;
634 ++t_n;
635 }
636 ++t_coords;
637 }
638
640}
#define N_MBQUAD3(x, y)
quad shape function
Definition fem_tools.h:60
#define N_MBHEX7(x, y, z)
Definition fem_tools.h:78
#define N_MBHEX3(x, y, z)
Definition fem_tools.h:74
#define N_MBHEX5(x, y, z)
Definition fem_tools.h:76
#define N_MBHEX4(x, y, z)
Definition fem_tools.h:75
#define N_MBHEX0(x, y, z)
Definition fem_tools.h:71
#define N_MBHEX6(x, y, z)
Definition fem_tools.h:77
#define N_MBHEX2(x, y, z)
Definition fem_tools.h:73
#define N_MBQUAD0(x, y)
quad shape function
Definition fem_tools.h:57
#define N_MBHEX1(x, y, z)
Definition fem_tools.h:72
#define N_MBQUAD2(x, y)
quad shape function
Definition fem_tools.h:59
#define N_MBQUAD1(x, y)
quad shape function
Definition fem_tools.h:58
double eta
FTensor::Index< 'i', SPACE_DIM > i
double zeta
Viscous hardening.
Definition plastic.cpp:130
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition Tools.hpp:747
static MoFEMErrorCode shapeFunMBTRI(double *shape, const double *ksi, const double *eta, const int nb)
Calculate shape functions on triangle.
Definition Tools.hpp:716
MatrixDouble shapeFunctions
Definition plot_base.cpp:45

Member Data Documentation

◆ refEleMap

ublas::matrix<int> MyPostProc::refEleMap
protected
Examples
plot_base.cpp.

Definition at line 44 of file plot_base.cpp.

◆ shapeFunctions

MatrixDouble MyPostProc::shapeFunctions
protected
Examples
plot_base.cpp.

Definition at line 45 of file plot_base.cpp.


The documentation for this struct was generated from the following file: