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 397 of file plot_base.cpp.

397  {
398  MoFEMFunctionBegin;
399  moab::Core core_ref;
400  moab::Interface &moab_ref = core_ref;
401 
402  char ref_mesh_file_name[255];
403 
404  if (SPACE_DIM == 2)
405  strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
406  else if (SPACE_DIM == 3)
407  strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
408  else
409  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
410  "Dimension not implemented");
411 
412  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-ref_file", ref_mesh_file_name,
413  255, PETSC_NULL);
414  CHKERR moab_ref.load_file(ref_mesh_file_name, 0, "");
415 
416  // Get elements
417  Range elems;
418  CHKERR moab_ref.get_entities_by_dimension(0, SPACE_DIM, elems);
419 
420  // Add mid-nodes on edges
421  EntityHandle meshset;
422  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
423  CHKERR moab_ref.add_entities(meshset, elems);
424  CHKERR moab_ref.convert_entities(meshset, true, false, false);
425  CHKERR moab_ref.delete_entities(&meshset, 1);
426 
427  // Get nodes on the mesh
428  Range elem_nodes;
429  CHKERR moab_ref.get_connectivity(elems, elem_nodes, false);
430 
431  // Map node entity and Gauss pint number
432  std::map<EntityHandle, int> nodes_pts_map;
433 
434  // Set gauss points coordinates from the reference mesh
435  gaussPts.resize(SPACE_DIM + 1, elem_nodes.size(), false);
436  gaussPts.clear();
437  Range::iterator nit = elem_nodes.begin();
438  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
439  double coords[3];
440  CHKERR moab_ref.get_coords(&*nit, 1, coords);
441  for (auto d : {0, 1, 2})
442  gaussPts(d, gg) = coords[d];
443  nodes_pts_map[*nit] = gg;
444  }
445 
446  if (SPACE_DIM == 2) {
447  // Set size of adjacency matrix (note ho order nodes 3 nodes and 3 nodes on
448  // edges)
449  refEleMap.resize(elems.size(), 3 + 3);
450  } else if (SPACE_DIM == 3) {
451  refEleMap.resize(elems.size(), 4 + 6);
452  }
453 
454  // Set adjacency matrix
455  Range::iterator tit = elems.begin();
456  for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
457  const EntityHandle *conn;
458  int num_nodes;
459  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
460  for (int nn = 0; nn != num_nodes; ++nn) {
461  refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
462  }
463  }
464 
465  MoFEMFunctionReturn(0);
466 }

◆ postProcess()

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

Definition at line 621 of file plot_base.cpp.

621  {
622  MoFEMFunctionBeginHot;
623 
624  auto resolve_shared_ents = [&]() {
625  MoFEMFunctionBegin;
626 
627  ParallelComm *pcomm_post_proc_mesh =
628  ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
629  if (pcomm_post_proc_mesh == NULL) {
630  // wrapRefMeshComm =
631  // boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
632  pcomm_post_proc_mesh = new ParallelComm(
633  &(getPostProcMesh()),
634  PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
635  }
636 
637  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
638 
639  MoFEMFunctionReturn(0);
640  };
641 
642  CHKERR resolve_shared_ents();
643 
644  MoFEMFunctionReturnHot(0);
645 }

◆ preProcess()

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

Definition at line 612 of file plot_base.cpp.

612  {
613  MoFEMFunctionBegin;
614  ParallelComm *pcomm_post_proc_mesh =
615  ParallelComm::get_pcomm(coreMeshPtr.get(), MYPCOMM_INDEX);
616  if (pcomm_post_proc_mesh != NULL)
617  delete pcomm_post_proc_mesh;
618  MoFEMFunctionReturn(0);
619 };

◆ 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 468 of file plot_base.cpp.

468  {
469  MoFEMFunctionBegin;
470 
471  const int num_nodes = gaussPts.size2();
472 
473  // Calculate shape functions
474 
475  switch (numeredEntFiniteElementPtr->getEntType()) {
476  case MBTRI:
477  shapeFunctions.resize(num_nodes, 3);
478  CHKERR Tools::shapeFunMBTRI(&*shapeFunctions.data().begin(),
479  &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
480  break;
481  case MBQUAD: {
482  shapeFunctions.resize(num_nodes, 4);
483  for (int gg = 0; gg != num_nodes; gg++) {
484  double ksi = gaussPts(0, gg);
485  double eta = gaussPts(1, gg);
486  shapeFunctions(gg, 0) = N_MBQUAD0(ksi, eta);
487  shapeFunctions(gg, 1) = N_MBQUAD1(ksi, eta);
488  shapeFunctions(gg, 2) = N_MBQUAD2(ksi, eta);
489  shapeFunctions(gg, 3) = N_MBQUAD3(ksi, eta);
490  }
491  } break;
492  case MBTET: {
493  shapeFunctions.resize(num_nodes, 8);
494  CHKERR Tools::shapeFunMBTET(&*shapeFunctions.data().begin(),
495  &gaussPts(0, 0), &gaussPts(1, 0),
496  &gaussPts(2, 0), num_nodes);
497  } break;
498  case MBHEX: {
499  shapeFunctions.resize(num_nodes, 8);
500  for (int gg = 0; gg != num_nodes; gg++) {
501  double ksi = gaussPts(0, gg);
502  double eta = gaussPts(1, gg);
503  double zeta = gaussPts(2, gg);
504  shapeFunctions(gg, 0) = N_MBHEX0(ksi, eta, zeta);
505  shapeFunctions(gg, 1) = N_MBHEX1(ksi, eta, zeta);
506  shapeFunctions(gg, 2) = N_MBHEX2(ksi, eta, zeta);
507  shapeFunctions(gg, 3) = N_MBHEX3(ksi, eta, zeta);
508  shapeFunctions(gg, 4) = N_MBHEX4(ksi, eta, zeta);
509  shapeFunctions(gg, 5) = N_MBHEX5(ksi, eta, zeta);
510  shapeFunctions(gg, 6) = N_MBHEX6(ksi, eta, zeta);
511  shapeFunctions(gg, 7) = N_MBHEX7(ksi, eta, zeta);
512  }
513  } break;
514  default:
515  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
516  "Not implemented element type");
517  }
518 
519  // Create physical nodes
520 
521  // MoAB interface allowing for creating nodes and elements in the bulk
522  ReadUtilIface *iface;
523  CHKERR getPostProcMesh().query_interface(iface);
524 
525  std::vector<double *> arrays; /// pointers to memory allocated by MoAB for
526  /// storing X, Y, and Z coordinates
527  EntityHandle startv; // Starting handle for vertex
528  // Allocate memory for num_nodes, and return starting handle, and access to
529  // memort.
530  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
531 
532  mapGaussPts.resize(gaussPts.size2());
533  for (int gg = 0; gg != num_nodes; ++gg)
534  mapGaussPts[gg] = startv + gg;
535 
536  Tag th;
537  int def_in_the_loop = -1;
538  CHKERR getPostProcMesh().tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
539  th, MB_TAG_CREAT | MB_TAG_SPARSE,
540  &def_in_the_loop);
541 
542  // Create physical elements
543 
544  const int num_el = refEleMap.size1();
545  const int num_nodes_on_ele = refEleMap.size2();
546 
547  EntityHandle starte; // Starting handle to first created element
548  EntityHandle *conn; // Access to MOAB memory with connectivity of elements
549 
550  // Create tris/tets in the bulk in MoAB database
551  if (SPACE_DIM == 2)
552  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
553  starte, conn);
554  else if (SPACE_DIM == 3)
555  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
556  starte, conn);
557  else
558  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
559  "Dimension not implemented");
560 
561  // At this point elements (memory for elements) is allocated, at code bellow
562  // actual connectivity of elements is set.
563  for (unsigned int tt = 0; tt != refEleMap.size1(); ++tt) {
564  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
565  conn[num_nodes_on_ele * tt + nn] = mapGaussPts[refEleMap(tt, nn)];
566  }
567 
568  // Finalise elements creation. At that point MOAB updates adjacency tables,
569  // and elements are ready to use.
570  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
571 
572  auto physical_elements = Range(starte, starte + num_el - 1);
573  CHKERR getPostProcMesh().tag_clear_data(th, physical_elements, &(nInTheLoop));
574 
575  EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
576  int fe_num_nodes;
577  {
578  const EntityHandle *conn;
579  mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes, true);
580  coords.resize(3 * fe_num_nodes, false);
581  CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
582  }
583 
584  // Set physical coordinates to physical nodes
585  FTensor::Index<'i', 3> i;
586  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_n(
587  &*shapeFunctions.data().begin());
588 
589  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords(
590  arrays[0], arrays[1], arrays[2]);
591  const double *t_coords_ele_x = &coords[0];
592  const double *t_coords_ele_y = &coords[1];
593  const double *t_coords_ele_z = &coords[2];
594  for (int gg = 0; gg != num_nodes; ++gg) {
595  FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_ele_coords(
596  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
597  t_coords(i) = 0;
598  for (int nn = 0; nn != fe_num_nodes; ++nn) {
599  t_coords(i) += t_n * t_ele_coords(i);
600  for (auto ii : {0, 1, 2})
601  if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
602  t_coords(ii) = 0;
603  ++t_ele_coords;
604  ++t_n;
605  }
606  ++t_coords;
607  }
608 
609  MoFEMFunctionReturn(0);
610 }

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:447
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
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:177
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
N_MBHEX3
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
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:1094
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
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:440
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:416
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:346
Generated by Doxygen 1.8.17 and hosted at University of Glasgow