v0.8.20
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, 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

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

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

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

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

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

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

838  {
840 
841  // PetscAttachDebugger();
842 
843  Range seed;
844  seed.insert(tris[0]);
845  tris.erase(tris[0]);
846  std::vector<Range> vec_seed;
847  vec_seed.push_back(seed);
848  sorted.push_back(vec_seed);
849 
850  for (;;) {
851  std::vector<Range> &vec = sorted.back();
852  CHKERR groupPlanar_Triangle(tris, vec, eps);
853  if (tris.empty()) {
855  } else {
856  Range seed;
857  seed.insert(tris[0]);
858  tris.erase(tris[0]);
859  std::vector<Range> vec_seed;
860  vec_seed.push_back(seed);
861  sorted.push_back(vec_seed);
862  }
863  }
864 
866 }
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:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
#define CHKERR
Inline error check.
Definition: definitions.h:594
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  //
68  ErrorCode rval;
69 
70  Tag th_marker;
71  int def_marker = 0;
72  CHKERR m_field.get_moab().tag_get_handle(
73  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
74  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
75 
76  // All indices start from 0
77  in.firstnumber = 0;
78 
79  Range points = ents.subset_by_dimension(0);
80  in.numberofpoints = points.size();
81  if (points.size() > 0) {
82  in.pointlist = new double[in.numberofpoints * 3];
83  in.pointmarkerlist = new int[in.numberofpoints];
84  if (th) {
85  CHKERR m_field.get_moab().tag_get_data(th, points, in.pointlist);
86  } else {
87  CHKERR m_field.get_moab().get_coords(points, in.pointlist);
88  }
89  CHKERR m_field.get_moab().tag_get_data(th_marker, points,
90  in.pointmarkerlist);
91  it = points.begin();
92  for (int ii = 0; it != points.end(); it++, ii++) {
93  unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
94  tetgen_moab_map[iii] = *it;
95  moab_tetgen_map[*it] = iii;
96  }
97  }
98 
99  Range tets = ents.subset_by_type(MBTET);
100  in.numberoftetrahedra = tets.size();
101  if (in.numberoftetrahedra > 0) {
102  in.tetrahedronlist = new int[4 * ents.subset_by_type(MBTET).size()];
103  it = tets.begin();
104  for (int ii = 0; it != tets.end(); it++, ii++) {
105  int num_nodes;
106  const EntityHandle *conn;
107  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
108 #ifdef DEBUG_TETGEN
109  {
110  double coords[12];
111  if (th) {
112  CHKERR m_field.tag_get_data(th, conn, num_nodes, coords);
113  } else {
114  CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords);
115  }
116  double diff_n[12];
117  ShapeDiffMBTET(diff_n);
118  FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1],
119  &diff_n[2], 3);
120  FTensor::Tensor1<double *, 3> t_coords(&coords[0], &coords[1],
121  &coords[2], 3);
125  jac(i, j) = 0;
126  for (int nn = 0; nn != 4; nn++) {
127  jac(i, j) += t_coords(i) * t_diff_n(j);
128  ++t_coords;
129  ++t_diff_n;
130  }
131  double det;
132  determinantTensor3by3(jac, det);
133  if (det <= 0) {
134  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
135  "Negative volume det = %6.4e", det);
136  }
137  }
138 #endif
139  tetgen_moab_map[MBTET | (ii << MB_TYPE_WIDTH)] = *it;
140  moab_tetgen_map[*it] = MBTET | (ii << MB_TYPE_WIDTH);
141  for (int nn = 0; nn != 4; nn++) {
142  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
143  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144  "data inconsistency between TetGen and MoAB");
145  }
146  in.tetrahedronlist[4 * ii + nn] =
147  moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
148  }
149  }
150  }
151 
152  Range tris = ents.subset_by_type(MBTRI);
153  in.numberoftrifaces = tris.size();
154  if (in.numberoftrifaces) {
155  in.trifacelist = new int[3 * in.numberoftrifaces];
156  in.trifacemarkerlist = new int[in.numberoftrifaces];
157  // std::fill(&in.trifacemarkerlist[0],&in.trifacemarkerlist[in.numberoftrifaces],1);
158  CHKERR m_field.get_moab().tag_get_data(th_marker, tris,
159  in.trifacemarkerlist);
160  it = tris.begin();
161  for (int ii = 0; it != tris.end(); it++, ii++) {
162  int num_nodes;
163  const EntityHandle *conn;
164  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
165  int order[] = {0, 1, 2};
166  Range adj_tets;
167  CHKERR m_field.get_moab().get_adjacencies(&*it, 1, 3, true, adj_tets);
168  adj_tets = intersect(adj_tets, tets);
169  if (adj_tets.size() == 1) {
170  int side_number;
171  int sense;
172  int offset;
173  CHKERR m_field.get_moab().side_number(adj_tets[0], *it, side_number,
174  sense, offset);
175  if (sense == -1) {
176  order[0] = 1;
177  order[1] = 0;
178  }
179  }
180  tetgen_moab_map[MBTRI | (ii << MB_TYPE_WIDTH)] = *it;
181  moab_tetgen_map[*it] = MBTRI | (ii << MB_TYPE_WIDTH);
182  for (int nn = 0; nn < 3; nn++) {
183  if (moab_tetgen_map.find(conn[order[nn]]) == moab_tetgen_map.end()) {
184  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
185  "data inconsistency between TetGen and MoAB");
186  }
187  in.trifacelist[3 * ii + nn] =
188  moab_tetgen_map[conn[order[nn]]] >> MB_TYPE_WIDTH;
189  }
190  }
191  }
192 
193  Range edges = ents.subset_by_type(MBEDGE);
194  in.numberofedges = edges.size();
195  if (in.numberofedges > 0) {
196  in.edgelist = new int[2 * in.numberofedges];
197  in.edgemarkerlist = new int[in.numberofedges];
198  // std::fill(&in.edgemarkerlist[0],&in.edgemarkerlist[in.numberofedges],1);
199  CHKERR m_field.get_moab().tag_get_data(th_marker, edges, in.edgemarkerlist);
200  it = edges.begin();
201  for (int ii = 0; it != edges.end(); it++, ii++) {
202  int num_nodes;
203  const EntityHandle *conn;
204  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
205  tetgen_moab_map[MBEDGE | (ii << MB_TYPE_WIDTH)] = *it;
206  moab_tetgen_map[*it] = MBEDGE | (ii << MB_TYPE_WIDTH);
207  for (int nn = 0; nn < 2; nn++) {
208  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
209  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
210  "data inconsistency between TetGen and MoAB");
211  }
212  in.edgelist[2 * ii + nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
213  }
214  }
215  }
216 
218 }
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define MB_TYPE_WIDTH
Definition: definitions.h:286
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:303
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ loadPoly()

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

load poly file

Definition at line 677 of file TetGenInterface.cpp.

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

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

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

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

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

467  {
469  Interface &m_field = cOre;
470  Range ents;
471  CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
472  error_if_created, assume_first_nodes_the_same,th);
473  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
474  ents.subset_by_type(MBTET), bit);
476 }
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:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ 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 }
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:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
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 477 of file TetGenInterface.cpp.

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

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

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

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

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

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

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: