Branch data Line data Source code
1 : : #include "meshkit/NGTetMesher.hpp"
2 : : #include "meshkit/MeshScheme.hpp"
3 : : #include "meshkit/ModelEnt.hpp"
4 : : #include "moab/Interface.hpp"
5 : : #include "moab/ReadUtilIface.hpp"
6 : : #include "moab/Range.hpp"
7 : : #include "meshkit/RegisterMeshOp.hpp"
8 : : namespace nglib
9 : : {
10 : : #include "nglib.h"
11 : : }
12 : :
13 : : #include <vector>
14 : :
15 : : namespace MeshKit
16 : : {
17 : :
18 : : moab::EntityType NGTetMesher::meshTps[] = {moab::MBVERTEX, moab::MBTET, moab::MBMAXTYPE};
19 : :
20 : 2 : NGTetMesher::NGTetMesher(MKCore *mk_core, const MEntVector &me_vec)
21 : 2 : : MeshScheme(mk_core, me_vec)
22 : : {
23 : 2 : }
24 : :
25 : 6 : NGTetMesher::~NGTetMesher()
26 : : {
27 [ - + ]: 4 : }
28 : :
29 : 2 : MeshOp *NGTetMesher::get_tri_mesher()
30 : : {
31 [ + - ]: 2 : std::vector<MeshOpProxy *> proxies;
32 [ + - ][ + - ]: 2 : mk_core()->meshop_by_mesh_type(moab::MBTRI, proxies);
33 [ - + ][ # # ]: 2 : if (proxies.empty()) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing triangles.");
34 [ + - ][ + - ]: 2 : return mk_core()->construct_meshop("AF2DfltTriangleMeshOp");
[ + - ][ + - ]
35 : : }
36 : :
37 : 2 : void NGTetMesher::setup_this()
38 : : {
39 : 2 : MeshOp *tri_mesher = NULL;
40 [ + - ]: 2 : MEntVector surfs;
41 : :
42 [ + - ][ + - ]: 4 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + + ]
43 [ + - ]: 2 : ModelEnt *me = (*sit).first;
44 [ + - ][ - + ]: 2 : if (me->dimension() != 3) throw Error(MK_BAD_INPUT, "Tet mesher assigned to an entity with dimension != 3.");
[ # # ]
45 : :
46 : : // get the bounding faces
47 : 2 : surfs.clear();
48 [ + - ]: 2 : me->get_adjacencies(2, surfs);
49 : :
50 : : // check the mesh scheme on each one; if there is one, verify it can generate tris; if there
51 : : // isn't one, make one
52 : 2 : bool inserted = false;
53 : :
54 [ + - ][ + - ]: 24 : for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) {
[ + + ]
55 [ + - ][ + - ]: 22 : if ((*vit)->is_meshops_list_empty()) {
[ + - ]
56 : : // get a tri mesher if we haven't already
57 [ + + ][ + - ]: 22 : if (!tri_mesher) tri_mesher = get_tri_mesher();
58 : :
59 : : // add this surface to it, and if first for the volume, make sure it's added upstream
60 [ + - ][ + - ]: 22 : tri_mesher->add_modelent(*vit);
61 [ + + ]: 22 : if (!inserted) {
62 [ + - ][ + - ]: 2 : mk_core()->insert_node(tri_mesher, this);
63 : 2 : inserted = true;
64 : : }
65 : : } // if no meshops
66 : : } // over surfs
67 : 2 : } // over vols
68 : 2 : }
69 : :
70 : 2 : void NGTetMesher::execute_this()
71 : : {
72 : 2 : nglib::Ng_Init();
73 : :
74 [ + - ][ + - ]: 4 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + + ]
75 : : // make a me, for convenience
76 [ + - ]: 2 : ModelEnt *me = (*sit).first;
77 : :
78 : : // assemble bounding mesh
79 [ + - ]: 2 : std::vector<moab::EntityHandle> bdy;
80 [ + - ][ + - ]: 4 : std::vector<int> group_sizes, senses, bdy_ids;
[ + - ]
81 [ + - ]: 2 : me->boundary(2, bdy, &senses, &group_sizes);
82 : :
83 : : // get connectivity
84 : :
85 : : // convert the handles to integers for input to CMLTetMesher
86 [ + - ]: 4 : moab::Range bdy_vrange;
87 [ + - ]: 4 : std::vector<double> coords;
88 [ + - ]: 2 : me->get_indexed_connect_coords(bdy, &senses, NULL, bdy_ids, coords, &bdy_vrange, 1);
89 : :
90 [ + - ]: 2 : nglib::Ng_Mesh *ngmesh = nglib::Ng_NewMesh();
91 [ + - ]: 2 : unsigned int numvs = bdy_vrange.size();
92 : :
93 : : // add the points
94 [ + + ]: 2264 : for (unsigned int i = 0; i < numvs; i++)
95 [ + - ][ + - ]: 2262 : nglib::Ng_AddPoint(ngmesh, &coords[3*i]);
96 : :
97 : : // add the tris
98 : 2 : unsigned int numtris = bdy.size();
99 [ + + ]: 4538 : for (unsigned int i = 0; i < numtris; i++)
100 [ + - ][ + - ]: 4536 : nglib::Ng_AddSurfaceElement(ngmesh, nglib::NG_TRIG, &bdy_ids[3*i]);
101 : :
102 [ + - ]: 2 : double my_size = me->mesh_interval_size();
103 [ - + ]: 2 : if (0 > my_size) my_size = 1.0;
104 [ + - ]: 2 : nglib::Ng_RestrictMeshSizeGlobal(ngmesh, my_size);
105 : :
106 [ + - ]: 2 : nglib::Ng_Meshing_Parameters ngp;
107 : 2 : ngp.maxh = my_size;
108 : 2 : ngp.fineness = 0.5;
109 : :
110 : : // Use variable name second_order if using newer version of NetGen
111 : : //ngp.secondorder = 0;
112 : 2 : ngp.second_order = 0;
113 : : // nglib.h Rev. 663 http://netgen-mesher.svn.sourceforge.net/viewvc/netgen-mesher/netgen/nglib/nglib.h?revision=663&view=markup
114 : : //uses: second_order instead of secondorder(4.9.13)
115 : :
116 : :
117 [ + - ]: 2 : nglib::Ng_Result result = nglib::Ng_GenerateVolumeMesh(ngmesh, &ngp);
118 [ - + ][ # # ]: 2 : if (nglib::NG_OK != result) ECERRCHK(MK_FAILURE, "Netgen mesher returned !ok.");
[ # # ]
119 : :
120 : : moab::ErrorCode rval;
121 [ + - ]: 4 : moab::Range new_ents;
122 [ + - ]: 2 : int num_pts = nglib::Ng_GetNP(ngmesh);
123 [ - + ]: 2 : assert(num_pts >= (int)numvs);
124 [ + - ]: 2 : if (num_pts > (int)numvs) {
125 [ + - ]: 2 : coords.resize(3*(num_pts-numvs));
126 : :
127 [ + + ]: 1373 : for (unsigned int i = (unsigned int) numvs; i < (unsigned int) num_pts; i++)
128 [ + - ][ + - ]: 1371 : nglib::Ng_GetPoint(ngmesh, i, &coords[3*(i-numvs)]);
129 : :
130 : : // create the new vertices' entities
131 [ + - ][ + - ]: 2 : rval = mk_core()->moab_instance()->create_vertices(&coords[0], num_pts-numvs, new_ents);
[ + - ][ + - ]
132 [ - + ][ # # ]: 2 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
133 : : }
134 : :
135 : : // for tets, pre-allocate connectivity
136 : : moab::ReadUtilIface *iface;
137 [ + - ][ + - ]: 2 : rval = mk_core()-> moab_instance() -> query_interface(iface);
[ + - ]
138 [ - + ][ # # ]: 2 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
139 : :
140 : : //create the tris, get a direct ptr to connectivity
141 [ + - ]: 2 : int num_tets = (int) nglib::Ng_GetNE(ngmesh);
142 : : moab::EntityHandle starth, *connect;
143 [ + - ]: 2 : rval = iface->get_element_connect(num_tets, 4, moab::MBTET, 1, starth, connect);
144 [ - + ][ # # ]: 2 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
145 : :
146 : : // read connectivity directly into that array, as int's
147 : 2 : int *connecti = (int*) connect;
148 [ + + ]: 14465 : for (unsigned int i = 0; i < (unsigned int) num_tets; i++)
149 [ + - ]: 14463 : nglib::Ng_GetVolumeElement(ngmesh, i+1, connecti+4*i);
150 : :
151 : : // put vertex handles into an indexible array
152 [ + - ]: 4 : std::vector<moab::EntityHandle> bdy_vvec;
153 [ + - ][ + - ]: 2 : std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
154 [ + - ][ + - ]: 2 : std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
155 : :
156 : : // now convert vertex indices into handles in-place, working from the back
157 [ + + ]: 57854 : for (int i = 4*num_tets-1; i >= 0; i--) {
158 [ + - ][ - + ]: 57852 : assert(connecti[i] > 0 && connecti[i] <= (int)bdy_vvec.size());
159 [ + - ]: 57852 : connect[i] = bdy_vvec[connecti[i]-1];
160 : : }
161 : :
162 : : // put new tris into new entity range, then commit the mesh
163 [ + - ]: 2 : new_ents.insert(starth, starth+num_tets-1);
164 [ + - ]: 2 : me->commit_mesh(new_ents, COMPLETE_MESH);
165 : :
166 [ + - ]: 2 : nglib::Ng_DeleteMesh(ngmesh);
167 : 2 : }
168 : :
169 : 2 : nglib::Ng_Exit();
170 : 2 : }
171 : :
172 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
