v0.8.16
Public Types | Public Member Functions | Public Attributes | List of all members
MoFEM::TetGenInterface Struct Reference

use TetGen to generate mesh More...

#include <src/interfaces/TetGenInterface.hpp>

Inheritance diagram for MoFEM::TetGenInterface:
[legend]
Collaboration diagram for MoFEM::TetGenInterface:
[legend]

Public Types

enum  tetGenNodesTypes { RIDGEVERTEX = 0, FREESEGVERTEX = 1, FREEFACETVERTEX = 2, FREEVOLVERTEX = 3 }
 
typedef std::map< EntityHandle, unsigned long > moabTetGen_Map
 
typedef std::map< unsigned long, EntityHandletetGenMoab_Map
 
typedef std::map< int, Range > idxRange_Map
 

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 TetGenInterface (const MoFEM::Core &core)
 
MoFEMErrorCode inData (Range &ents, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Tag th=NULL)
 create TetGen data structure form range of moab entities More...
 
MoFEMErrorCode setGeomData (tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, std::map< int, Range > &type_ents)
 set point tags and type More...
 
MoFEMErrorCode outData (tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false)
 get entities for TetGen data structure More...
 
MoFEMErrorCode outData (tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, BitRefLevel bit, bool id_in_tags=false, bool error_if_created=false)
 get entities for TetGen data structure More...
 
MoFEMErrorCode setFaceData (std::vector< std::pair< Range, int > > &markers, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map)
 set markers to faces More...
 
MoFEMErrorCode getTriangleMarkers (tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL, bool only_non_zero=true)
 get markers to faces More...
 
MoFEMErrorCode setRegionData (std::vector< std::pair< EntityHandle, int > > &regions, tetgenio &in, Tag th=NULL)
 set region data to tetrahedral More...
 
MoFEMErrorCode getRegionData (tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL)
 get region data to tetrahedral More...
 
MoFEMErrorCode tetRahedralize (char switches[], tetgenio &in, tetgenio &out)
 run tetgen More...
 
MoFEMErrorCode loadPoly (char file_name[], tetgenio &in)
 load poly file More...
 
MoFEMErrorCode checkPlanar_Trinagle (double coords[], bool *result, const double eps=1e-9)
 
MoFEMErrorCode groupPlanar_Triangle (Range &tris, std::vector< Range > &sorted, const double eps=1e-9, Tag th=NULL)
 
MoFEMErrorCode groupRegion_Triangle (Range &tris, std::vector< std::vector< Range > > &sorted, const double eps=1e-9)
 Group surface triangles in planar regions. More...
 
MoFEMErrorCode makePolygonFacet (Range &ents, Range &polygons, bool reduce_edges=false, Range *not_reducable_nodes=NULL, const double eps=1e-9, Tag th=NULL)
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

use TetGen to generate mesh

Definition at line 36 of file TetGenInterface.hpp.

Member Typedef Documentation

◆ idxRange_Map

typedef std::map<int, Range> MoFEM::TetGenInterface::idxRange_Map

Definition at line 47 of file TetGenInterface.hpp.

◆ moabTetGen_Map

typedef std::map<EntityHandle, unsigned long> MoFEM::TetGenInterface::moabTetGen_Map

Definition at line 45 of file TetGenInterface.hpp.

◆ tetGenMoab_Map

typedef std::map<unsigned long, EntityHandle> MoFEM::TetGenInterface::tetGenMoab_Map

Definition at line 46 of file TetGenInterface.hpp.

Member Enumeration Documentation

◆ tetGenNodesTypes

Enumerator
RIDGEVERTEX 
FREESEGVERTEX 
FREEFACETVERTEX 
FREEVOLVERTEX 

Definition at line 62 of file TetGenInterface.hpp.

Constructor & Destructor Documentation

◆ TetGenInterface()

MoFEM::TetGenInterface::TetGenInterface ( const MoFEM::Core core)

Definition at line 42 of file TetGenInterface.hpp.

43  : cOre(const_cast<MoFEM::Core &>(core)) {}

Member Function Documentation

◆ checkPlanar_Trinagle()

MoFEMErrorCode MoFEM::TetGenInterface::checkPlanar_Trinagle ( double  coords[],
bool result,
const double  eps = 1e-9 
)

Definition at line 707 of file TetGenInterface.cpp.

709  {
711  double *pa = &coords[0];
712  double *pb = &coords[3];
713  double *pc = &coords[6];
714  double *pd = &coords[9];
715  double adx = pa[0] - pd[0];
716  double bdx = pb[0] - pd[0];
717  double cdx = pc[0] - pd[0];
718  double ady = pa[1] - pd[1];
719  double bdy = pb[1] - pd[1];
720  double cdy = pc[1] - pd[1];
721  double adz = pa[2] - pd[2];
722  double bdz = pb[2] - pd[2];
723  double cdz = pc[2] - pd[2];
724  double v = adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
725  cdx * (ady * bdz - adz * bdy);
726  double l = sqrt(pow(pa[0] - pb[0], 2) + pow(pa[1] - pb[1], 2) +
727  pow(pa[2] - pb[2], 2)) +
728  sqrt(pow(pa[0] - pc[0], 2) + pow(pa[1] - pc[1], 2) +
729  pow(pa[2] - pc[2], 2)) +
730  sqrt(pow(pa[0] - pd[0], 2) + pow(pa[1] - pd[1], 2) +
731  pow(pa[2] - pd[2], 2)) +
732  sqrt(pow(pb[0] - pc[0], 2) + pow(pb[1] - pc[1], 2) +
733  pow(pb[2] - pc[2], 2)) +
734  sqrt(pow(pb[0] - pd[0], 2) + pow(pb[1] - pd[1], 2) +
735  pow(pb[2] - pd[2], 2)) +
736  sqrt(pow(pc[0] - pd[0], 2) + pow(pc[1] - pd[1], 2) +
737  pow(pc[2] - pd[2], 2));
738  // std::cerr << fabs(v/pow(l,3)) << " ";
739  *result = fabs(v / pow(l, 3)) < eps ? true : false;
741  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static const double eps

◆ getRegionData()

MoFEMErrorCode MoFEM::TetGenInterface::getRegionData ( tetGenMoab_Map tetgen_moab_map,
tetgenio &  out,
Range *  ents = NULL,
idxRange_Map ents_map = NULL 
)

get region data to tetrahedral

Definition at line 648 of file TetGenInterface.cpp.

650  {
652  ErrorCode rval;
653  Interface &m_field = cOre;
654  int nbattributes = out.numberoftetrahedronattributes;
655  if (nbattributes == 0) {
656  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
657  "tetgen has no regions attributes");
658  }
659  Tag th_region;
660  rval = m_field.get_moab().tag_get_handle("TETGEN_REGION", th_region);
661  if (rval == MB_SUCCESS) {
662  rval = m_field.get_moab().tag_delete(th_region);
664  }
665  double def_marker = 0;
666  rval = m_field.get_moab().tag_get_handle(
667  "TETGEN_REGION", nbattributes, MB_TYPE_DOUBLE, th_region,
668  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
670  int ii = 0;
671  for (; ii < out.numberoftetrahedra; ii++) {
672  int jj = 0;
673  for (; jj < nbattributes; jj++) {
674  double id = out.tetrahedronattributelist[ii * nbattributes + jj];
675  int iii = MBTET | (ii << MB_TYPE_WIDTH);
676  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
677  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
678  "data inconsistency between TetGen and MoAB");
679  }
680  EntityHandle ent = tetgen_moab_map[iii];
681  rval = m_field.get_moab().tag_set_data(th_region, &ent, 1, &id);
683  if (ents != NULL)
684  ents->insert(ent);
685  if (ents_map != NULL)
686  (*ents_map)[id].insert(ent);
687  }
688  }
690  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MB_TYPE_WIDTH
Definition: definitions.h:285
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ getTriangleMarkers()

MoFEMErrorCode MoFEM::TetGenInterface::getTriangleMarkers ( tetGenMoab_Map tetgen_moab_map,
tetgenio &  out,
Range *  ents = NULL,
idxRange_Map ents_map = NULL,
bool  only_non_zero = true 
)

get markers to faces

Parameters
markersdata structure with markers
intetgen data structure (look to TetGen user manual)
moab_tetgen_mapmapping MoAB to TetGen entities
tetgen_moab_mapmapping TetGen to MoAB entities

Definition at line 550 of file TetGenInterface.cpp.

552  {
554  ErrorCode rval;
555  Interface &m_field = cOre;
556  Tag th_marker;
557  int def_marker = 0;
558  rval = m_field.get_moab().tag_get_handle(
559  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
560  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
562  int ii = 0;
563  for (; ii < out.numberoftrifaces; ii++) {
564  if (only_non_zero) {
565  if (out.trifacemarkerlist[ii] == 0) {
566  continue;
567  }
568  }
569  EntityHandle conn[3];
570  for (int nn = 0; nn < 3; nn++) {
571  int iii = MBVERTEX | (out.trifacelist[3 * ii + nn] << MB_TYPE_WIDTH);
572  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
573  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
574  "data inconsistency between TetGen and MoAB");
575  } else {
576  conn[nn] = tetgen_moab_map[iii];
577  }
578  }
579  Range face;
580  rval = m_field.get_moab().get_adjacencies(conn, 3, 2, false, face);
582  face = face.subset_by_type(MBTRI);
583  if (face.size() != 1) {
584  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
585  "data inconsistency between TetGen and MoAB, %u", face.size());
586  }
587  if (ents != NULL)
588  ents->merge(face);
589  if (ents_map != NULL)
590  (*ents_map)[out.trifacemarkerlist[ii]].merge(face);
591  rval = m_field.get_moab().tag_set_data(th_marker, &*face.begin(), 1,
592  &out.trifacemarkerlist[ii]);
594  }
596  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MB_TYPE_WIDTH
Definition: definitions.h:285
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ groupPlanar_Triangle()

MoFEMErrorCode MoFEM::TetGenInterface::groupPlanar_Triangle ( Range &  tris,
std::vector< Range > &  sorted,
const double  eps = 1e-9,
Tag  th = NULL 
)

Definition at line 743 of file TetGenInterface.cpp.

744  {
746 
747  Interface &m_field = cOre;
748 
749  ErrorCode rval;
750  Skinner skin(&m_field.get_moab());
751 
752  for (;;) {
753 
754  Range noplanar_to_anyother;
755  std::vector<Range>::iterator vit = sorted.begin();
756 
757  do {
758 
759  bool repeat = false;
760 
761  // get edges on vit skin
762  Range skin_edges;
763  rval = skin.find_skin(0, *vit, false, skin_edges);
765 
766  // skin edge nodes
767  Range skin_edges_nodes;
768  rval = m_field.get_moab().get_connectivity(skin_edges, skin_edges_nodes,
769  true);
771 
772  // get tris adjacent to vit skin edges
773  Range skin_edges_tris;
774  rval = m_field.get_moab().get_adjacencies(
775  skin_edges, 2, false, skin_edges_tris, moab::Interface::UNION);
777  // get tris which are part of facet
778  Range inner_tris = intersect(skin_edges_tris, *vit);
779  Range outer_tris = intersect(skin_edges_tris, tris);
780 
781  // tris coplanar with vit tris
782  Range coplanar_tris;
783 
784  Range::iterator tit = outer_tris.begin();
785  for (; tit != outer_tris.end(); tit++) {
786  Range tit_conn;
787  rval = m_field.get_moab().get_connectivity(&*tit, 1, tit_conn, true);
789  tit_conn = subtract(tit_conn, skin_edges_nodes);
790  if (tit_conn.empty()) {
791  coplanar_tris.insert(*tit);
792  noplanar_to_anyother.erase(*tit);
793  repeat = true;
794  } else {
795  Range tit_edges;
796  rval = m_field.get_moab().get_adjacencies(
797  &*tit, 1, 1, false, tit_edges, moab::Interface::UNION);
799  tit_edges = intersect(tit_edges, skin_edges);
800  if (tit_edges.size() != 1) {
801  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
802  "data inconsistency");
803  }
804  Range inner_tri;
805  rval = m_field.get_moab().get_adjacencies(
806  tit_edges, 2, false, inner_tri, moab::Interface::UNION);
808  inner_tri = intersect(inner_tri, inner_tris);
809  if (inner_tri.size() != 1) {
810  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
811  "data inconsistency");
812  }
813  // get connectivity
814  int num_nodes;
815  const EntityHandle *inner_tri_conn;
816  rval = m_field.get_moab().get_connectivity(
817  *inner_tri.begin(), inner_tri_conn, num_nodes, true);
819  double coords[12];
820  if(th) {
821  rval = m_field.get_moab().tag_get_data(th, inner_tri_conn, 3,
822  coords);
824  rval = m_field.get_moab().tag_get_data(th,&*tit_conn.begin(), 1,
825  &coords[9]);
827  } else {
828  rval = m_field.get_moab().get_coords(inner_tri_conn, 3, coords);
830  rval = m_field.get_moab().get_coords(&*tit_conn.begin(), 1,
831  &coords[9]);
833  }
834  bool coplanar;
835  ierr = checkPlanar_Trinagle(coords, &coplanar, eps);
836  CHKERRG(ierr);
837  if (coplanar) {
838  coplanar_tris.insert(*tit);
839  noplanar_to_anyother.erase(*tit);
840  repeat = true;
841  } else {
842  noplanar_to_anyother.insert(*tit);
843  }
844  }
845  }
846 
847  vit->merge(coplanar_tris);
848  tris = subtract(tris, *vit);
849 
850  if (repeat) {
851  vit = sorted.begin();
852  } else {
853  vit++;
854  }
855 
856  } while (vit != sorted.end());
857 
858  if (noplanar_to_anyother.empty()) {
860  } else {
861  Range seed;
862  seed.insert(noplanar_to_anyother[0]);
863  tris.erase(noplanar_to_anyother[0]);
864  sorted.push_back(seed);
865  }
866  }
867 
869  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:526
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
MoFEMErrorCode checkPlanar_Trinagle(double coords[], bool *result, const double eps=1e-9)
static const double eps

◆ groupRegion_Triangle()

MoFEMErrorCode MoFEM::TetGenInterface::groupRegion_Triangle ( Range &  tris,
std::vector< std::vector< Range > > &  sorted,
const double  eps = 1e-9 
)

Group surface triangles in planar regions.

Parameters
trisinput triangles
sortedoutput sorted planar faces
epstolerance

Definition at line 871 of file TetGenInterface.cpp.

872  {
874 
875  // PetscAttachDebugger();
876 
877  Range seed;
878  seed.insert(tris[0]);
879  tris.erase(tris[0]);
880  std::vector<Range> vec_seed;
881  vec_seed.push_back(seed);
882  sorted.push_back(vec_seed);
883 
884  for (;;) {
885  std::vector<Range> &vec = sorted.back();
886  ierr = groupPlanar_Triangle(tris, vec, eps);
887  CHKERRG(ierr);
888  if (tris.empty()) {
890  } else {
891  Range seed;
892  seed.insert(tris[0]);
893  tris.erase(tris[0]);
894  std::vector<Range> vec_seed;
895  vec_seed.push_back(seed);
896  sorted.push_back(vec_seed);
897  }
898  }
899 
901  }
MoFEMErrorCode groupPlanar_Triangle(Range &tris, std::vector< Range > &sorted, const double eps=1e-9, Tag th=NULL)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:526
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static const double eps

◆ inData()

MoFEMErrorCode MoFEM::TetGenInterface::inData ( Range &  ents,
tetgenio &  in,
moabTetGen_Map moab_tetgen_map,
tetGenMoab_Map tetgen_moab_map,
Tag  th = NULL 
)

create TetGen data structure form range of moab entities

Move mesh to TetGen data structures

Parameters
entsrange of entities (tetrahedrons or nodes)
intetgen data structure (look to TetGen user manual)
moab_tetgen_mapmapping moab to TetGen entities
tetgen_moab_mapmapping tetgen to moab entities

Definition at line 60 of file TetGenInterface.cpp.

63  {
65 
66  Interface &m_field = cOre;
67  Range::iterator it;
68 
69  //
70  ErrorCode rval;
71 
72  Tag th_marker;
73  int def_marker = 0;
74  rval = m_field.get_moab().tag_get_handle(
75  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
76  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
78 
79  // All indices start from 0
80  in.firstnumber = 0;
81 
82  Range points = ents.subset_by_dimension(0);
83  in.numberofpoints = points.size();
84  if (points.size() > 0) {
85  in.pointlist = new double[in.numberofpoints * 3];
86  in.pointmarkerlist = new int[in.numberofpoints];
87  if (th) {
88  rval = m_field.get_moab().tag_get_data(th, points, in.pointlist);
90  } else {
91  rval = m_field.get_moab().get_coords(points, in.pointlist);
93  }
94  rval = m_field.get_moab().tag_get_data(th_marker, points,
95  in.pointmarkerlist);
97  it = points.begin();
98  for (int ii = 0; it != points.end(); it++, ii++) {
99  unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
100  tetgen_moab_map[iii] = *it;
101  moab_tetgen_map[*it] = iii;
102  }
103  }
104 
105  Range tets = ents.subset_by_type(MBTET);
106  in.numberoftetrahedra = tets.size();
107  if (in.numberoftetrahedra > 0) {
108  in.tetrahedronlist = new int[4 * ents.subset_by_type(MBTET).size()];
109  it = tets.begin();
110  for (int ii = 0; it != tets.end(); it++, ii++) {
111  int num_nodes;
112  const EntityHandle *conn;
113  rval = m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
115 #ifdef DEBUG_TETGEN
116  {
117  double coords[12];
118  if (th) {
119  rval = m_field.tag_get_data(th, conn, num_nodes, coords);
121  } else {
122  rval = m_field.get_moab().get_coords(conn, num_nodes, coords);
124  }
125  double diff_n[12];
126  ShapeDiffMBTET(diff_n);
127  FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1],
128  &diff_n[2], 3);
129  FTensor::Tensor1<double *, 3> t_coords(&coords[0], &coords[1],
130  &coords[2], 3);
134  jac(i, j) = 0;
135  for (int nn = 0; nn != 4; nn++) {
136  jac(i, j) += t_coords(i) * t_diff_n(j);
137  ++t_coords;
138  ++t_diff_n;
139  }
140  double det;
141  determinantTensor3by3(jac, det);
142  if (det <= 0) {
143  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144  "Negative volume det = %6.4e", det);
145  }
146  }
147 #endif
148  tetgen_moab_map[MBTET | (ii << MB_TYPE_WIDTH)] = *it;
149  moab_tetgen_map[*it] = MBTET | (ii << MB_TYPE_WIDTH);
150  for (int nn = 0; nn != 4; nn++) {
151  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
152  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
153  "data inconsistency between TetGen and MoAB");
154  }
155  in.tetrahedronlist[4 * ii + nn] =
156  moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
157  }
158  }
159  }
160 
161  Range tris = ents.subset_by_type(MBTRI);
162  in.numberoftrifaces = tris.size();
163  if (in.numberoftrifaces) {
164  in.trifacelist = new int[3 * in.numberoftrifaces];
165  in.trifacemarkerlist = new int[in.numberoftrifaces];
166  // std::fill(&in.trifacemarkerlist[0],&in.trifacemarkerlist[in.numberoftrifaces],1);
167  rval = m_field.get_moab().tag_get_data(th_marker, tris,
168  in.trifacemarkerlist);
170  it = tris.begin();
171  for (int ii = 0; it != tris.end(); it++, ii++) {
172  int num_nodes;
173  const EntityHandle *conn;
174  rval = m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
176  int order[] = {0, 1, 2};
177  Range adj_tets;
178  rval = m_field.get_moab().get_adjacencies(&*it, 1, 3, true, adj_tets);
180  adj_tets = intersect(adj_tets, tets);
181  if (adj_tets.size() == 1) {
182  int side_number;
183  int sense;
184  int offset;
185  rval = m_field.get_moab().side_number(adj_tets[0], *it, side_number,
186  sense, offset);
188  if (sense == -1) {
189  order[0] = 1;
190  order[1] = 0;
191  }
192  }
193  tetgen_moab_map[MBTRI | (ii << MB_TYPE_WIDTH)] = *it;
194  moab_tetgen_map[*it] = MBTRI | (ii << MB_TYPE_WIDTH);
195  for (int nn = 0; nn < 3; nn++) {
196  if (moab_tetgen_map.find(conn[order[nn]]) == moab_tetgen_map.end()) {
197  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
198  "data inconsistency between TetGen and MoAB");
199  }
200  in.trifacelist[3 * ii + nn] =
201  moab_tetgen_map[conn[order[nn]]] >> MB_TYPE_WIDTH;
202  }
203  }
204  }
205 
206  Range edges = ents.subset_by_type(MBEDGE);
207  in.numberofedges = edges.size();
208  if (in.numberofedges > 0) {
209  in.edgelist = new int[2 * in.numberofedges];
210  in.edgemarkerlist = new int[in.numberofedges];
211  // std::fill(&in.edgemarkerlist[0],&in.edgemarkerlist[in.numberofedges],1);
212  rval =
213  m_field.get_moab().tag_get_data(th_marker, edges, in.edgemarkerlist);
215  it = edges.begin();
216  for (int ii = 0; it != edges.end(); it++, ii++) {
217  int num_nodes;
218  const EntityHandle *conn;
219  rval = m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
221  tetgen_moab_map[MBEDGE | (ii << MB_TYPE_WIDTH)] = *it;
222  moab_tetgen_map[*it] = MBEDGE | (ii << MB_TYPE_WIDTH);
223  for (int nn = 0; nn < 2; nn++) {
224  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
225  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
226  "data inconsistency between TetGen and MoAB");
227  }
228  in.edgelist[2 * ii + nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
229  }
230  }
231  }
232 
234  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant.
#define MB_TYPE_WIDTH
Definition: definitions.h:285
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:303
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ loadPoly()

MoFEMErrorCode MoFEM::TetGenInterface::loadPoly ( char  file_name[],
tetgenio &  in 
)

load poly file

Definition at line 699 of file TetGenInterface.cpp.

699  {
701  if (!in.load_poly(file_name)) {
702  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
703  "can not read TetGen poly file");
704  }
706  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ makePolygonFacet()

MoFEMErrorCode MoFEM::TetGenInterface::makePolygonFacet ( Range &  ents,
Range &  polygons,
bool  reduce_edges = false,
Range *  not_reducable_nodes = NULL,
const double  eps = 1e-9,
Tag  th = NULL 
)

make planar polygon facet

Parameters
entssurface triangles
polygonsoutput list of polygons
reduce_edgesreduce edges if on the line
not_reducable_nodesdo not reduce node on edge if in this range
epstolerance
Bug:
assumes that are no holes

Definition at line 902 of file TetGenInterface.cpp.

905  {
907  // FIXME: assumes that are no holes
908 
909  if (ents.empty()) {
910  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
911  "no ents to build polygon");
912  }
913 
914  Interface &m_field = cOre;
915 
916  //
917  ErrorCode rval;
918  Skinner skin(&m_field.get_moab());
919 
920  Range skin_edges;
921  rval = skin.find_skin(0, ents, false, skin_edges);
923 
924  std::vector<EntityHandle> polygon_nodes;
925  EntityHandle seed = skin_edges[0];
926  Range seen_edges;
927  seen_edges.insert(seed);
928  skin_edges.erase(seed);
929  int num_nodes;
930  const EntityHandle *conn;
931  rval = m_field.get_moab().get_connectivity(seed, conn, num_nodes, true);
933  polygon_nodes.push_back(conn[0]);
934  // std::cerr << std::endl;
935  // std::cerr << conn[0] << " " << conn[1] << std::endl;
936  do {
937  EntityHandle last_node = polygon_nodes.back();
938  Range adj_edges;
939  rval = m_field.get_moab().get_adjacencies(&last_node, 1, 1, false,
940  adj_edges);
942  adj_edges = intersect(adj_edges, skin_edges);
943  if (adj_edges.size() == 0) {
944  break;
945  }
946  if (adj_edges.size() != 1) {
947  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
948  "should be only one edge");
949  }
950  seen_edges.insert(adj_edges[0]);
951  skin_edges.erase(adj_edges[0]);
952  rval = m_field.get_moab().get_connectivity(adj_edges[0], conn, num_nodes,
953  true);
955  EntityHandle add_node = (last_node == conn[0]) ? conn[1] : conn[0];
956  polygon_nodes.push_back(add_node);
957  // std::cerr << "\t" << add_node << std::endl;
958  } while (1);
959 
960  if (reduce_edges) {
961  // std::cerr << "polygon " << polygon_nodes.size();
962  std::vector<EntityHandle>::iterator pit = polygon_nodes.begin();
963  // std::cerr << std::endl;
964  for (; pit != polygon_nodes.end();) {
965  if (not_reducable_nodes != NULL) {
966  if (not_reducable_nodes->find(*pit) != not_reducable_nodes->end()) {
967  pit++;
968  continue;
969  }
970  }
971  EntityHandle mm;
972  if (pit == polygon_nodes.begin()) {
973  mm = polygon_nodes.back();
974  } else {
975  mm = *(pit - 1);
976  }
977  EntityHandle mc = *pit;
978  EntityHandle mp;
979  if (polygon_nodes.back() == mc) {
980  mp = polygon_nodes[0];
981  } else {
982  mp = *(pit + 1);
983  }
984  double coords[9];
985  if (th) {
986  rval = m_field.get_moab().tag_get_data(th, &mm, 1, &coords[3 * 0]);
988  rval = m_field.get_moab().tag_get_data(th, &mc, 1, &coords[3 * 1]);
990  rval = m_field.get_moab().tag_get_data(th, &mp, 1, &coords[3 * 2]);
992  } else {
993  rval = m_field.get_moab().get_coords(&mm, 1, &coords[3 * 0]);
995  rval = m_field.get_moab().get_coords(&mc, 1, &coords[3 * 1]);
997  rval = m_field.get_moab().get_coords(&mp, 1, &coords[3 * 2]);
999  }
1000  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 0], 1); // mm = mm -
1001  // mc
1002  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 2], 1); // mp = mp -
1003  // mc
1004  double spin[9];
1005  ierr = Spin(spin, &coords[3 * 0]);
1006  CHKERRG(ierr);
1007  double l0 = cblas_dnrm2(3, &coords[3 * 0], 1);
1008  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1. / l0, spin, 3,
1009  &coords[3 * 2], 1, 0., &coords[3 * 1], 1);
1010  double dot = cblas_dnrm2(3, &coords[3 * 1], 1);
1011  // std::cerr << mm << " " << mc << " " << mp << " " << dot << std::endl;
1012  if (dot < eps) {
1013  polygon_nodes.erase(pit);
1014  pit = polygon_nodes.begin();
1015  // std::cerr << std::endl;
1016  } else {
1017  pit++;
1018  }
1019  }
1020  }
1021  // std::cerr << " " << polygon_nodes.size() << std::endl;
1022  /*pit = polygon_nodes.begin();
1023  for(;pit!=polygon_nodes.end();pit++) {
1024  double coords[3];
1025  rval = m_field.get_moab().get_coords(&*pit,1,coords); CHKERRQ_MOAB(rval);
1026  std::cerr << *pit << " " << coords[0] << " " << coords[1] << " " <<
1027  coords[2] << std::endl;
1028  }*/
1029 
1030  Range existing_polygon;
1031  rval = m_field.get_moab().get_adjacencies(
1032  &polygon_nodes[0], polygon_nodes.size(), 2, true, existing_polygon);
1033  CHKERRQ_MOAB(rval);
1034  if (existing_polygon.empty()) {
1035  EntityHandle polygon;
1036  rval = m_field.get_moab().create_element(MBPOLYGON, &polygon_nodes[0],
1037  polygon_nodes.size(), polygon);
1038  CHKERRQ_MOAB(rval);
1039  polygons.insert(polygon);
1040  } else {
1041  polygons.merge(existing_polygon);
1042  }
1043 
1045  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:480
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:526
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
static const double eps

◆ outData() [1/2]

MoFEMErrorCode MoFEM::TetGenInterface::outData ( tetgenio &  in,
tetgenio &  out,
moabTetGen_Map moab_tetgen_map,
tetGenMoab_Map tetgen_moab_map,
Range *  ents = NULL,
bool  id_in_tags = false,
bool  error_if_created = false 
)

get entities for TetGen data structure

Parameters
entsrange of entities (tetrahedrons or nodes)
intetgen data structure (look to TetGen user manual)
moab_tetgen_mapmapping MoAB to TetGen entities
tetgen_moab_mapmapping TetGen to moab entities
entsrerun entities which are in TetGen Data structure
id_in_tagsuse tags as entity handles, if that is a case use tag to find MoAB vertex id
error_if_createdthrow error if node need to be created

Definition at line 284 of file TetGenInterface.cpp.

288  {
290 
291  Interface &m_field = cOre;
292 
293  Tag th_marker;
294  int def_marker = 0;
295  CHKERR m_field.get_moab().tag_get_handle(
296  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
297  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
298 
299  int num_nodes = 0;
300  // std::vector<int> new_node_markes;
301  // new_node_markes.reserve(out.numberofpoints);
302  // std::vector<int> new_node_ii;
303  // new_node_ii.reserve(out.numberofpoints);
304 
305  int ii = 0;
306  for (; ii < out.numberofpoints; ii++) {
307  if (ii < in.numberofpoints) {
308  if (memcmp(&in.pointlist[3 * ii], &out.pointlist[3 * ii],
309  3 * sizeof(double)) == 0) {
310  unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
311  if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
312  CHKERR m_field.get_moab().tag_set_data(
313  th_marker, &tetgen_moab_map[iii], 1, &out.pointmarkerlist[ii]);
314  continue;
315  } else {
316  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
317  "data inconsistency between TetGen and MoAB");
318  }
319  }
320  }
321  if (id_in_tags) {
322  if (out.pointparamlist[ii].tag > 0) {
323  EntityHandle node;
324  CHKERR m_field.get_moab().handle_from_id(
325  MBVERTEX, in.pointparamlist[ii].tag - 1, node);
326  if (moab_tetgen_map.find(node) != moab_tetgen_map.end()) {
327  continue;
328  }
329  }
330  }
331  if (error_if_created) {
332  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
333  "node should not be created");
334  }
335  num_nodes++;
336  // new_node_ii.push_back(ii);
337  // new_node_markes.push_back(out.pointmarkerlist[ii]);
338  // EntityHandle node;
339  // rval = m_field.get_moab().create_vertex(&out.pointlist[3*ii],node);
340  // CHKERRQ_MOAB(rval);
341  // rval =
342  // m_field.get_moab().tag_set_data(th_marker,&node,1,&out.pointmarkerlist[ii]);
343  // CHKERRQ_MOAB(rval);
344  // moab_tetgen_map[node] = MBVERTEX|(ii<<MB_TYPE_WIDTH);
345  // tetgen_moab_map[MBVERTEX|(ii<<MB_TYPE_WIDTH)] = node;
346  // if(ents!=NULL) ents->insert(node);
347  }
348 
349  ReadUtilIface *iface;
350  CHKERR m_field.get_moab().query_interface(iface);
351 
352  if (num_nodes) {
353  vector<double *> arrays;
354  EntityHandle startv;
355  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
356  Range verts(startv, startv + num_nodes - 1);
357  int ii = in.numberofpoints;
358  for (Range::iterator vit = verts.begin(); vit != verts.end();
359  vit++, ii++) {
360  for (int nn = 0; nn != 3; nn++) {
361  arrays[nn][ii - in.numberofpoints] = out.pointlist[3 * ii + nn];
362  }
363  moab_tetgen_map[*vit] = MBVERTEX | (ii << MB_TYPE_WIDTH);
364  tetgen_moab_map[MBVERTEX | (ii << MB_TYPE_WIDTH)] = *vit;
365  if (ents != NULL)
366  ents->insert(*vit);
367  }
368  CHKERR m_field.get_moab().tag_set_data(
369  th_marker, verts, &out.pointmarkerlist[in.numberofpoints]);
370  }
371 
372  std::vector<int> tetgen_ii;
373 
374  // Build tets
375  std::vector<EntityHandle> conn_seq_tet;
376  conn_seq_tet.reserve(4 * out.numberoftetrahedra);
377  tetgen_ii.reserve(out.numberoftetrahedra);
378  conn_seq_tet.clear();
379  tetgen_ii.clear();
380  ii = 0;
381  for (; ii < out.numberoftetrahedra; ii++) {
382  unsigned long iii = MBTET | (ii << MB_TYPE_WIDTH);
383  if (ii < in.numberoftetrahedra) {
384  if (memcmp(&in.tetrahedronlist[4 * ii], &out.tetrahedronlist[4 * ii],
385  4 * sizeof(int)) == 0) {
386  if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
387  if (ents != NULL)
388  ents->insert(tetgen_moab_map[iii]);
389  continue;
390  } else {
391  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
392  "data inconsistency between TetGen and MoAB");
393  }
394  }
395  }
396  EntityHandle conn[4];
397  for (int nn = 0; nn < 4; nn++) {
398  int nnn = out.tetrahedronlist[4 * ii + nn];
399  if (tetgen_moab_map.find(MBVERTEX | (nnn << MB_TYPE_WIDTH)) ==
400  tetgen_moab_map.end()) {
401  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
402  "data inconsistency between TetGen and MoAB");
403  }
404  conn[nn] = tetgen_moab_map.at(MBVERTEX | (nnn << MB_TYPE_WIDTH));
405  }
406  auto get_coords = [&m_field](const EntityHandle *conn) {
407  VectorDouble12 coords(12);
408  CHKERR m_field.get_moab().get_coords(conn, 4, &coords[0]);
409  return coords;
410  };
411  auto get_volume = [](VectorDouble12 &&coords) {
412  return Tools::tetVolume(&coords[0]);
413  };
414  if (get_volume(get_coords(conn)) < 0) {
415  EntityHandle n0 = conn[0];
416  conn[0] = conn[1];
417  conn[1] = n0;
418  }
419  Range tets;
420  CHKERR m_field.get_moab().get_adjacencies(conn, 4, 3, false, tets);
421  bool tet_found = false;
422  for (auto tet : tets) {
423  const EntityHandle *tet_conn;
424  int num_nodes;
425  CHKERR m_field.get_moab().get_connectivity(tet, tet_conn, num_nodes,
426  true);
427  const EntityHandle *p = std::find(tet_conn, &tet_conn[4], conn[0]);
428  if (p != &tet_conn[4]) {
429  int s = std::distance(tet_conn, p);
430  int n = 0;
431  for (; n != 4; ++n) {
432  const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
433  if (tet_conn[idx[s + n]] != conn[n])
434  break;
435  }
436  if(n == 4 && !tet_found) {
437  moab_tetgen_map[tet] = iii;
438  tetgen_moab_map[iii] = tet;
439  if (ents != NULL)
440  ents->insert(tet);
441  tet_found = true;
442  } else if (n == 4) {
443  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
444  "More that one tet with the same connectivity");
445  }
446  } else {
447  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Imposible case");
448  }
449  }
450  if(!tet_found) {
451  for (int nn = 0; nn != 4; nn++) {
452  conn_seq_tet.push_back(conn[nn]);
453  }
454  tetgen_ii.push_back(ii);
455  }
456  }
457 
458  Range new_tets;
459  if (!conn_seq_tet.empty()) {
460  EntityHandle starte; // Connectivity
461  EntityHandle *conn;
462  int num_el = tetgen_ii.size();
463  CHKERR iface->get_element_connect(num_el, 4, MBTET, 0, starte, conn);
464  std::copy(conn_seq_tet.begin(), conn_seq_tet.end(), conn);
465  CHKERR iface->update_adjacencies(starte, num_el, 4, conn);
466  new_tets = Range(starte, starte + num_el - 1);
467  std::vector<int>::iterator ii_it = tetgen_ii.begin();
468  int ii = 0;
469  for (Range::iterator tit = new_tets.begin(); tit != new_tets.end();
470  tit++, ii_it++, ii++) {
471  unsigned long iii = MBTET | ((*ii_it) << MB_TYPE_WIDTH);
472  moab_tetgen_map[*tit] = iii;
473  tetgen_moab_map[iii] = *tit;
474  }
475  if (ents != NULL)
476  ents->merge(new_tets);
477  }
478 
480  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define MB_TYPE_WIDTH
Definition: definitions.h:285
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Common.hpp:223
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:47
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ outData() [2/2]

MoFEMErrorCode MoFEM::TetGenInterface::outData ( tetgenio &  in,
tetgenio &  out,
moabTetGen_Map moab_tetgen_map,
tetGenMoab_Map tetgen_moab_map,
BitRefLevel  bit,
bool  id_in_tags = false,
bool  error_if_created = false 
)

get entities for TetGen data structure

Parameters
entsrange of entities (tetrahedrons or nodes)
intetgen data structure (look to TetGen user manual)
moab_tetgen_mapmapping MoAB to TetGen entities
tetgen_moab_mapmapping TetGen to MoAB entities
entsrerun entities which are in TetGen data structure
bitset level to created entities
error_if_createdthrow error if node need to be created

Definition at line 481 of file TetGenInterface.cpp.

485  {
487  Interface &m_field = cOre;
488  Range ents;
489  CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
490  error_if_created);
491  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
492  ents.subset_by_type(MBTET), bit);
494  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define CHKERR
Inline error check.
Definition: definitions.h:578
MoFEMErrorCode outData(tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false)
get entities for TetGen data structure
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ query_interface()

MoFEMErrorCode MoFEM::TetGenInterface::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 47 of file TetGenInterface.cpp.

48  {
50  *iface = NULL;
51  if (uuid == IDD_MOFEMTetGegInterface) {
52  *iface = const_cast<TetGenInterface *>(this);
54  }
55  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
56 
58  }
TetGenInterface(const MoFEM::Core &core)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static const MOFEMuuid IDD_MOFEMTetGegInterface

◆ setFaceData()

MoFEMErrorCode MoFEM::TetGenInterface::setFaceData ( std::vector< std::pair< Range, int > > &  markers,
tetgenio &  in,
moabTetGen_Map moab_tetgen_map,
tetGenMoab_Map tetgen_moab_map 
)

set markers to faces

Parameters
markersdata structure with markers
intetgen data structure (look to TetGen user manual)
moab_tetgen_mapmapping MoAB to TetGen entities
tetgen_moab_mapmapping TetGen to MoAB entities

Definition at line 495 of file TetGenInterface.cpp.

498  {
500  ErrorCode rval;
501  Interface &m_field = cOre;
502  in.numberoffacets = markers.size();
503  in.facetlist = new tetgenio::facet[in.numberoffacets];
504  in.facetmarkerlist = new int[in.numberoffacets];
505  std::vector<std::pair<Range, int> >::iterator mit = markers.begin();
506  for (int ii = 0; mit != markers.end(); mit++, ii++) {
507  in.facetmarkerlist[ii] = mit->second;
508  Range &faces = mit->first;
509  tetgenio::facet *f = &(in.facetlist[ii]);
510  f->numberofpolygons = faces.size();
511  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
512  int jj = 0;
513  for (int dd = 3; dd >= 0; dd--) {
514  Range dd_faces = faces.subset_by_dimension(dd);
515  if (dd_faces.empty())
516  continue;
517  Range::iterator it = dd_faces.begin();
518  for (; it != dd_faces.end(); it++, jj++) {
519  int num_nodes;
520  const EntityHandle *conn;
521  tetgenio::polygon *p = &(f->polygonlist[jj]);
522  switch (m_field.get_moab().type_from_handle(*it)) {
523  case MBVERTEX: {
524  p->numberofvertices = 1;
525  conn = &*it;
526  } break;
527  default: {
528  rval =
529  m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
531  p->numberofvertices = num_nodes;
532  }
533  }
534  p->vertexlist = new int[p->numberofvertices];
535  for (int nn = 0; nn < p->numberofvertices; nn++) {
536  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
537  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
538  "data inconsistency between TetGen and MoAB");
539  }
540  p->vertexlist[nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
541  }
542  }
543  }
544  // holes
545  f->numberofholes = 0;
546  f->holelist = NULL;
547  }
549  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T *> &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MB_TYPE_WIDTH
Definition: definitions.h:285
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ setGeomData()

MoFEMErrorCode MoFEM::TetGenInterface::setGeomData ( tetgenio &  in,
moabTetGen_Map moab_tetgen_map,
tetGenMoab_Map tetgen_moab_map,
std::map< int, Range > &  type_ents 
)

set point tags and type

Set type of entity, look in TetGen manual for details

std::map<int,Range> types_ents;
//RIDGEVERTEX
types_ents[TetGenInterface::RIDGEVERTEX].merge(region_tets_skin_without_boundary_nodes);
//FREESEGVERTEX
types_ents[TetGenInterface::FREESEGVERTEX].merge(crack_surface_tris_skin_nodes);
//FREEFACETVERTEX
types_ents[TetGenInterface::FREEFACETVERTEX].merge(region_tets_skin_nodes);
//FREEVOLVERTEX
types_ents[TetGenInterface::FREEVOLVERTEX].merge(region_nodes);

Definition at line 236 of file TetGenInterface.cpp.

239  {
241 
242  Interface &m_field = cOre;
243  //
244  // ErrorCode rval;
245  in.pointparamlist = new tetgenio::pointparam[in.numberofpoints];
246  // std::vector<bool> points_is_set(in.numberofpoints,false);
247  std::map<int, Range>::iterator mit = type_ents.begin();
248  for (; mit != type_ents.end(); mit++) {
249  if (mit->first < 0 && mit->first > 3) {
250  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
251  "Wrong TetGen point type");
252  }
253  Range::iterator it = mit->second.begin();
254  for (; it != mit->second.end(); it++) {
255  moabTetGen_Map::iterator miit = moab_tetgen_map.find(*it);
256  if (miit == moab_tetgen_map.end()) {
257  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
258  "Data inconsistency between TetGen and MoAB");
259  continue;
260  }
261  int id = miit->second >> MB_TYPE_WIDTH;
262  in.pointparamlist[id].uv[0] = 0;
263  in.pointparamlist[id].uv[1] = 0;
264  in.pointparamlist[id].type = mit->first;
265  in.pointparamlist[id].tag = m_field.get_moab().id_from_handle(*it) + 1;
266  // points_is_set[id] = true;
267  }
268  }
269 
270  // // Check only if tag and type is set to all points
271  // for(
272  // std::vector<bool>::iterator bit = points_is_set.begin();
273  // bit!=points_is_set.end();bit++
274  // ) {
275  // if(!*bit) {
276  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Point type for
277  // TetGen is not set");
278  // }
279  // }
280 
282  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MB_TYPE_WIDTH
Definition: definitions.h:285
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ setRegionData()

MoFEMErrorCode MoFEM::TetGenInterface::setRegionData ( std::vector< std::pair< EntityHandle, int > > &  regions,
tetgenio &  in,
Tag  th = NULL 
)

set region data to tetrahedral

Definition at line 597 of file TetGenInterface.cpp.

599  {
601  ErrorCode rval;
602  Interface &m_field = cOre;
603  in.numberofregions = regions.size();
604  in.regionlist = new double[5 * in.numberofregions];
605  int kk = 0;
606  std::vector<std::pair<EntityHandle, int> >::iterator it = regions.begin();
607  for (int ii = 0; it != regions.end(); it++, ii++) {
608  double coords[3];
609  switch (m_field.get_moab().type_from_handle(it->first)) {
610  case MBVERTEX: {
611  if(th) {
612  rval = m_field.get_moab().tag_get_data(th, &it->first, 1, coords);
614  } else {
615  rval = m_field.get_moab().get_coords(&it->first, 1, coords);
617  }
618  } break;
619  case MBTET: {
620  int num_nodes;
621  const EntityHandle *conn;
622  rval = m_field.get_moab().get_connectivity(it->first, conn, num_nodes,
623  true);
625  double _coords[12];
626  if (th) {
627  rval = m_field.get_moab().tag_get_data(th, conn, num_nodes, _coords);
629  } else {
630  rval = m_field.get_moab().get_coords(conn, num_nodes, _coords);
632  }
633  coords[0] = (_coords[0] + _coords[3] + _coords[6] + _coords[9]) / 4.;
634  coords[1] = (_coords[1] + _coords[4] + _coords[7] + _coords[10]) / 4.;
635  coords[2] = (_coords[2] + _coords[5] + _coords[8] + _coords[11]) / 4.;
636  } break;
637  default:
638  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
639  }
640  for (int nn = 0; nn < 3; nn++) {
641  in.regionlist[kk++] = coords[nn];
642  }
643  in.regionlist[kk++] = it->second;
644  in.regionlist[kk++] = it->second;
645  }
647  }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ tetRahedralize()

MoFEMErrorCode MoFEM::TetGenInterface::tetRahedralize ( char  switches[],
tetgenio &  in,
tetgenio &  out 
)

run tetgen

Definition at line 691 of file TetGenInterface.cpp.

692  {
694  tetgenbehavior a;
695  a.parse_commandline(switches);
696  tetrahedralize(&a, &in, &out);
698  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

Member Data Documentation

◆ cOre

MoFEM::Core& MoFEM::TetGenInterface::cOre

Definition at line 41 of file TetGenInterface.hpp.


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