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

TetGen interface. 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, bool assume_first_nodes_the_same=false, Tag th=nullptr)
 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, bool assume_first_nodes_the_same=false, Tag th=nullptr)
 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

TetGen interface.

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 682 of file TetGenInterface.cpp.

684  {
686  double *pa = &coords[0];
687  double *pb = &coords[3];
688  double *pc = &coords[6];
689  double *pd = &coords[9];
690  double adx = pa[0] - pd[0];
691  double bdx = pb[0] - pd[0];
692  double cdx = pc[0] - pd[0];
693  double ady = pa[1] - pd[1];
694  double bdy = pb[1] - pd[1];
695  double cdy = pc[1] - pd[1];
696  double adz = pa[2] - pd[2];
697  double bdz = pb[2] - pd[2];
698  double cdz = pc[2] - pd[2];
699  double v = adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
700  cdx * (ady * bdz - adz * bdy);
701  double l = sqrt(pow(pa[0] - pb[0], 2) + pow(pa[1] - pb[1], 2) +
702  pow(pa[2] - pb[2], 2)) +
703  sqrt(pow(pa[0] - pc[0], 2) + pow(pa[1] - pc[1], 2) +
704  pow(pa[2] - pc[2], 2)) +
705  sqrt(pow(pa[0] - pd[0], 2) + pow(pa[1] - pd[1], 2) +
706  pow(pa[2] - pd[2], 2)) +
707  sqrt(pow(pb[0] - pc[0], 2) + pow(pb[1] - pc[1], 2) +
708  pow(pb[2] - pc[2], 2)) +
709  sqrt(pow(pb[0] - pd[0], 2) + pow(pb[1] - pd[1], 2) +
710  pow(pb[2] - pd[2], 2)) +
711  sqrt(pow(pc[0] - pd[0], 2) + pow(pc[1] - pd[1], 2) +
712  pow(pc[2] - pd[2], 2));
713  // std::cerr << fabs(v/pow(l,3)) << " ";
714  *result = fabs(v / pow(l, 3)) < eps ? true : false;
716 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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 626 of file TetGenInterface.cpp.

628  {
630  Interface &m_field = cOre;
631  int nbattributes = out.numberoftetrahedronattributes;
632  if (nbattributes == 0) {
633  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
634  "tetgen has no regions attributes");
635  }
636  Tag th_region;
637  rval = m_field.get_moab().tag_get_handle("TETGEN_REGION", th_region);
638  if (rval == MB_SUCCESS) {
639  rval = m_field.get_moab().tag_delete(th_region);
641  }
642  double def_marker = 0;
643  CHKERR m_field.get_moab().tag_get_handle(
644  "TETGEN_REGION", nbattributes, MB_TYPE_DOUBLE, th_region,
645  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
646  int ii = 0;
647  for (; ii < out.numberoftetrahedra; ii++) {
648  int jj = 0;
649  for (; jj < nbattributes; jj++) {
650  double id = out.tetrahedronattributelist[ii * nbattributes + jj];
651  int iii = MBTET | (ii << MB_TYPE_WIDTH);
652  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
653  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
654  "data inconsistency between TetGen and MoAB");
655  }
656  EntityHandle ent = tetgen_moab_map[iii];
657  CHKERR m_field.get_moab().tag_set_data(th_region, &ent, 1, &id);
658  if (ents != NULL)
659  ents->insert(ent);
660  if (ents_map != NULL)
661  (*ents_map)[id].insert(ent);
662  }
663  }
665 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:515
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MB_TYPE_WIDTH
Definition: definitions.h:288
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
#define CHKERR
Inline error check.
Definition: definitions.h:596

◆ 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 529 of file TetGenInterface.cpp.

531  {
533  ErrorCode rval;
534  Interface &m_field = cOre;
535  Tag th_marker;
536  int def_marker = 0;
537  rval = m_field.get_moab().tag_get_handle(
538  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
539  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
541  int ii = 0;
542  for (; ii < out.numberoftrifaces; ii++) {
543  if (only_non_zero) {
544  if (out.trifacemarkerlist[ii] == 0) {
545  continue;
546  }
547  }
548  EntityHandle conn[3];
549  for (int nn = 0; nn < 3; nn++) {
550  int iii = MBVERTEX | (out.trifacelist[3 * ii + nn] << MB_TYPE_WIDTH);
551  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
552  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
553  "data inconsistency between TetGen and MoAB");
554  } else {
555  conn[nn] = tetgen_moab_map[iii];
556  }
557  }
558  Range face;
559  rval = m_field.get_moab().get_adjacencies(conn, 3, 2, false, face);
561  face = face.subset_by_type(MBTRI);
562  if (face.size() != 1) {
563  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
564  "data inconsistency between TetGen and MoAB, %u", face.size());
565  }
566  if (ents != NULL)
567  ents->merge(face);
568  if (ents_map != NULL)
569  (*ents_map)[out.trifacemarkerlist[ii]].merge(face);
570  rval = m_field.get_moab().tag_set_data(th_marker, &*face.begin(), 1,
571  &out.trifacemarkerlist[ii]);
573  }
575 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:515
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MB_TYPE_WIDTH
Definition: definitions.h:288
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84

◆ groupPlanar_Triangle()

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

Definition at line 717 of file TetGenInterface.cpp.

719  {
721 
722  Interface &m_field = cOre;
723  Skinner skin(&m_field.get_moab());
724 
725  for (;;) {
726 
727  Range noplanar_to_anyother;
728  std::vector<Range>::iterator vit = sorted.begin();
729 
730  do {
731 
732  bool repeat = false;
733 
734  // get edges on vit skin
735  Range skin_edges;
736  CHKERR skin.find_skin(0, *vit, false, skin_edges);
737 
738  // skin edge nodes
739  Range skin_edges_nodes;
740  CHKERR m_field.get_moab().get_connectivity(skin_edges, skin_edges_nodes,
741  true);
742 
743  // get tris adjacent to vit skin edges
744  Range skin_edges_tris;
745  CHKERR m_field.get_moab().get_adjacencies(
746  skin_edges, 2, false, skin_edges_tris, moab::Interface::UNION);
747  // get tris which are part of facet
748  Range inner_tris = intersect(skin_edges_tris, *vit);
749  Range outer_tris = intersect(skin_edges_tris, tris);
750 
751  // tris coplanar with vit tris
752  Range coplanar_tris;
753 
754  Range::iterator tit = outer_tris.begin();
755  for (; tit != outer_tris.end(); tit++) {
756  Range tit_conn;
757  CHKERR m_field.get_moab().get_connectivity(&*tit, 1, tit_conn, true);
758  tit_conn = subtract(tit_conn, skin_edges_nodes);
759  if (tit_conn.empty()) {
760  coplanar_tris.insert(*tit);
761  noplanar_to_anyother.erase(*tit);
762  repeat = true;
763  } else {
764  Range tit_edges;
765  CHKERR m_field.get_moab().get_adjacencies(
766  &*tit, 1, 1, false, tit_edges, moab::Interface::UNION);
767  tit_edges = intersect(tit_edges, skin_edges);
768  if (tit_edges.size() != 1) {
769  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
770  "data inconsistency");
771  }
772  Range inner_tri;
773  CHKERR m_field.get_moab().get_adjacencies(
774  tit_edges, 2, false, inner_tri, moab::Interface::UNION);
775  inner_tri = intersect(inner_tri, inner_tris);
776  if (inner_tri.size() != 1) {
777  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
778  "data inconsistency");
779  }
780  // get connectivity
781  int num_nodes;
782  const EntityHandle *inner_tri_conn;
783  CHKERR m_field.get_moab().get_connectivity(
784  *inner_tri.begin(), inner_tri_conn, num_nodes, true);
785  double coords[12];
786  if (th) {
787  CHKERR m_field.get_moab().tag_get_data(th, inner_tri_conn, 3,
788  coords);
789  CHKERR m_field.get_moab().tag_get_data(th, &*tit_conn.begin(), 1,
790  &coords[9]);
791  } else {
792  CHKERR m_field.get_moab().get_coords(inner_tri_conn, 3, coords);
793  CHKERR m_field.get_moab().get_coords(&*tit_conn.begin(), 1,
794  &coords[9]);
795  }
796  bool coplanar;
797  CHKERR checkPlanar_Trinagle(coords, &coplanar, eps);
798  if (coplanar) {
799  coplanar_tris.insert(*tit);
800  noplanar_to_anyother.erase(*tit);
801  repeat = true;
802  } else {
803  noplanar_to_anyother.insert(*tit);
804  }
805  }
806  }
807 
808  vit->merge(coplanar_tris);
809  tris = subtract(tris, *vit);
810 
811  if (repeat) {
812  vit = sorted.begin();
813  } else {
814  vit++;
815  }
816 
817  } while (vit != sorted.end());
818 
819  if (noplanar_to_anyother.empty()) {
821  } else {
822  Range seed;
823  seed.insert(noplanar_to_anyother[0]);
824  tris.erase(noplanar_to_anyother[0]);
825  sorted.push_back(seed);
826  }
827  }
828 
830 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode checkPlanar_Trinagle(double coords[], bool *result, const double eps=1e-9)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
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 832 of file TetGenInterface.cpp.

833  {
835 
836  // PetscAttachDebugger();
837 
838  Range seed;
839  seed.insert(tris[0]);
840  tris.erase(tris[0]);
841  std::vector<Range> vec_seed;
842  vec_seed.push_back(seed);
843  sorted.push_back(vec_seed);
844 
845  for (;;) {
846  std::vector<Range> &vec = sorted.back();
847  CHKERR groupPlanar_Triangle(tris, vec, eps);
848  if (tris.empty()) {
850  } else {
851  Range seed;
852  seed.insert(tris[0]);
853  tris.erase(tris[0]);
854  std::vector<Range> vec_seed;
855  vec_seed.push_back(seed);
856  sorted.push_back(vec_seed);
857  }
858  }
859 
861 }
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:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define CHKERR
Inline error check.
Definition: definitions.h:596
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 58 of file TetGenInterface.cpp.

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

◆ loadPoly()

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

load poly file

Definition at line 674 of file TetGenInterface.cpp.

674  {
676  if (!in.load_poly(file_name)) {
677  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
678  "can not read TetGen poly file");
679  }
681 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ 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 862 of file TetGenInterface.cpp.

865  {
867  // FIXME: assumes that are no holes
868 
869  if (ents.empty()) {
870  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
871  "no ents to build polygon");
872  }
873 
874  Interface &m_field = cOre;
875 
876  Skinner skin(&m_field.get_moab());
877 
878  Range skin_edges;
879  CHKERR skin.find_skin(0, ents, false, skin_edges);
880 
881  std::vector<EntityHandle> polygon_nodes;
882  EntityHandle seed = skin_edges[0];
883  Range seen_edges;
884  seen_edges.insert(seed);
885  skin_edges.erase(seed);
886  int num_nodes;
887  const EntityHandle *conn;
888  CHKERR m_field.get_moab().get_connectivity(seed, conn, num_nodes, true);
889  polygon_nodes.push_back(conn[0]);
890  // std::cerr << std::endl;
891  // std::cerr << conn[0] << " " << conn[1] << std::endl;
892  do {
893  EntityHandle last_node = polygon_nodes.back();
894  Range adj_edges;
895  CHKERR m_field.get_moab().get_adjacencies(&last_node, 1, 1, false,
896  adj_edges);
897  adj_edges = intersect(adj_edges, skin_edges);
898  if (adj_edges.size() == 0) {
899  break;
900  }
901  if (adj_edges.size() != 1) {
902  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
903  "should be only one edge");
904  }
905  seen_edges.insert(adj_edges[0]);
906  skin_edges.erase(adj_edges[0]);
907  CHKERR m_field.get_moab().get_connectivity(adj_edges[0], conn, num_nodes,
908  true);
909  EntityHandle add_node = (last_node == conn[0]) ? conn[1] : conn[0];
910  polygon_nodes.push_back(add_node);
911  // std::cerr << "\t" << add_node << std::endl;
912  } while (1);
913 
914  if (reduce_edges) {
915  // std::cerr << "polygon " << polygon_nodes.size();
916  std::vector<EntityHandle>::iterator pit = polygon_nodes.begin();
917  // std::cerr << std::endl;
918  for (; pit != polygon_nodes.end();) {
919  if (not_reducable_nodes != NULL) {
920  if (not_reducable_nodes->find(*pit) != not_reducable_nodes->end()) {
921  pit++;
922  continue;
923  }
924  }
925  EntityHandle mm;
926  if (pit == polygon_nodes.begin()) {
927  mm = polygon_nodes.back();
928  } else {
929  mm = *(pit - 1);
930  }
931  EntityHandle mc = *pit;
932  EntityHandle mp;
933  if (polygon_nodes.back() == mc) {
934  mp = polygon_nodes[0];
935  } else {
936  mp = *(pit + 1);
937  }
938  double coords[9];
939  if (th) {
940  CHKERR m_field.get_moab().tag_get_data(th, &mm, 1, &coords[3 * 0]);
941  CHKERR m_field.get_moab().tag_get_data(th, &mc, 1, &coords[3 * 1]);
942  CHKERR m_field.get_moab().tag_get_data(th, &mp, 1, &coords[3 * 2]);
943  } else {
944  CHKERR m_field.get_moab().get_coords(&mm, 1, &coords[3 * 0]);
945  CHKERR m_field.get_moab().get_coords(&mc, 1, &coords[3 * 1]);
946  CHKERR m_field.get_moab().get_coords(&mp, 1, &coords[3 * 2]);
947  }
948  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 0], 1); // mm = mm -
949  // mc
950  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 2], 1); // mp = mp -
951  // mc
952  double spin[9];
953  CHKERR Spin(spin, &coords[3 * 0]);
954  double l0 = cblas_dnrm2(3, &coords[3 * 0], 1);
955  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1. / l0, spin, 3,
956  &coords[3 * 2], 1, 0., &coords[3 * 1], 1);
957  double dot = cblas_dnrm2(3, &coords[3 * 1], 1);
958  // std::cerr << mm << " " << mc << " " << mp << " " << dot << std::endl;
959  if (dot < eps) {
960  polygon_nodes.erase(pit);
961  pit = polygon_nodes.begin();
962  // std::cerr << std::endl;
963  } else {
964  pit++;
965  }
966  }
967  }
968 
969  Range existing_polygon;
970  CHKERR m_field.get_moab().get_adjacencies(
971  &polygon_nodes[0], polygon_nodes.size(), 2, true, existing_polygon);
972  if (existing_polygon.empty()) {
973  EntityHandle polygon;
974  CHKERR m_field.get_moab().create_element(MBPOLYGON, &polygon_nodes[0],
975  polygon_nodes.size(), polygon);
976  polygons.insert(polygon);
977  } else {
978  polygons.merge(existing_polygon);
979  }
980 
982 }
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:558
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 MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
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,
bool  assume_first_nodes_the_same = false,
Tag  th = nullptr 
)

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 266 of file TetGenInterface.cpp.

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

◆ 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,
bool  assume_first_nodes_the_same = false,
Tag  th = nullptr 
)

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 459 of file TetGenInterface.cpp.

464  {
466  Interface &m_field = cOre;
467  Range ents;
468  CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
469  error_if_created, assume_first_nodes_the_same,th);
470  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
471  ents.subset_by_type(MBTET), bit);
473 }
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, bool assume_first_nodes_the_same=false, Tag th=nullptr)
get entities for TetGen data structure
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 44 of file TetGenInterface.cpp.

45  {
47  *iface = NULL;
48  if (uuid == IDD_MOFEMTetGegInterface) {
49  *iface = const_cast<TetGenInterface *>(this);
51  }
52  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
53 
55 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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 474 of file TetGenInterface.cpp.

477  {
479  ErrorCode rval;
480  Interface &m_field = cOre;
481  in.numberoffacets = markers.size();
482  in.facetlist = new tetgenio::facet[in.numberoffacets];
483  in.facetmarkerlist = new int[in.numberoffacets];
484  std::vector<std::pair<Range, int>>::iterator mit = markers.begin();
485  for (int ii = 0; mit != markers.end(); mit++, ii++) {
486  in.facetmarkerlist[ii] = mit->second;
487  Range &faces = mit->first;
488  tetgenio::facet *f = &(in.facetlist[ii]);
489  f->numberofpolygons = faces.size();
490  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
491  int jj = 0;
492  for (int dd = 3; dd >= 0; dd--) {
493  Range dd_faces = faces.subset_by_dimension(dd);
494  if (dd_faces.empty())
495  continue;
496  Range::iterator it = dd_faces.begin();
497  for (; it != dd_faces.end(); it++, jj++) {
498  int num_nodes;
499  const EntityHandle *conn;
500  tetgenio::polygon *p = &(f->polygonlist[jj]);
501  switch (m_field.get_moab().type_from_handle(*it)) {
502  case MBVERTEX: {
503  p->numberofvertices = 1;
504  conn = &*it;
505  } break;
506  default: {
507  rval =
508  m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
510  p->numberofvertices = num_nodes;
511  }
512  }
513  p->vertexlist = new int[p->numberofvertices];
514  for (int nn = 0; nn < p->numberofvertices; nn++) {
515  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
516  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
517  "data inconsistency between TetGen and MoAB");
518  }
519  p->vertexlist[nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
520  }
521  }
522  }
523  // holes
524  f->numberofholes = 0;
525  f->holelist = NULL;
526  }
528 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:515
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MB_TYPE_WIDTH
Definition: definitions.h:288
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
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

◆ 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 217 of file TetGenInterface.cpp.

220  {
222 
223  Interface &m_field = cOre;
224  //
225  // ErrorCode rval;
226  in.pointparamlist = new tetgenio::pointparam[in.numberofpoints];
227  // std::vector<bool> points_is_set(in.numberofpoints,false);
228  std::map<int, Range>::iterator mit = type_ents.begin();
229  for (; mit != type_ents.end(); mit++) {
230  if (mit->first < 0 && mit->first > 3) {
231  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
232  "Wrong TetGen point type");
233  }
234  Range::iterator it = mit->second.begin();
235  for (; it != mit->second.end(); it++) {
236  moabTetGen_Map::iterator miit = moab_tetgen_map.find(*it);
237  if (miit == moab_tetgen_map.end()) {
238  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
239  "Data inconsistency between TetGen and MoAB");
240  continue;
241  }
242  int id = miit->second >> MB_TYPE_WIDTH;
243  in.pointparamlist[id].uv[0] = 0;
244  in.pointparamlist[id].uv[1] = 0;
245  in.pointparamlist[id].type = mit->first;
246  in.pointparamlist[id].tag = m_field.get_moab().id_from_handle(*it) + 1;
247  // points_is_set[id] = true;
248  }
249  }
250 
251  // // Check only if tag and type is set to all points
252  // for(
253  // std::vector<bool>::iterator bit = points_is_set.begin();
254  // bit!=points_is_set.end();bit++
255  // ) {
256  // if(!*bit) {
257  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Point type for
258  // TetGen is not set");
259  // }
260  // }
261 
263 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MB_TYPE_WIDTH
Definition: definitions.h:288
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ 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 576 of file TetGenInterface.cpp.

577  {
579  ErrorCode rval;
580  Interface &m_field = cOre;
581  in.numberofregions = regions.size();
582  in.regionlist = new double[5 * in.numberofregions];
583  int kk = 0;
584  std::vector<std::pair<EntityHandle, int>>::iterator it = regions.begin();
585  for (int ii = 0; it != regions.end(); it++, ii++) {
586  double coords[3];
587  switch (m_field.get_moab().type_from_handle(it->first)) {
588  case MBVERTEX: {
589  if (th) {
590  rval = m_field.get_moab().tag_get_data(th, &it->first, 1, coords);
592  } else {
593  rval = m_field.get_moab().get_coords(&it->first, 1, coords);
595  }
596  } break;
597  case MBTET: {
598  int num_nodes;
599  const EntityHandle *conn;
600  rval =
601  m_field.get_moab().get_connectivity(it->first, conn, num_nodes, true);
603  double _coords[12];
604  if (th) {
605  rval = m_field.get_moab().tag_get_data(th, conn, num_nodes, _coords);
607  } else {
608  rval = m_field.get_moab().get_coords(conn, num_nodes, _coords);
610  }
611  coords[0] = (_coords[0] + _coords[3] + _coords[6] + _coords[9]) / 4.;
612  coords[1] = (_coords[1] + _coords[4] + _coords[7] + _coords[10]) / 4.;
613  coords[2] = (_coords[2] + _coords[5] + _coords[8] + _coords[11]) / 4.;
614  } break;
615  default:
616  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
617  }
618  for (int nn = 0; nn < 3; nn++) {
619  in.regionlist[kk++] = coords[nn];
620  }
621  in.regionlist[kk++] = it->second;
622  in.regionlist[kk++] = it->second;
623  }
625 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:515
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84

◆ tetRahedralize()

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

run tetgen

Definition at line 666 of file TetGenInterface.cpp.

667  {
669  tetgenbehavior a;
670  a.parse_commandline(switches);
671  tetrahedralize(&a, &in, &out);
673 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

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: