v0.13.0
Classes | Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
PostProcFaceOnRefinedMesh Struct Reference

Postprocess on face. More...

#include <users_modules/basic_finite_elements/src/PostProcOnRefMesh.hpp>

Inheritance diagram for PostProcFaceOnRefinedMesh:
[legend]
Collaboration diagram for PostProcFaceOnRefinedMesh:
[legend]

Classes

struct  CommonData
 
struct  OpGetFieldValuesOnSkinImpl
 

Public Types

typedef struct OpGetFieldValuesOnSkinImpl< 3 > OpGetFieldGradientValuesOnSkin
 
typedef struct OpGetFieldValuesOnSkinImpl< 1 > OpGetFieldValuesOnSkin
 

Public Member Functions

 PostProcFaceOnRefinedMesh (MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
 
int getRule (int order)
 
MoFEMErrorCode generateReferenceElementMesh ()
 
MoFEMErrorCode setGaussPts (int order)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode postProcess ()
 
virtual PostProcCommonOnRefMesh::CommonDatagetCommonData ()
 
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin (const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
 
MoFEMErrorCode addFieldValuesPostProcOnSkin (const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > mat_ptr=nullptr, bool save_on_tag=true)
 
- Public Member Functions inherited from PostProcTemplateOnRefineMesh< MoFEM::FaceElementForcesAndSourcesCore >
 PostProcTemplateOnRefineMesh (MoFEM::Interface &m_field)
 
virtual ~PostProcTemplateOnRefineMesh ()
 
MoFEMErrorCode addFieldValuesPostProc (const std::string field_name, Vec v=PETSC_NULL)
 Add operator to post-process L2, H1, Hdiv, Hcurl field value. More...
 
MoFEMErrorCode addFieldValuesPostProc (const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
 Add operator to post-process L2 or H1 field value. More...
 
MoFEMErrorCode addFieldValuesGradientPostProc (const std::string field_name, Vec v=PETSC_NULL)
 Add operator to post-process L2 or H1 field gradient. More...
 
MoFEMErrorCode addFieldValuesGradientPostProc (const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
 Add operator to post-process L2 or H1 field gradient. More...
 
MoFEMErrorCode addFieldValuesGradientPostProc (const std::string field_name, int space_dim, Vec v=PETSC_NULL)
 Add operator to post-process L2 or H1 field gradient. More...
 
MoFEMErrorCode writeFile (const std::string file_name)
 wrote results in (MOAB) format, use "file_name.h5m" More...
 

Public Attributes

bool sixNodePostProcTris
 
std::map< EntityHandle, EntityHandleelementsMap
 
CommonData commonData
 
- Public Attributes inherited from PostProcTemplateOnRefineMesh< MoFEM::FaceElementForcesAndSourcesCore >
moab::Core coreMesh
 
moab::Interface & postProcMesh
 
boost::shared_ptr< WrapMPIComm > wrapRefMeshComm
 
std::vector< EntityHandlemapGaussPts
 

Private Attributes

MatrixDouble gaussPtsTri
 Gauss points coordinates on reference triangle. More...
 
MatrixDouble gaussPtsQuad
 Gauss points coordinates on reference quad. More...
 
EntityHandletriConn
 Connectivity for created tri elements. More...
 
EntityHandlequadConn
 Connectivity for created quad elements. More...
 
EntityHandle startingVertTriHandle
 Starting handle for vertices on triangles. More...
 
EntityHandle startingVertQuadHandle
 Starting handle for vertices on quads. More...
 
std::vector< double * > verticesOnTriArrays
 
std::vector< double * > verticesOnQuadArrays
 
EntityHandle startingEleTriHandle
 Starting handle for triangles post proc. More...
 
EntityHandle startingEleQuadHandle
 Starting handle for quads post proc. More...
 
int numberOfTriangles
 Number of triangles to create. More...
 
int numberOfQuads
 NUmber of quads to create. More...
 
int counterTris
 
int counterQuads
 

Detailed Description

Postprocess on face.

Examples
approx_sphere.cpp, cell_forces.cpp, contact.cpp, eigen_elastic.cpp, elasticity.cpp, free_surface.cpp, heat_equation.cpp, minimal_surface_area.cpp, mortar_contact_thermal.cpp, nonlinear_elastic.cpp, plastic.cpp, plate.cpp, plot_base.cpp, poisson_2d_dis_galerkin.cpp, reaction_diffusion.cpp, shallow_wave.cpp, simple_contact.cpp, simple_contact_thermal.cpp, thermo_plastic.cpp, and wave_equation.cpp.

Definition at line 1052 of file PostProcOnRefMesh.hpp.

Member Typedef Documentation

◆ OpGetFieldGradientValuesOnSkin

Definition at line 1078 of file PostProcOnRefMesh.hpp.

◆ OpGetFieldValuesOnSkin

Definition at line 1078 of file PostProcOnRefMesh.hpp.

Constructor & Destructor Documentation

◆ PostProcFaceOnRefinedMesh()

PostProcFaceOnRefinedMesh::PostProcFaceOnRefinedMesh ( MoFEM::Interface m_field,
bool  six_node_post_proc_tris = true 
)

Member Function Documentation

◆ addFieldValuesGradientPostProcOnSkin()

MoFEMErrorCode PostProcFaceOnRefinedMesh::addFieldValuesGradientPostProcOnSkin ( const std::string  field_name,
const std::string  vol_fe_name = "dFE",
boost::shared_ptr< MatrixDouble >  grad_mat_ptr = nullptr,
bool  save_on_tag = true 
)

Definition at line 1004 of file PostProcOnRefMesh.cpp.

1006  {
1008 
1009  if (!grad_mat_ptr)
1010  grad_mat_ptr = boost::make_shared<MatrixDouble>();
1011 
1012  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
1013  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1014 
1015  if (mField.check_field("MESH_NODE_POSITIONS"))
1016  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *my_side_fe, true, false, false,
1017  false);
1018 
1019  // check number of coefficients
1020  auto field_ptr = mField.get_field_structure(field_name);
1021  const int nb_coefficients = field_ptr->getNbOfCoeffs();
1022 
1023  switch (nb_coefficients) {
1024  case 1:
1025  my_side_fe->getOpPtrVector().push_back(
1026  new OpCalculateScalarFieldGradient<3>(field_name, grad_mat_ptr));
1027  break;
1028  case 3:
1029  my_side_fe->getOpPtrVector().push_back(
1030  new OpCalculateVectorFieldGradient<3, 3>(field_name, grad_mat_ptr));
1031  break;
1032  case 6:
1033  my_side_fe->getOpPtrVector().push_back(
1035  grad_mat_ptr));
1036  break;
1037  default:
1038  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1039  "field with that number of coefficients is not implemented");
1040  }
1041 
1042  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1043  new OpGetFieldValuesOnSkinImpl<3>(
1044  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
1045  my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
1046 
1048 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get field gradients at integration pts for symmetric tensorial field rank 2.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh

◆ addFieldValuesPostProcOnSkin()

MoFEMErrorCode PostProcFaceOnRefinedMesh::addFieldValuesPostProcOnSkin ( const std::string  field_name,
const std::string  vol_fe_name = "dFE",
boost::shared_ptr< MatrixDouble >  mat_ptr = nullptr,
bool  save_on_tag = true 
)

Definition at line 1050 of file PostProcOnRefMesh.cpp.

1052  {
1054 
1055  if (!mat_ptr)
1056  mat_ptr = boost::make_shared<MatrixDouble>();
1057 
1058  auto my_side_fe =
1059  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1060 
1061  // check number of coefficients
1062  auto field_ptr = mField.get_field_structure(field_name);
1063  const int nb_coefficients = field_ptr->getNbOfCoeffs();
1064  FieldSpace space = field_ptr->getSpace();
1065 
1066  switch (space) {
1067  case L2:
1068  case H1:
1069  switch (nb_coefficients) {
1070  case 1:
1071  my_side_fe->getOpPtrVector().push_back(
1072  new OpCalculateVectorFieldValues<1>(field_name, mat_ptr));
1073  break;
1074  case 3:
1075  my_side_fe->getOpPtrVector().push_back(
1076  new OpCalculateVectorFieldValues<3>(field_name, mat_ptr));
1077  break;
1078  case 6:
1079  my_side_fe->getOpPtrVector().push_back(
1080  new OpCalculateTensor2SymmetricFieldValues<3>(field_name, mat_ptr));
1081  break;
1082  case 9:
1083  my_side_fe->getOpPtrVector().push_back(
1084  new OpCalculateTensor2FieldValues<3, 3>(field_name, mat_ptr));
1085  break;
1086  default:
1087  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1088  "field with that number of coefficients is not implemented");
1089  }
1090  break;
1091  case HDIV:
1092  switch (nb_coefficients) {
1093  case 1:
1094  my_side_fe->getOpPtrVector().push_back(
1095  new OpCalculateHVecVectorField<3>(field_name, mat_ptr));
1096  break;
1097  case 3:
1098  my_side_fe->getOpPtrVector().push_back(
1099  new OpCalculateHVecTensorField<3, 3>(field_name, mat_ptr));
1100  break;
1101  default:
1102  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1103  "field with that number of coefficients is not implemented");
1104  }
1105  break;
1106  default:
1107  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1108  "field with that space is not implemented.");
1109  }
1110 
1111  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1112  new OpGetFieldValuesOnSkinImpl<1>(postProcMesh, mapGaussPts, field_name,
1113  field_name, my_side_fe, vol_fe_name,
1114  mat_ptr, save_on_tag));
1115 
1117 }
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get vector field for H-div approximation.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Calculate symmetric tensor field values at integration pts.
Get values at integration pts for tensor filed rank 1, i.e. vector field.

◆ generateReferenceElementMesh()

MoFEMErrorCode PostProcFaceOnRefinedMesh::generateReferenceElementMesh ( )
Examples
minimal_surface_area.cpp.

Definition at line 552 of file PostProcOnRefMesh.cpp.

552  {
554 
555  auto generate_tri = [&](auto &gauss_pts) {
557  gauss_pts.resize(3, 3, false);
558  gauss_pts.clear();
559  gauss_pts(0, 0) = 0;
560  gauss_pts(1, 0) = 0;
561  gauss_pts(0, 1) = 1;
562  gauss_pts(1, 1) = 0;
563  gauss_pts(0, 2) = 0;
564  gauss_pts(1, 2) = 1;
565 
566  moab::Core core_ref;
567  moab::Interface &moab_ref = core_ref;
568  const EntityHandle *conn;
569  int num_nodes;
570  EntityHandle tri_conn[3];
571  MatrixDouble coords(6, 3);
572  for (int gg = 0; gg != 3; gg++) {
573  coords(gg, 0) = gauss_pts(0, gg);
574  coords(gg, 1) = gauss_pts(1, gg);
575  coords(gg, 2) = 0;
576  CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
577  }
578 
579  EntityHandle tri;
580  CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
581 
582  if (sixNodePostProcTris) {
583  Range edges;
584  CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges);
585  EntityHandle meshset;
586  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
587  CHKERR moab_ref.add_entities(meshset, &tri, 1);
588  CHKERR moab_ref.add_entities(meshset, edges);
589  CHKERR moab_ref.convert_entities(meshset, true, false, false);
590  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
591  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
592  gauss_pts.resize(3, num_nodes, false);
593  gauss_pts.clear();
594  for (int nn = 0; nn < num_nodes; nn++) {
595  gauss_pts(0, nn) = coords(nn, 0);
596  gauss_pts(1, nn) = coords(nn, 1);
597  }
598  } else {
599  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
600  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
601  gauss_pts.resize(3, num_nodes, false);
602  gauss_pts.clear();
603  for (int nn = 0; nn < 3; nn++) {
604  gauss_pts(0, nn) = coords(nn, 0);
605  gauss_pts(1, nn) = coords(nn, 1);
606  }
607  }
609  };
610 
611  auto generate_quad = [&](auto &gauss_pts) {
613 
614  gauss_pts.resize(3, 4, false);
615  gauss_pts.clear();
616 
617  gauss_pts(0, 0) = 0;
618  gauss_pts(1, 0) = 0;
619  gauss_pts(0, 1) = 1;
620  gauss_pts(1, 1) = 0;
621  gauss_pts(0, 2) = 1;
622  gauss_pts(1, 2) = 1;
623  gauss_pts(0, 3) = 0;
624  gauss_pts(1, 3) = 1;
625 
626  moab::Core core_ref;
627  moab::Interface &moab_ref = core_ref;
628 
629  const EntityHandle *conn;
630  int num_nodes;
631  std::array<EntityHandle, 4> quad_conn;
632  MatrixDouble coords(8, 3);
633  for (int gg = 0; gg != 4; gg++) {
634  coords(gg, 0) = gauss_pts(0, gg);
635  coords(gg, 1) = gauss_pts(1, gg);
636  coords(gg, 2) = 0;
637  CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
638  }
639 
640  EntityHandle quad;
641  CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
642 
643  if (sixNodePostProcTris) {
644  Range edges;
645  CHKERR moab_ref.get_adjacencies(&quad, 1, 1, true, edges);
646  EntityHandle meshset;
647  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
648  CHKERR moab_ref.add_entities(meshset, &quad, 1);
649  CHKERR moab_ref.add_entities(meshset, edges);
650  CHKERR moab_ref.convert_entities(meshset, true, false, false);
651  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
652  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
653  gauss_pts.resize(3, num_nodes, false);
654  gauss_pts.clear();
655  for (int nn = 0; nn != num_nodes; nn++) {
656  gauss_pts(0, nn) = coords(nn, 0);
657  gauss_pts(1, nn) = coords(nn, 1);
658  }
659  } else {
660  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
661  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
662  gauss_pts.resize(3, num_nodes, false);
663  gauss_pts.clear();
664  for (int nn = 0; nn != 4; nn++) {
665  gauss_pts(0, nn) = coords(nn, 0);
666  gauss_pts(1, nn) = coords(nn, 1);
667  }
668  }
669 
671  };
672 
673  CHKERR generate_tri(gaussPtsTri);
674  CHKERR generate_quad(gaussPtsQuad);
675 
677 }
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MatrixDouble gaussPtsQuad
Gauss points coordinates on reference quad.
MatrixDouble gaussPtsTri
Gauss points coordinates on reference triangle.

◆ getCommonData()

virtual PostProcCommonOnRefMesh::CommonData& PostProcFaceOnRefinedMesh::getCommonData ( )
virtual

Reimplemented from PostProcTemplateOnRefineMesh< MoFEM::FaceElementForcesAndSourcesCore >.

Definition at line 1078 of file PostProcOnRefMesh.hpp.

1078  {
1079  return commonData;
1080  }
CommonData commonData

◆ getRule()

int PostProcFaceOnRefinedMesh::getRule ( int  order)

Definition at line 1065 of file PostProcOnRefMesh.hpp.

1065 { return -1; };

◆ postProcess()

MoFEMErrorCode PostProcFaceOnRefinedMesh::postProcess ( )

Definition at line 890 of file PostProcOnRefMesh.cpp.

890  {
892 
893  auto update_elements = [&]() {
895  ReadUtilIface *iface;
896  CHKERR postProcMesh.query_interface(iface);
897 
898  if (counterTris)
899  CHKERR iface->update_adjacencies(startingEleTriHandle, counterTris,
900  gaussPtsTri.size2(), triConn);
901  if (counterQuads)
902  CHKERR iface->update_adjacencies(startingEleQuadHandle, counterQuads,
903  gaussPtsQuad.size2(), quadConn);
905  };
906 
907  auto resolve_shared = [&]() {
909  ParallelComm *pcomm_post_proc_mesh =
910  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
911  if (pcomm_post_proc_mesh == NULL) {
913  boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
914  pcomm_post_proc_mesh =
915  new ParallelComm(&postProcMesh, wrapRefMeshComm->get_comm());
916  }
917 
918  Range faces;
919  CHKERR postProcMesh.get_entities_by_dimension(0, 2, faces, false);
920  int rank = mField.get_comm_rank();
921  CHKERR postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(), faces,
922  &rank);
923  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
925  };
926 
927  CHKERR update_elements();
928  CHKERR resolve_shared();
929 
931 }
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
EntityHandle * quadConn
Connectivity for created quad elements.
EntityHandle startingEleQuadHandle
Starting handle for quads post proc.
EntityHandle startingEleTriHandle
Starting handle for triangles post proc.
EntityHandle * triConn
Connectivity for created tri elements.
boost::shared_ptr< WrapMPIComm > wrapRefMeshComm

◆ preProcess()

MoFEMErrorCode PostProcFaceOnRefinedMesh::preProcess ( )

Definition at line 813 of file PostProcOnRefMesh.cpp.

813  {
815 
816  ReadUtilIface *iface;
817  CHKERR postProcMesh.query_interface(iface);
818 
819  const int number_of_ents_in_the_loop = getLoopSize();
820  if (elementsMap.size() != number_of_ents_in_the_loop) {
821 
822  elementsMap.clear();
823  postProcMesh.delete_mesh();
824 
825  auto get_number_of_computational_elements = [&]() {
826  auto fe_name = this->feName;
827  auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
828 
829  auto miit =
830  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
831  boost::make_tuple(fe_name, this->getLoFERank()));
832  auto hi_miit =
833  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
834  boost::make_tuple(fe_name, this->getHiFERank()));
835 
836  const int number_of_ents_in_the_loop = this->getLoopSize();
837  if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
839  "Wrong size of indicices. Inconsistent size number of iterated "
840  "elements iterated by problem and from range.");
841  }
842 
843  std::array<int, MBMAXTYPE> nb_elemms_by_type;
844  std::fill(nb_elemms_by_type.begin(), nb_elemms_by_type.end(), 0);
845 
846  for (; miit != hi_miit; ++miit) {
847  auto type = (*miit)->getEntType();
848  ++nb_elemms_by_type[type];
849  }
850 
851  return nb_elemms_by_type;
852  };
853 
854  auto nb_computational_elements_by_type =
855  get_number_of_computational_elements();
856 
857  const int numberOfTriangles = nb_computational_elements_by_type[MBTRI];
858  const int numberOfQuads = nb_computational_elements_by_type[MBQUAD];
859 
860  // Here we create vertices using ReadUtilface
861  const int total_number_of_nodes_on_tri =
862  numberOfTriangles * gaussPtsTri.size2();
863  const int total_number_of_nodes_on_quad =
864  numberOfQuads * gaussPtsQuad.size2();
865 
866  if (total_number_of_nodes_on_tri) {
867  CHKERR iface->get_node_coords(3, total_number_of_nodes_on_tri, 0,
869  CHKERR iface->get_element_connect(numberOfTriangles, gaussPtsTri.size2(),
870  MBTRI, 0, startingEleTriHandle,
871  triConn);
872  }
873 
874  if (total_number_of_nodes_on_quad) {
875  CHKERR iface->get_node_coords(3, total_number_of_nodes_on_quad, 0,
878  CHKERR iface->get_element_connect(numberOfQuads, gaussPtsQuad.size2(),
879  MBQUAD, 0, startingEleQuadHandle,
880  quadConn);
881  }
882  }
883 
884  counterTris = 0;
885  counterQuads = 0;
886 
888 }
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
int numberOfTriangles
Number of triangles to create.
std::vector< double * > verticesOnTriArrays
int numberOfQuads
NUmber of quads to create.
std::map< EntityHandle, EntityHandle > elementsMap
std::vector< double * > verticesOnQuadArrays
EntityHandle startingVertQuadHandle
Starting handle for vertices on quads.
EntityHandle startingVertTriHandle
Starting handle for vertices on triangles.

◆ setGaussPts()

MoFEMErrorCode PostProcFaceOnRefinedMesh::setGaussPts ( int  order)

Definition at line 679 of file PostProcOnRefMesh.cpp.

679  {
681  if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
682  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
683  "post-process mesh not generated");
684 
685  auto create_tri = [&](auto &gauss_pts) {
686  std::array<EntityHandle, 3> tri_conn;
687  MatrixDouble3by3 coords_tri(3, 3);
688  CHKERR mField.get_moab().get_connectivity(
689  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
690  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_tri(0, 0));
691  const int num_nodes_on_ele = gauss_pts.size2();
692 
693  for (int gg = 0; gg != num_nodes_on_ele; gg++) {
694 
695  const double ksi = gauss_pts(0, gg);
696  const double eta = gauss_pts(1, gg);
697  const double n0 = N_MBTRI0(ksi, eta);
698  const double n1 = N_MBTRI1(ksi, eta);
699  const double n2 = N_MBTRI2(ksi, eta);
700  const double x =
701  n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
702  const double y =
703  n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
704  const double z =
705  n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
706 
707  verticesOnTriArrays[0][counterTris * num_nodes_on_ele + gg] = x;
708  verticesOnTriArrays[1][counterTris * num_nodes_on_ele + gg] = y;
709  verticesOnTriArrays[2][counterTris * num_nodes_on_ele + gg] = z;
710  }
711 
712  mapGaussPts.resize(num_nodes_on_ele);
713  for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
714  triConn[num_nodes_on_ele * counterTris + nn] =
715  num_nodes_on_ele * counterTris + nn + startingVertTriHandle;
716  mapGaussPts[nn] = triConn[num_nodes_on_ele * counterTris + nn];
717  }
718 
719  const auto tri = startingEleTriHandle + counterTris;
720 
721  return tri;
722  };
723 
724  auto create_quad = [&](auto &gauss_pts) {
725  std::array<EntityHandle, 4> quad_conn;
726  MatrixDouble coords_quad(4, 3);
727  CHKERR mField.get_moab().get_connectivity(
728  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
729  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_quad(0, 0));
730  const int num_nodes_on_ele = gauss_pts.size2();
731  for (int gg = 0; gg != num_nodes_on_ele; ++gg) {
732  double ksi = gauss_pts(0, gg);
733  double eta = gauss_pts(1, gg);
734  double n0 = N_MBQUAD0(ksi, eta);
735  double n1 = N_MBQUAD1(ksi, eta);
736  double n2 = N_MBQUAD2(ksi, eta);
737  double n3 = N_MBQUAD3(ksi, eta);
738  double x = n0 * coords_quad(0, 0) + n1 * coords_quad(1, 0) +
739  n2 * coords_quad(2, 0) + n3 * coords_quad(3, 0);
740  double y = n0 * coords_quad(0, 1) + n1 * coords_quad(1, 1) +
741  n2 * coords_quad(2, 1) + n3 * coords_quad(3, 1);
742  double z = n0 * coords_quad(0, 2) + n1 * coords_quad(1, 2) +
743  n2 * coords_quad(2, 2) + n3 * coords_quad(3, 2);
744  verticesOnQuadArrays[0][counterQuads * num_nodes_on_ele + gg] = x;
745  verticesOnQuadArrays[1][counterQuads * num_nodes_on_ele + gg] = y;
746  verticesOnQuadArrays[2][counterQuads * num_nodes_on_ele + gg] = z;
747  }
748 
749  mapGaussPts.resize(num_nodes_on_ele);
750 
751  for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
752  quadConn[num_nodes_on_ele * counterQuads + nn] =
753  num_nodes_on_ele * counterQuads + nn + startingVertQuadHandle;
754  mapGaussPts[nn] = quadConn[num_nodes_on_ele * counterQuads + nn];
755  }
756 
757  const auto quad = startingEleQuadHandle + counterQuads;
758  return quad;
759  };
760 
761  EntityHandle tri;
762 
763  if (elementsMap.size() == getLoopSize()) {
764  // Note "at" that will trigger error if element is not there.
765  tri = elementsMap.at(numeredEntFiniteElementPtr->getEnt());
766  switch (numeredEntFiniteElementPtr->getEntType()) {
767  case MBTRI:
768  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
769  noalias(gaussPts) = gaussPtsTri;
770  break;
771  case MBQUAD:
772  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
773  noalias(gaussPts) = gaussPtsQuad;
774  break;
775  default:
776  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
777  "Not implemented element type");
778  }
779 
780  // Set values which map nodes with integration points on the prism
781  const EntityHandle *conn;
782  int num_nodes;
783  CHKERR postProcMesh.get_connectivity(tri, conn, num_nodes, false);
784  mapGaussPts.resize(num_nodes);
785  for (int nn = 0; nn != num_nodes; nn++)
786  mapGaussPts[nn] = conn[nn];
787 
788  } else {
789  switch (numeredEntFiniteElementPtr->getEntType()) {
790  case MBTRI:
791  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
792  noalias(gaussPts) = gaussPtsTri;
793  tri = create_tri(gaussPtsTri);
794  ++counterTris;
795  break;
796  case MBQUAD:
797  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
798  noalias(gaussPts) = gaussPtsQuad;
799  tri = create_quad(gaussPtsQuad);
800  ++counterQuads;
801  break;
802  default:
803  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
804  "Not implemented element type");
805  }
806 
807  elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
808  }
809 
811 }
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:57
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:70
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:67
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:56
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:69
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:58
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:68
constexpr double eta
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:116

Member Data Documentation

◆ commonData

CommonData PostProcFaceOnRefinedMesh::commonData

Definition at line 1076 of file PostProcOnRefMesh.hpp.

◆ counterQuads

int PostProcFaceOnRefinedMesh::counterQuads
private

Definition at line 1149 of file PostProcOnRefMesh.hpp.

◆ counterTris

int PostProcFaceOnRefinedMesh::counterTris
private

Definition at line 1148 of file PostProcOnRefMesh.hpp.

◆ elementsMap

std::map<EntityHandle, EntityHandle> PostProcFaceOnRefinedMesh::elementsMap

Definition at line 1070 of file PostProcOnRefMesh.hpp.

◆ gaussPtsQuad

MatrixDouble PostProcFaceOnRefinedMesh::gaussPtsQuad
private

Gauss points coordinates on reference quad.

Definition at line 1127 of file PostProcOnRefMesh.hpp.

◆ gaussPtsTri

MatrixDouble PostProcFaceOnRefinedMesh::gaussPtsTri
private

Gauss points coordinates on reference triangle.

Definition at line 1126 of file PostProcOnRefMesh.hpp.

◆ numberOfQuads

int PostProcFaceOnRefinedMesh::numberOfQuads
private

NUmber of quads to create.

Definition at line 1147 of file PostProcOnRefMesh.hpp.

◆ numberOfTriangles

int PostProcFaceOnRefinedMesh::numberOfTriangles
private

Number of triangles to create.

Definition at line 1146 of file PostProcOnRefMesh.hpp.

◆ quadConn

EntityHandle* PostProcFaceOnRefinedMesh::quadConn
private

Connectivity for created quad elements.

Definition at line 1130 of file PostProcOnRefMesh.hpp.

◆ sixNodePostProcTris

bool PostProcFaceOnRefinedMesh::sixNodePostProcTris

Definition at line 1055 of file PostProcOnRefMesh.hpp.

◆ startingEleQuadHandle

EntityHandle PostProcFaceOnRefinedMesh::startingEleQuadHandle
private

Starting handle for quads post proc.

Definition at line 1144 of file PostProcOnRefMesh.hpp.

◆ startingEleTriHandle

EntityHandle PostProcFaceOnRefinedMesh::startingEleTriHandle
private

Starting handle for triangles post proc.

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

Definition at line 1143 of file PostProcOnRefMesh.hpp.

◆ startingVertQuadHandle

EntityHandle PostProcFaceOnRefinedMesh::startingVertQuadHandle
private

Starting handle for vertices on quads.

Definition at line 1134 of file PostProcOnRefMesh.hpp.

◆ startingVertTriHandle

EntityHandle PostProcFaceOnRefinedMesh::startingVertTriHandle
private

Starting handle for vertices on triangles.

Definition at line 1132 of file PostProcOnRefMesh.hpp.

◆ triConn

EntityHandle* PostProcFaceOnRefinedMesh::triConn
private

Connectivity for created tri elements.

Definition at line 1129 of file PostProcOnRefMesh.hpp.

◆ verticesOnQuadArrays

std::vector<double *> PostProcFaceOnRefinedMesh::verticesOnQuadArrays
private

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

Definition at line 1139 of file PostProcOnRefMesh.hpp.

◆ verticesOnTriArrays

std::vector<double *> PostProcFaceOnRefinedMesh::verticesOnTriArrays
private

Definition at line 1136 of file PostProcOnRefMesh.hpp.


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