v0.13.2
Loading...
Searching...
No Matches
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, RangeidxRange_Map
 

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, 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
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register 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 ()=default
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

TetGen interface.

Definition at line 23 of file TetGenInterface.hpp.

Member Typedef Documentation

◆ idxRange_Map

Definition at line 34 of file TetGenInterface.hpp.

◆ moabTetGen_Map

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

Definition at line 32 of file TetGenInterface.hpp.

◆ tetGenMoab_Map

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

Definition at line 33 of file TetGenInterface.hpp.

Member Enumeration Documentation

◆ tetGenNodesTypes

Enumerator
RIDGEVERTEX 
FREESEGVERTEX 
FREEFACETVERTEX 
FREEVOLVERTEX 

Definition at line 49 of file TetGenInterface.hpp.

Constructor & Destructor Documentation

◆ TetGenInterface()

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

Definition at line 29 of file TetGenInterface.hpp.

30 : cOre(const_cast<MoFEM::Core &>(core)) {}
Core (interface) class.
Definition: Core.hpp:82

Member Function Documentation

◆ checkPlanar_Trinagle()

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

Definition at line 662 of file TetGenInterface.cpp.

664 {
666 double *pa = &coords[0];
667 double *pb = &coords[3];
668 double *pc = &coords[6];
669 double *pd = &coords[9];
670 double adx = pa[0] - pd[0];
671 double bdx = pb[0] - pd[0];
672 double cdx = pc[0] - pd[0];
673 double ady = pa[1] - pd[1];
674 double bdy = pb[1] - pd[1];
675 double cdy = pc[1] - pd[1];
676 double adz = pa[2] - pd[2];
677 double bdz = pb[2] - pd[2];
678 double cdz = pc[2] - pd[2];
679 double v = adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
680 cdx * (ady * bdz - adz * bdy);
681 double l = sqrt(pow(pa[0] - pb[0], 2) + pow(pa[1] - pb[1], 2) +
682 pow(pa[2] - pb[2], 2)) +
683 sqrt(pow(pa[0] - pc[0], 2) + pow(pa[1] - pc[1], 2) +
684 pow(pa[2] - pc[2], 2)) +
685 sqrt(pow(pa[0] - pd[0], 2) + pow(pa[1] - pd[1], 2) +
686 pow(pa[2] - pd[2], 2)) +
687 sqrt(pow(pb[0] - pc[0], 2) + pow(pb[1] - pc[1], 2) +
688 pow(pb[2] - pc[2], 2)) +
689 sqrt(pow(pb[0] - pd[0], 2) + pow(pb[1] - pd[1], 2) +
690 pow(pb[2] - pd[2], 2)) +
691 sqrt(pow(pc[0] - pd[0], 2) + pow(pc[1] - pd[1], 2) +
692 pow(pc[2] - pd[2], 2));
693 // std::cerr << fabs(v/pow(l,3)) << " ";
694 *result = fabs(v / pow(l, 3)) < eps ? true : false;
696}
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l

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

608 {
610 Interface &m_field = cOre;
611 int nbattributes = out.numberoftetrahedronattributes;
612 if (nbattributes == 0) {
613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614 "tetgen has no regions attributes");
615 }
616 Tag th_region;
617 rval = m_field.get_moab().tag_get_handle("TETGEN_REGION", th_region);
618 if (rval == MB_SUCCESS) {
619 rval = m_field.get_moab().tag_delete(th_region);
621 }
622 double def_marker = 0;
623 CHKERR m_field.get_moab().tag_get_handle(
624 "TETGEN_REGION", nbattributes, MB_TYPE_DOUBLE, th_region,
625 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
626 int ii = 0;
627 for (; ii < out.numberoftetrahedra; ii++) {
628 int jj = 0;
629 for (; jj < nbattributes; jj++) {
630 double id = out.tetrahedronattributelist[ii * nbattributes + jj];
631 int iii = MBTET | (ii << MB_TYPE_WIDTH);
632 if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
633 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
634 "data inconsistency between TetGen and MoAB");
635 }
636 EntityHandle ent = tetgen_moab_map[iii];
637 CHKERR m_field.get_moab().tag_set_data(th_region, &ent, 1, &id);
638 if (ents != NULL)
639 ents->insert(ent);
640 if (ents_map != NULL)
641 (*ents_map)[id].insert(ent);
642 }
643 }
645}
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MB_TYPE_WIDTH
Definition: definitions.h:226
#define CHKERR
Inline error check.
Definition: definitions.h:535
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975

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

511 {
513 ErrorCode rval;
514 Interface &m_field = cOre;
515 Tag th_marker;
516 int def_marker = 0;
517 rval = m_field.get_moab().tag_get_handle(
518 "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
519 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
521 int ii = 0;
522 for (; ii < out.numberoftrifaces; ii++) {
523 if (only_non_zero) {
524 if (out.trifacemarkerlist[ii] == 0) {
525 continue;
526 }
527 }
528 EntityHandle conn[3];
529 for (int nn = 0; nn < 3; nn++) {
530 int iii = MBVERTEX | (out.trifacelist[3 * ii + nn] << MB_TYPE_WIDTH);
531 if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
532 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533 "data inconsistency between TetGen and MoAB");
534 } else {
535 conn[nn] = tetgen_moab_map[iii];
536 }
537 }
538 Range face;
539 rval = m_field.get_moab().get_adjacencies(conn, 3, 2, false, face);
541 face = face.subset_by_type(MBTRI);
542 if (face.size() != 1) {
543 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
544 "data inconsistency between TetGen and MoAB, %u", face.size());
545 }
546 if (ents != NULL)
547 ents->merge(face);
548 if (ents_map != NULL)
549 (*ents_map)[out.trifacemarkerlist[ii]].merge(face);
550 rval = m_field.get_moab().tag_set_data(th_marker, &*face.begin(), 1,
551 &out.trifacemarkerlist[ii]);
553 }
555}

◆ groupPlanar_Triangle()

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

Definition at line 697 of file TetGenInterface.cpp.

699 {
701
702 Interface &m_field = cOre;
703 Skinner skin(&m_field.get_moab());
704
705 for (;;) {
706
707 Range noplanar_to_anyother;
708 std::vector<Range>::iterator vit = sorted.begin();
709
710 do {
711
712 bool repeat = false;
713
714 // get edges on vit skin
715 Range skin_edges;
716 CHKERR skin.find_skin(0, *vit, false, skin_edges);
717
718 // skin edge nodes
719 Range skin_edges_nodes;
720 CHKERR m_field.get_moab().get_connectivity(skin_edges, skin_edges_nodes,
721 true);
722
723 // get tris adjacent to vit skin edges
724 Range skin_edges_tris;
725 CHKERR m_field.get_moab().get_adjacencies(
726 skin_edges, 2, false, skin_edges_tris, moab::Interface::UNION);
727 // get tris which are part of facet
728 Range inner_tris = intersect(skin_edges_tris, *vit);
729 Range outer_tris = intersect(skin_edges_tris, tris);
730
731 // tris coplanar with vit tris
732 Range coplanar_tris;
733
734 Range::iterator tit = outer_tris.begin();
735 for (; tit != outer_tris.end(); tit++) {
736 Range tit_conn;
737 CHKERR m_field.get_moab().get_connectivity(&*tit, 1, tit_conn, true);
738 tit_conn = subtract(tit_conn, skin_edges_nodes);
739 if (tit_conn.empty()) {
740 coplanar_tris.insert(*tit);
741 noplanar_to_anyother.erase(*tit);
742 repeat = true;
743 } else {
744 Range tit_edges;
745 CHKERR m_field.get_moab().get_adjacencies(
746 &*tit, 1, 1, false, tit_edges, moab::Interface::UNION);
747 tit_edges = intersect(tit_edges, skin_edges);
748 if (tit_edges.size() != 1) {
749 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
750 "data inconsistency");
751 }
752 Range inner_tri;
753 CHKERR m_field.get_moab().get_adjacencies(
754 tit_edges, 2, false, inner_tri, moab::Interface::UNION);
755 inner_tri = intersect(inner_tri, inner_tris);
756 if (inner_tri.size() != 1) {
757 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
758 "data inconsistency");
759 }
760 // get connectivity
761 int num_nodes;
762 const EntityHandle *inner_tri_conn;
763 CHKERR m_field.get_moab().get_connectivity(
764 *inner_tri.begin(), inner_tri_conn, num_nodes, true);
765 double coords[12];
766 if (th) {
767 CHKERR m_field.get_moab().tag_get_data(th, inner_tri_conn, 3,
768 coords);
769 CHKERR m_field.get_moab().tag_get_data(th, &*tit_conn.begin(), 1,
770 &coords[9]);
771 } else {
772 CHKERR m_field.get_moab().get_coords(inner_tri_conn, 3, coords);
773 CHKERR m_field.get_moab().get_coords(&*tit_conn.begin(), 1,
774 &coords[9]);
775 }
776 bool coplanar;
777 CHKERR checkPlanar_Trinagle(coords, &coplanar, eps);
778 if (coplanar) {
779 coplanar_tris.insert(*tit);
780 noplanar_to_anyother.erase(*tit);
781 repeat = true;
782 } else {
783 noplanar_to_anyother.insert(*tit);
784 }
785 }
786 }
787
788 vit->merge(coplanar_tris);
789 tris = subtract(tris, *vit);
790
791 if (repeat) {
792 vit = sorted.begin();
793 } else {
794 vit++;
795 }
796
797 } while (vit != sorted.end());
798
799 if (noplanar_to_anyother.empty()) {
801 } else {
802 Range seed;
803 seed.insert(noplanar_to_anyother[0]);
804 tris.erase(noplanar_to_anyother[0]);
805 sorted.push_back(seed);
806 }
807 }
808
810}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMErrorCode checkPlanar_Trinagle(double coords[], bool *result, const double eps=1e-9)

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

813 {
815
816 // PetscAttachDebugger();
817
818 Range seed;
819 seed.insert(tris[0]);
820 tris.erase(tris[0]);
821 std::vector<Range> vec_seed;
822 vec_seed.push_back(seed);
823 sorted.push_back(vec_seed);
824
825 for (;;) {
826 std::vector<Range> &vec = sorted.back();
827 CHKERR groupPlanar_Triangle(tris, vec, eps);
828 if (tris.empty()) {
830 } else {
831 Range seed;
832 seed.insert(tris[0]);
833 tris.erase(tris[0]);
834 std::vector<Range> vec_seed;
835 vec_seed.push_back(seed);
836 sorted.push_back(vec_seed);
837 }
838 }
839
841}
MoFEMErrorCode groupPlanar_Triangle(Range &tris, std::vector< Range > &sorted, const double eps=1e-9, Tag th=NULL)

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

41 {
43
44 Interface &m_field = cOre;
45 Range::iterator it;
46
47 Tag th_marker;
48 int def_marker = 0;
49 CHKERR m_field.get_moab().tag_get_handle(
50 "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
51 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
52
53 // All indices start from 0
54 in.firstnumber = 0;
55
56 Range points = ents.subset_by_dimension(0);
57 in.numberofpoints = points.size();
58 if (points.size() > 0) {
59 in.pointlist = new double[in.numberofpoints * 3];
60 in.pointmarkerlist = new int[in.numberofpoints];
61 if (th) {
62 CHKERR m_field.get_moab().tag_get_data(th, points, in.pointlist);
63 } else {
64 CHKERR m_field.get_moab().get_coords(points, in.pointlist);
65 }
66 CHKERR m_field.get_moab().tag_get_data(th_marker, points,
67 in.pointmarkerlist);
68 it = points.begin();
69 for (int ii = 0; it != points.end(); it++, ii++) {
70 unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
71 tetgen_moab_map[iii] = *it;
72 moab_tetgen_map[*it] = iii;
73 }
74 }
75
76 Range tets = ents.subset_by_type(MBTET);
77 in.numberoftetrahedra = tets.size();
78 if (in.numberoftetrahedra > 0) {
79 in.tetrahedronlist = new int[4 * ents.subset_by_type(MBTET).size()];
80 it = tets.begin();
81 for (int ii = 0; it != tets.end(); it++, ii++) {
82 int num_nodes;
83 const EntityHandle *conn;
84 CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
85#ifdef DEBUG_TETGEN
86 {
87 double coords[12];
88 if (th) {
89 CHKERR m_field.tag_get_data(th, conn, num_nodes, coords);
90 } else {
91 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords);
92 }
93 double diff_n[12];
94 ShapeDiffMBTET(diff_n);
95 FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1],
96 &diff_n[2], 3);
97 FTensor::Tensor1<double *, 3> t_coords(&coords[0], &coords[1],
98 &coords[2], 3);
102 jac(i, j) = 0;
103 for (int nn = 0; nn != 4; nn++) {
104 jac(i, j) += t_coords(i) * t_diff_n(j);
105 ++t_coords;
106 ++t_diff_n;
107 }
108 double det;
109 determinantTensor3by3(jac, det);
110 if (det <= 0) {
111 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
112 "Negative volume det = %6.4e", det);
113 }
114 }
115#endif
116 tetgen_moab_map[MBTET | (ii << MB_TYPE_WIDTH)] = *it;
117 moab_tetgen_map[*it] = MBTET | (ii << MB_TYPE_WIDTH);
118 for (int nn = 0; nn != 4; nn++) {
119 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
120 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
121 "data inconsistency between TetGen and MoAB");
122 }
123 in.tetrahedronlist[4 * ii + nn] =
124 moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
125 }
126 }
127 }
128
129 Range tris = ents.subset_by_type(MBTRI);
130 in.numberoftrifaces = tris.size();
131 if (in.numberoftrifaces) {
132 in.trifacelist = new int[3 * in.numberoftrifaces];
133 in.trifacemarkerlist = new int[in.numberoftrifaces];
134 // std::fill(&in.trifacemarkerlist[0],&in.trifacemarkerlist[in.numberoftrifaces],1);
135 CHKERR m_field.get_moab().tag_get_data(th_marker, tris,
136 in.trifacemarkerlist);
137 it = tris.begin();
138 for (int ii = 0; it != tris.end(); it++, ii++) {
139 int num_nodes;
140 const EntityHandle *conn;
141 CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
142 int order[] = {0, 1, 2};
143 Range adj_tets;
144 CHKERR m_field.get_moab().get_adjacencies(&*it, 1, 3, true, adj_tets);
145 adj_tets = intersect(adj_tets, tets);
146 if (adj_tets.size() == 1) {
147 int side_number;
148 int sense;
149 int offset;
150 CHKERR m_field.get_moab().side_number(adj_tets[0], *it, side_number,
151 sense, offset);
152 if (sense == -1) {
153 order[0] = 1;
154 order[1] = 0;
155 }
156 }
157 tetgen_moab_map[MBTRI | (ii << MB_TYPE_WIDTH)] = *it;
158 moab_tetgen_map[*it] = MBTRI | (ii << MB_TYPE_WIDTH);
159 for (int nn = 0; nn < 3; nn++) {
160 if (moab_tetgen_map.find(conn[order[nn]]) == moab_tetgen_map.end()) {
161 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
162 "data inconsistency between TetGen and MoAB");
163 }
164 in.trifacelist[3 * ii + nn] =
165 moab_tetgen_map[conn[order[nn]]] >> MB_TYPE_WIDTH;
166 }
167 }
168 }
169
170 Range edges = ents.subset_by_type(MBEDGE);
171 in.numberofedges = edges.size();
172 if (in.numberofedges > 0) {
173 in.edgelist = new int[2 * in.numberofedges];
174 in.edgemarkerlist = new int[in.numberofedges];
175 // std::fill(&in.edgemarkerlist[0],&in.edgemarkerlist[in.numberofedges],1);
176 CHKERR m_field.get_moab().tag_get_data(th_marker, edges, in.edgemarkerlist);
177 it = edges.begin();
178 for (int ii = 0; it != edges.end(); it++, ii++) {
179 int num_nodes;
180 const EntityHandle *conn;
181 CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
182 tetgen_moab_map[MBEDGE | (ii << MB_TYPE_WIDTH)] = *it;
183 moab_tetgen_map[*it] = MBEDGE | (ii << MB_TYPE_WIDTH);
184 for (int nn = 0; nn < 2; nn++) {
185 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
186 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
187 "data inconsistency between TetGen and MoAB");
188 }
189 in.edgelist[2 * ii + nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
190 }
191 }
192 }
193
195}
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352

◆ loadPoly()

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

load poly file

Definition at line 654 of file TetGenInterface.cpp.

654 {
656 if (!in.load_poly(file_name)) {
657 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
658 "can not read TetGen poly file");
659 }
661}
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34

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

845 {
847 // FIXME: assumes that are no holes
848
849 if (ents.empty()) {
850 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
851 "no ents to build polygon");
852 }
853
854 Interface &m_field = cOre;
855
856 Skinner skin(&m_field.get_moab());
857
858 Range skin_edges;
859 CHKERR skin.find_skin(0, ents, false, skin_edges);
860
861 std::vector<EntityHandle> polygon_nodes;
862 EntityHandle seed = skin_edges[0];
863 Range seen_edges;
864 seen_edges.insert(seed);
865 skin_edges.erase(seed);
866 int num_nodes;
867 const EntityHandle *conn;
868 CHKERR m_field.get_moab().get_connectivity(seed, conn, num_nodes, true);
869 polygon_nodes.push_back(conn[0]);
870 // std::cerr << std::endl;
871 // std::cerr << conn[0] << " " << conn[1] << std::endl;
872 do {
873 EntityHandle last_node = polygon_nodes.back();
874 Range adj_edges;
875 CHKERR m_field.get_moab().get_adjacencies(&last_node, 1, 1, false,
876 adj_edges);
877 adj_edges = intersect(adj_edges, skin_edges);
878 if (adj_edges.size() == 0) {
879 break;
880 }
881 if (adj_edges.size() != 1) {
882 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
883 "should be only one edge");
884 }
885 seen_edges.insert(adj_edges[0]);
886 skin_edges.erase(adj_edges[0]);
887 CHKERR m_field.get_moab().get_connectivity(adj_edges[0], conn, num_nodes,
888 true);
889 EntityHandle add_node = (last_node == conn[0]) ? conn[1] : conn[0];
890 polygon_nodes.push_back(add_node);
891 // std::cerr << "\t" << add_node << std::endl;
892 } while (1);
893
894 if (reduce_edges) {
895 // std::cerr << "polygon " << polygon_nodes.size();
896 std::vector<EntityHandle>::iterator pit = polygon_nodes.begin();
897 // std::cerr << std::endl;
898 for (; pit != polygon_nodes.end();) {
899 if (not_reducable_nodes != NULL) {
900 if (not_reducable_nodes->find(*pit) != not_reducable_nodes->end()) {
901 pit++;
902 continue;
903 }
904 }
905 EntityHandle mm;
906 if (pit == polygon_nodes.begin()) {
907 mm = polygon_nodes.back();
908 } else {
909 mm = *(pit - 1);
910 }
911 EntityHandle mc = *pit;
912 EntityHandle mp;
913 if (polygon_nodes.back() == mc) {
914 mp = polygon_nodes[0];
915 } else {
916 mp = *(pit + 1);
917 }
918 double coords[9];
919 if (th) {
920 CHKERR m_field.get_moab().tag_get_data(th, &mm, 1, &coords[3 * 0]);
921 CHKERR m_field.get_moab().tag_get_data(th, &mc, 1, &coords[3 * 1]);
922 CHKERR m_field.get_moab().tag_get_data(th, &mp, 1, &coords[3 * 2]);
923 } else {
924 CHKERR m_field.get_moab().get_coords(&mm, 1, &coords[3 * 0]);
925 CHKERR m_field.get_moab().get_coords(&mc, 1, &coords[3 * 1]);
926 CHKERR m_field.get_moab().get_coords(&mp, 1, &coords[3 * 2]);
927 }
928 cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 0], 1); // mm = mm -
929 // mc
930 cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 2], 1); // mp = mp -
931 // mc
932 double spin[9];
933 CHKERR Spin(spin, &coords[3 * 0]);
934 double l0 = cblas_dnrm2(3, &coords[3 * 0], 1);
935 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1. / l0, spin, 3,
936 &coords[3 * 2], 1, 0., &coords[3 * 1], 1);
937 double dot = cblas_dnrm2(3, &coords[3 * 1], 1);
938 // std::cerr << mm << " " << mc << " " << mp << " " << dot << std::endl;
939 if (dot < eps) {
940 polygon_nodes.erase(pit);
941 pit = polygon_nodes.begin();
942 // std::cerr << std::endl;
943 } else {
944 pit++;
945 }
946 }
947 }
948
949 Range existing_polygon;
950 CHKERR m_field.get_moab().get_adjacencies(
951 &polygon_nodes[0], polygon_nodes.size(), 2, true, existing_polygon);
952 if (existing_polygon.empty()) {
953 EntityHandle polygon;
954 CHKERR m_field.get_moab().create_element(MBPOLYGON, &polygon_nodes[0],
955 polygon_nodes.size(), polygon);
956 polygons.insert(polygon);
957 } else {
958 polygons.merge(existing_polygon);
959 }
960
962}
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:546

◆ outData() [1/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 439 of file TetGenInterface.cpp.

444 {
446 Interface &m_field = cOre;
447 Range ents;
448 CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
449 error_if_created, assume_first_nodes_the_same,th);
450 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
451 ents.subset_by_type(MBTET), bit);
453}
auto bit
set bit
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

◆ outData() [2/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 246 of file TetGenInterface.cpp.

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

◆ query_interface()

MoFEMErrorCode MoFEM::TetGenInterface::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 31 of file TetGenInterface.cpp.

32 {
33 *iface = const_cast<TetGenInterface *>(this);
34 return 0;
35}
TetGenInterface(const MoFEM::Core &core)

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

457 {
459 ErrorCode rval;
460 Interface &m_field = cOre;
461 in.numberoffacets = markers.size();
462 in.facetlist = new tetgenio::facet[in.numberoffacets];
463 in.facetmarkerlist = new int[in.numberoffacets];
464 std::vector<std::pair<Range, int>>::iterator mit = markers.begin();
465 for (int ii = 0; mit != markers.end(); mit++, ii++) {
466 in.facetmarkerlist[ii] = mit->second;
467 Range &faces = mit->first;
468 tetgenio::facet *f = &(in.facetlist[ii]);
469 f->numberofpolygons = faces.size();
470 f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
471 int jj = 0;
472 for (int dd = 3; dd >= 0; dd--) {
473 Range dd_faces = faces.subset_by_dimension(dd);
474 if (dd_faces.empty())
475 continue;
476 Range::iterator it = dd_faces.begin();
477 for (; it != dd_faces.end(); it++, jj++) {
478 int num_nodes;
479 const EntityHandle *conn;
480 tetgenio::polygon *p = &(f->polygonlist[jj]);
481 switch (type_from_handle(*it)) {
482 case MBVERTEX: {
483 p->numberofvertices = 1;
484 conn = &*it;
485 } break;
486 default: {
487 rval =
488 m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
490 p->numberofvertices = num_nodes;
491 }
492 }
493 p->vertexlist = new int[p->numberofvertices];
494 for (int nn = 0; nn < p->numberofvertices; nn++) {
495 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
496 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
497 "data inconsistency between TetGen and MoAB");
498 }
499 p->vertexlist[nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
500 }
501 }
502 }
503 // holes
504 f->numberofholes = 0;
505 f->holelist = NULL;
506 }
508}
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
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1634

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

200 {
202
203 Interface &m_field = cOre;
204 //
205 // ErrorCode rval;
206 in.pointparamlist = new tetgenio::pointparam[in.numberofpoints];
207 // std::vector<bool> points_is_set(in.numberofpoints,false);
208 std::map<int, Range>::iterator mit = type_ents.begin();
209 for (; mit != type_ents.end(); mit++) {
210 if (mit->first < 0 && mit->first > 3) {
211 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
212 "Wrong TetGen point type");
213 }
214 Range::iterator it = mit->second.begin();
215 for (; it != mit->second.end(); it++) {
216 moabTetGen_Map::iterator miit = moab_tetgen_map.find(*it);
217 if (miit == moab_tetgen_map.end()) {
218 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
219 "Data inconsistency between TetGen and MoAB");
220 continue;
221 }
222 int id = miit->second >> MB_TYPE_WIDTH;
223 in.pointparamlist[id].uv[0] = 0;
224 in.pointparamlist[id].uv[1] = 0;
225 in.pointparamlist[id].type = mit->first;
226 in.pointparamlist[id].tag = m_field.get_moab().id_from_handle(*it) + 1;
227 // points_is_set[id] = true;
228 }
229 }
230
231 // // Check only if tag and type is set to all points
232 // for(
233 // std::vector<bool>::iterator bit = points_is_set.begin();
234 // bit!=points_is_set.end();bit++
235 // ) {
236 // if(!*bit) {
237 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Point type for
238 // TetGen is not set");
239 // }
240 // }
241
243}

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

557 {
559 ErrorCode rval;
560 Interface &m_field = cOre;
561 in.numberofregions = regions.size();
562 in.regionlist = new double[5 * in.numberofregions];
563 int kk = 0;
564 std::vector<std::pair<EntityHandle, int>>::iterator it = regions.begin();
565 for (int ii = 0; it != regions.end(); it++, ii++) {
566 double coords[3];
567 switch (type_from_handle(it->first)) {
568 case MBVERTEX: {
569 if (th) {
570 rval = m_field.get_moab().tag_get_data(th, &it->first, 1, coords);
572 } else {
573 rval = m_field.get_moab().get_coords(&it->first, 1, coords);
575 }
576 } break;
577 case MBTET: {
578 int num_nodes;
579 const EntityHandle *conn;
580 rval =
581 m_field.get_moab().get_connectivity(it->first, conn, num_nodes, true);
583 double _coords[12];
584 if (th) {
585 rval = m_field.get_moab().tag_get_data(th, conn, num_nodes, _coords);
587 } else {
588 rval = m_field.get_moab().get_coords(conn, num_nodes, _coords);
590 }
591 coords[0] = (_coords[0] + _coords[3] + _coords[6] + _coords[9]) / 4.;
592 coords[1] = (_coords[1] + _coords[4] + _coords[7] + _coords[10]) / 4.;
593 coords[2] = (_coords[2] + _coords[5] + _coords[8] + _coords[11]) / 4.;
594 } break;
595 default:
596 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
597 }
598 for (int nn = 0; nn < 3; nn++) {
599 in.regionlist[kk++] = coords[nn];
600 }
601 in.regionlist[kk++] = it->second;
602 in.regionlist[kk++] = it->second;
603 }
605}
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32

◆ tetRahedralize()

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

run tetgen

Definition at line 646 of file TetGenInterface.cpp.

647 {
649 tetgenbehavior a;
650 a.parse_commandline(switches);
651 tetrahedralize(&a, &in, &out);
653}
constexpr double a

Member Data Documentation

◆ cOre

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

Definition at line 28 of file TetGenInterface.hpp.


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