Branch data Line data Source code
1 : :
2 : : #include "meshkit/MKCore.hpp"
3 : : #include "meshkit/CurveFacetMeshReader.hpp"
4 : : #include "meshkit/ModelEnt.hpp"
5 : : #include <iostream>
6 : : #include <iGeom.h>
7 : : #include <cmath>
8 : : #include "MBTagConventions.hpp"
9 : : #include "moab/Core.hpp"
10 : : #include "moab/GeomTopoTool.hpp"
11 : :
12 : : namespace MeshKit
13 : : {
14 : : // Construction Function for CurveFacetMeshReader
15 : :
16 : : moab::EntityType CurveFacetMeshReader_types[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBMAXTYPE};
17 : 40 : const moab::EntityType* CurveFacetMeshReader::output_types()
18 : 40 : { return CurveFacetMeshReader_types; }
19 : :
20 : :
21 : 2 : CurveFacetMeshReader::CurveFacetMeshReader(MKCore *mk_core, const MEntVector &ments)
22 : 2 : : MeshScheme(mk_core, ments)
23 : : {
24 : 2 : mk = mk_core;
25 : 2 : }
26 : :
27 : : // Destructor Function for CurveFacetMeshReader
28 : 0 : CurveFacetMeshReader::~CurveFacetMeshReader()
29 : : {
30 [ # # ]: 0 : }
31 : :
32 : 2 : void CurveFacetMeshReader::setup_this()
33 : : {
34 : :
35 : : //create a vertex mesher
36 : 2 : int meshop_dim = 0;
37 [ + - ][ + - ]: 2 : MeshOp *vm = (MeshOp*) mk_core()->construct_meshop(meshop_dim);
38 : :
39 [ + - ][ + - ]: 32 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
40 : : {
41 [ + - ]: 30 : ModelEnt *me = mit->first;
42 : : //do an initial check that the model entity is of the correct dimension
43 [ + - ][ - + ]: 30 : if(me->dimension() != 1)
44 : : {
45 [ # # ][ # # ]: 0 : std::cout << "Found an entity that is of an incorrect dimension" << std::endl;
46 [ # # ]: 0 : mentSelection.erase(mit);
47 : 0 : continue;
48 : : }
49 : :
50 : :
51 : : //make sure that its children (vertices) have been meshed
52 [ + - ]: 30 : MEntVector children;
53 [ + - ]: 30 : me->get_adjacencies(0, children);
54 : :
55 [ + - ][ + - ]: 90 : for(MEntVector::iterator j = children.begin(); j != children.end(); j++)
[ + + ]
56 : : {
57 [ + - ]: 60 : ModelEnt *child = *j;
58 [ + - ][ + + ]: 60 : if(child->is_meshops_list_empty()) vm->add_modelent(child);
[ + - ]
59 : : }
60 : :
61 : 30 : }
62 : :
63 : : //make sure we're meshing the veritces first
64 : 2 : mk_core()->insert_node(vm, this, mk_core()->root_node());
65 : :
66 : : //set faceting parameters to default values if the are not already set.
67 : :
68 : : //set the faceting_tolerance
69 [ - + ]: 2 : if(!facet_tol) facet_tol = 1e-4;
70 : :
71 : : //set the geom_reabs value
72 [ - + ]: 2 : if(!geom_res) geom_res = 1e-6;
73 : :
74 : 2 : }
75 : :
76 : 2 : void CurveFacetMeshReader::execute_this()
77 : : {
78 : :
79 [ + - ][ + - ]: 32 : for(MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
80 : : {
81 : :
82 [ + - ]: 30 : ModelEnt *me = mit->first;
83 : :
84 : :
85 : : //---------------------------------//
86 : : // FACET //
87 : : //---------------------------------//
88 [ + - ]: 30 : facet(me);
89 : :
90 : : // Not sure if this needs to be done yet
91 : : //set_senses(me);
92 : :
93 [ + - ]: 30 : me->set_meshed_state(COMPLETE_MESH);
94 : : }
95 : 2 : }
96 : :
97 : 30 : void CurveFacetMeshReader::facet(ModelEnt *curve)
98 : : {
99 [ + - ]: 30 : iGeom::EntityHandle h = curve->geom_handle();
100 [ + - ]: 30 : iMesh::EntitySetHandle sh = IBSH(curve->mesh_handle());
101 : :
102 : : //get the facets for this curve/ref_edge
103 [ + - ]: 30 : std::vector<double> pnts;
104 [ + - ][ + - ]: 60 : std::vector<int> conn;
105 [ + - ][ + - ]: 30 : iGeom::Error errval = mk_core()->igeom_instance()->getFacets(h,facet_tol,pnts,conn);
[ + - ]
106 [ - + ][ # # ]: 30 : if (errval != iBase_SUCCESS) throw Error(MK_FAILURE, "Unable get the facets of the curve.");
107 : : //create vector for keeping track of the vertices
108 [ + - ][ + - ]: 60 : std::vector<iBase_EntityHandle> verts;
109 : :
110 : :
111 : : // loop over the facet points
112 [ + + ]: 1078 : for(unsigned int j = 0; j < pnts.size(); j+=3)
113 : : {
114 : : //create a new vertex handle
115 : : iBase_EntityHandle v;
116 [ + - ][ + - ]: 1048 : mk_core()->imesh_instance()->createVtx(pnts[j],pnts[j+1],pnts[j+2],v);
[ + - ][ + - ]
[ + - ][ + - ]
117 [ + - ]: 1048 : verts.push_back(v);
118 : : }
119 : :
120 : : //--------------------------------------------------//
121 : : //check that the start and end vertices are within //
122 : : //the geom resabs of the facet verts //
123 : : //--------------------------------------------------//
124 : :
125 : : //get the curve vertices
126 [ + - ][ + - ]: 60 : MEntVector end_verts;
127 [ + - ]: 30 : curve->get_adjacencies(0, end_verts);
128 : :
129 : : //get the entity handle of the start vertex
130 [ + - ][ + - ]: 60 : std::vector<moab::EntityHandle> dum_handles;
131 [ + - ][ + - ]: 30 : end_verts.front()->get_mesh(0,dum_handles);
132 [ - + ]: 30 : assert( 1 == dum_handles.size() ); // should only be one vertex in a vertex entity set
133 [ + - ]: 30 : iMesh::EntityHandle start_vtx = IBEH(dum_handles.front());
134 : :
135 : 30 : dum_handles.clear();
136 : :
137 : : // do the same for the end vertex
138 [ + - ][ + - ]: 30 : end_verts.back()->get_mesh(0, dum_handles);
139 [ - + ]: 30 : assert( 1 == dum_handles.size() ); // should only be one vertex in a vertex entity set
140 [ + - ]: 30 : iMesh::EntityHandle end_vtx = IBEH(dum_handles.front());
141 : :
142 : 30 : dum_handles.clear();
143 : :
144 : : //std::cout << "Number of end_verts: " << end_verts.size() << std::endl;
145 : :
146 : : //special case for a point curve
147 [ - + ]: 30 : if(verts.size() < 2)
148 : : {
149 [ # # ][ # # ]: 0 : if( start_vtx != end_vtx || curve->measure() > geom_res)
[ # # ][ # # ]
150 : : {
151 [ # # ][ # # ]: 0 : std::cout << "Warning: no facetting for curve " << curve->id() << std::endl;
[ # # ][ # # ]
152 : 0 : return;
153 : : }
154 : : else
155 : : {
156 [ # # ][ # # ]: 0 : mk_core()->imesh_instance()->addEntToSet(start_vtx, sh);
[ # # ]
157 : 0 : return;
158 : : }
159 : : }
160 : :
161 : : //check for a closed curve
162 [ + - ][ + - ]: 30 : bool closed = (mvtx2mvtx_dist(verts.front(),verts.back()) < geom_res);
[ + - ]
163 : :
164 : : // check that geometry and facetting agree that the curve is closed
165 [ - + ]: 30 : if ( closed != (start_vtx == end_vtx))
166 [ # # ]: 0 : std::cout << "Warning: topology and geometry inconsistent for possibly closed curve "
167 [ # # ][ # # ]: 0 : << curve->id() << std::endl;
[ # # ]
168 : :
169 : : //check the proximity of the front and end vertices to the start and end points, respectively
170 [ + - ][ + - ]: 60 : if(vtx2vtx_dist(end_verts.front()->geom_handle(), verts.front()) > geom_res ||
[ + - ][ + - ]
[ + - ][ - + ]
[ - + ]
171 [ + - ][ + - ]: 30 : vtx2vtx_dist(end_verts.back()->geom_handle(), verts.back()) > geom_res)
[ + - ][ + - ]
172 : : {
173 : : //try reversing the points
174 : : //std::reverse(verts.begin(), verts.end());
175 : :
176 : : //check again, if this time it fails, give a warning
177 [ # # ][ # # ]: 0 : if(vtx2vtx_dist(end_verts.front()->geom_handle(), verts.front()) > geom_res ||
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
178 [ # # ][ # # ]: 0 : vtx2vtx_dist(end_verts.back()->geom_handle(), verts.back()) > geom_res)
[ # # ][ # # ]
179 : : {
180 [ # # ][ # # ]: 0 : std::cerr << "Warning: vertices not at the ends of the curve " << curve->id() << std::endl;
[ # # ][ # # ]
181 : : }
182 : : }
183 : :
184 : :
185 : : //now replace start and end facet points with the curve's child vert(s)
186 : :
187 : : //capture the beginning and end vertex handles for deletion
188 [ + - ]: 30 : iMesh::EntityHandle front_vert_to_del = verts.front();
189 [ + - ]: 30 : iMesh::EntityHandle back_vert_to_del = verts.back();
190 : :
191 [ + - ]: 30 : verts.front() = start_vtx;
192 [ + - ]: 30 : verts.back() = end_vtx;
193 : :
194 : : //delete the old vertices
195 [ + - ][ + - ]: 30 : mk_core()->imesh_instance()->deleteEnt(front_vert_to_del);
[ + - ]
196 [ + - ][ + - ]: 30 : mk_core()->imesh_instance()->deleteEnt(back_vert_to_del);
[ + - ]
197 : :
198 : : //now create the edges with the vertices we've chosen
199 : :
200 : : //vector to contain the edges
201 [ + - ][ + - ]: 60 : std::vector<iBase_EntityHandle> edges;
202 : : //loop over vertices and create edges
203 [ + + ]: 1048 : for (unsigned int j = 0; j < verts.size()-1; j++)
204 : : {
205 : : //create edges
206 : : iBase_EntityHandle e;
207 [ + - ][ + - ]: 1018 : mk_core()->imesh_instance()->createEnt(iMesh_LINE_SEGMENT, &verts[j], 2, e);
[ + - ][ + - ]
208 [ + - ]: 1018 : edges.push_back(e);
209 : : }
210 : :
211 : : //remove duplicate vertex if the curve is closed
212 [ - + ][ # # ]: 30 : if(closed && start_vtx == end_vtx) verts.pop_back();
[ # # ]
213 : :
214 : : //add vertices and edges to the entity set
215 [ + - ][ + - ]: 30 : mk_core()->imesh_instance()->addEntArrToSet(&verts[0], verts.size(), sh);
[ + - ][ + - ]
216 [ + - ][ + - ]: 60 : mk_core()->imesh_instance()->addEntArrToSet(&edges[0], edges.size(), sh);
[ + - ][ + - ]
217 : : }
218 : :
219 : 0 : void CurveFacetMeshReader::set_senses(ModelEnt *curve)
220 : : {
221 : :
222 : : //get the geom_handle for this curve
223 [ # # ]: 0 : iGeom::EntityHandle gh = curve->geom_handle();
224 : :
225 : : //get all surfaces adjacent to this curve
226 [ # # ]: 0 : MEntVector adj_surfs;
227 [ # # ]: 0 : curve->get_adjacencies(2, adj_surfs);
228 : :
229 : 0 : MEntVector::iterator i;
230 [ # # ]: 0 : std::vector<int> senses;
231 [ # # ]: 0 : std::vector<moab::EntityHandle> meshsets;
232 [ # # ][ # # ]: 0 : for(i = adj_surfs.begin(); i !=adj_surfs.end(); i++)
[ # # ]
233 : : {
234 : : int sense;
235 : : //create a model ent for the current adjacent surface
236 [ # # ]: 0 : ModelEnt *adj_ent = *i;
237 : : //get the sense of the curve wrt the surface
238 [ # # ][ # # ]: 0 : mk_core()->igeom_instance()->getEgFcSense(gh, adj_ent->geom_handle(),sense);
[ # # ][ # # ]
239 : : //std::cout << "Sense: " << sense << std::endl;
240 : : // add the sense and the entityhandle to the appropriate vectors
241 [ # # ]: 0 : senses.push_back(sense);
242 [ # # ][ # # ]: 0 : meshsets.push_back(adj_ent->mesh_handle());
243 : : }
244 : :
245 : : // use moab to set the senses of the entities
246 : : // using moab to do this allows the use of variable length tags
247 : : // which are not currently available in iGeom
248 [ # # ][ # # ]: 0 : moab::GeomTopoTool gt(mk_core()->moab_instance());
[ # # ]
249 [ # # ][ # # ]: 0 : moab::ErrorCode rval = gt.set_senses(curve->mesh_handle(),meshsets,senses);
250 [ # # ][ # # ]: 0 : if(rval != moab::MB_SUCCESS) std::cout << "Error setting curve senses!" << std::endl;
[ # # ]
251 : :
252 : 0 : }
253 : :
254 : : // should I be using an iMesh entity to do this comparison??
255 : 60 : double CurveFacetMeshReader::vtx2vtx_dist(iGeom::EntityHandle vtx1, iMesh::EntityHandle vtx2)
256 : : {
257 : :
258 : : double x1,y1,z1;
259 : : double x2,y2,z2;
260 : :
261 [ + - ][ + - ]: 60 : mk_core()->igeom_instance()->getVtxCoord(vtx1, x1, y1, z1);
[ + - ]
262 [ + - ][ + - ]: 60 : mk_core()->imesh_instance()->getVtxCoord(vtx2, x2, y2, z2);
[ + - ]
263 : :
264 [ + - ][ + - ]: 60 : double dist = pow((x1-x2),2) + pow((y1-y2),2) + pow((z1-z2),2);
[ + - ]
265 : :
266 : 60 : return sqrt(dist);
267 : : }
268 : :
269 : 30 : double CurveFacetMeshReader::mvtx2mvtx_dist(iMesh::EntityHandle vtx1, iMesh::EntityHandle vtx2)
270 : : {
271 : :
272 : : double x1,y1,z1;
273 : : double x2,y2,z2;
274 : :
275 [ + - ][ + - ]: 30 : mk_core()->imesh_instance()->getVtxCoord(vtx1, x1, y1, z1);
[ + - ]
276 [ + - ][ + - ]: 30 : mk_core()->imesh_instance()->getVtxCoord(vtx2, x2, y2, z2);
[ + - ]
277 : :
278 [ + - ][ + - ]: 30 : double dist = pow((x1-x2),2) + pow((y1-y2),2) + pow((z1-z2),2);
[ + - ]
279 : :
280 : 30 : return sqrt(dist);
281 : : }
282 : :
283 : 2 : void CurveFacetMeshReader::set_mesh_params(double faceting_tolerance, double geom_resabs)
284 : : {
285 : : //Assign the faceting values if they are passed into the function
286 [ + - ]: 2 : if(faceting_tolerance) facet_tol = faceting_tolerance;
287 [ + - ]: 2 : if(geom_resabs) geom_res = geom_resabs;
288 : 2 : }
289 : :
290 [ + - ][ + - ]: 156 : }
|
