v0.14.0
Public Member Functions | Protected Attributes | List of all members
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  {
428  MoFEMFunctionBegin;
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_NULL, "", "-ref_file", ref_mesh_file_name,
443  255, PETSC_NULL);
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 
495  MoFEMFunctionReturn(0);
496 }

◆ postProcess()

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

Definition at line 651 of file plot_base.cpp.

651  {
652  MoFEMFunctionBeginHot;
653 
654  auto resolve_shared_ents = [&]() {
655  MoFEMFunctionBegin;
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 
669  MoFEMFunctionReturn(0);
670  };
671 
672  CHKERR resolve_shared_ents();
673 
674  MoFEMFunctionReturnHot(0);
675 }

◆ preProcess()

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

Definition at line 642 of file plot_base.cpp.

642  {
643  MoFEMFunctionBegin;
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;
648  MoFEMFunctionReturn(0);
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  {
499  MoFEMFunctionBegin;
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);
508  CHKERR Tools::shapeFunMBTRI(&*shapeFunctions.data().begin(),
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);
524  CHKERR Tools::shapeFunMBTET(&*shapeFunctions.data().begin(),
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;
616  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_n(
617  &*shapeFunctions.data().begin());
618 
619  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords(
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) {
625  FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_ele_coords(
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 
639  MoFEMFunctionReturn(0);
640 }

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:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
N_MBHEX5
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:76
N_MBHEX0
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
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:126
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
N_MBHEX3
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
SPACE_DIM
constexpr int SPACE_DIM
Definition: plot_base.cpp:26
Range
FTensor::Tensor0
Definition: Tensor0.hpp:16
MyPostProc::refEleMap
ublas::matrix< int > refEleMap
Definition: plot_base.cpp:44
N_MBHEX2
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
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
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
N_MBHEX6
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
MyPostProc::shapeFunctions
MatrixDouble shapeFunctions
Definition: plot_base.cpp:45
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Generated by Doxygen 1.8.17 and hosted at University of Glasgow