Branch data Line data Source code
1 : : #include <iostream>
2 : : #include <map>
3 : :
4 : : #include "moab/FBEngine.hpp"
5 : : #include "moab/Interface.hpp"
6 : : #include "moab/GeomTopoTool.hpp"
7 : : #include "moab/GeomUtil.hpp"
8 : : #include "moab/OrientedBoxTreeTool.hpp"
9 : :
10 : : #include <stdlib.h>
11 : : #include <cstring>
12 : : #include <map>
13 : : #include <set>
14 : : #include <queue>
15 : : #include <algorithm>
16 : : #include "assert.h"
17 : :
18 : : #include "SmoothCurve.hpp"
19 : : #include "SmoothFace.hpp"
20 : :
21 : : // this is just to replace MBI with moab interface, which is _mbImpl in this class
22 : : #define MBI _mbImpl
23 : : #define MBERRORR( rval, STR ) \
24 : : { \
25 : : if( MB_SUCCESS != rval ) \
26 : : { \
27 : : std::cout << STR << std::endl; \
28 : : return rval; \
29 : : } \
30 : : }
31 : :
32 : : namespace moab
33 : : {
34 : :
35 : : // some tolerances for ray tracing and geometry intersections
36 : : // these are involved in ray tracing, at least
37 : :
38 : : double tolerance = 0.01; // TODO: how is this used ????
39 : : double tolerance_segment = 0.01; // for segments intersection, points collapse, area coordinates for triangles
40 : : // it should be a relative, not an absolute value
41 : : // as this splitting operation can create small edges, it should be relatively high
42 : : // or, it should be coordinated with the decimation errors
43 : : // we really just want to preserve the integrity of the mesh, we should avoid creating small edges
44 : : // or angles
45 : : const bool Debug_surf_eval = false;
46 : : bool debug_splits = false;
47 : :
48 : : // will compute intersection between a segment and slice of a plane
49 : : // output is the intersection point
50 : 2834 : bool intersect_segment_and_plane_slice( CartVect& from, CartVect& to, CartVect& p1, CartVect& p2, CartVect&,
51 : : CartVect& normPlane, CartVect& intx_point, double& parPos )
52 : : {
53 : : //
54 : : // plane eq is normPlane % r + d = 0, or normPlane % r - normPlane%p1 = 0
55 [ + - ]: 2834 : double dd = -normPlane % p1;
56 : 2834 : double valFrom = normPlane % from + dd;
57 : 2834 : double valTo = normPlane % to + dd;
58 : :
59 [ - + ]: 2834 : if( fabs( valFrom ) < tolerance_segment )
60 : : {
61 : 0 : intx_point = from;
62 : 0 : parPos = 0.;
63 [ # # ][ # # ]: 0 : double proj1 = ( intx_point - p1 ) % ( p2 - p1 );
64 [ # # ][ # # ]: 0 : double proj2 = ( intx_point - p2 ) % ( p1 - p2 );
65 [ # # ][ # # ]: 0 : if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false;
66 [ # # ]: 0 : if( debug_splits ) std::cout << "intx : " << intx_point << "\n";
67 : 0 : return true;
68 : : }
69 [ + + ]: 2834 : if( fabs( valTo ) < tolerance_segment )
70 : : {
71 : 15 : intx_point = to;
72 : 15 : parPos = 1;
73 [ + - ][ + - ]: 15 : double proj1 = ( intx_point - p1 ) % ( p2 - p1 );
74 [ + - ][ + - ]: 15 : double proj2 = ( intx_point - p2 ) % ( p1 - p2 );
75 [ + - ][ - + ]: 15 : if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false;
76 [ - + ]: 15 : if( debug_splits ) std::cout << "intx : " << intx_point << "\n";
77 : 15 : return true;
78 : : }
79 [ + + ]: 2819 : if( valFrom * valTo > 0 ) return false; // no intersection, although it could be very close
80 : : // else, it could intersect the plane; check for the slice too.
81 : 1755 : parPos = valFrom / ( valFrom - valTo ); // this is 0 for valFrom 0, 1 for valTo 0
82 [ + - ][ + - ]: 1755 : intx_point = from + ( to - from ) * parPos;
83 : : // now check if the intx_point is indeed between p1 and p2 in the slice.
84 [ + - ][ + - ]: 1755 : double proj1 = ( intx_point - p1 ) % ( p2 - p1 );
85 [ + - ][ + - ]: 1755 : double proj2 = ( intx_point - p2 ) % ( p1 - p2 );
86 [ + + ][ - + ]: 1755 : if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false;
87 : :
88 [ - + ]: 1747 : if( debug_splits ) std::cout << "intx : " << intx_point << "\n";
89 : 2834 : return true;
90 : : }
91 : :
92 : 15 : ErrorCode area_coordinates( Interface* mbi, EntityHandle tri, CartVect& pnt, double* area_coord,
93 : : EntityHandle& boundary_handle )
94 : : {
95 : :
96 : : int nnodes;
97 : : const EntityHandle* conn3;
98 [ + - ]: 15 : ErrorCode rval = mbi->get_connectivity( tri, conn3, nnodes );
99 [ - + ][ # # ]: 15 : MBERRORR( rval, "Failed to get connectivity" );
[ # # ]
100 [ - + ]: 15 : assert( 3 == nnodes );
101 [ + - ][ + + ]: 60 : CartVect P[3];
102 [ + - ]: 15 : rval = mbi->get_coords( conn3, nnodes, (double*)&P[0] );
103 [ - + ][ # # ]: 15 : MBERRORR( rval, "Failed to get coordinates" );
[ # # ]
104 : :
105 [ + - ]: 15 : CartVect r0( P[0] - pnt );
106 [ + - ]: 15 : CartVect r1( P[1] - pnt );
107 [ + - ]: 15 : CartVect r2( P[2] - pnt );
108 [ - + ]: 15 : if( debug_splits )
109 : : {
110 [ # # ][ # # ]: 0 : std::cout << " nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
111 [ # # ][ # # ]: 0 : std::cout << " distances: " << r0.length() << " " << r1.length() << " " << r2.length() << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
112 : : }
113 [ + - ][ + + ]: 15 : if( r0.length() < tolerance_segment )
114 : : {
115 : 1 : area_coord[0] = 1.;
116 : 1 : area_coord[1] = 0.;
117 : 1 : area_coord[2] = 0.;
118 : 1 : boundary_handle = conn3[0];
119 : 1 : return MB_SUCCESS;
120 : : }
121 [ + - ][ + + ]: 14 : if( r1.length() < tolerance_segment )
122 : : {
123 : 2 : area_coord[0] = 0.;
124 : 2 : area_coord[1] = 1.;
125 : 2 : area_coord[2] = 0.;
126 : 2 : boundary_handle = conn3[1];
127 : 2 : return MB_SUCCESS;
128 : : }
129 [ + - ][ + + ]: 12 : if( r2.length() < tolerance_segment )
130 : : {
131 : 3 : area_coord[0] = 0.;
132 : 3 : area_coord[1] = 0.;
133 : 3 : area_coord[2] = 1.;
134 : 3 : boundary_handle = conn3[2];
135 : 3 : return MB_SUCCESS;
136 : : }
137 : :
138 [ + - ]: 9 : CartVect v1( P[1] - P[0] );
139 [ + - ]: 9 : CartVect v2( P[2] - P[0] );
140 : :
141 [ + - ][ + - ]: 9 : double areaDouble = ( v1 * v2 ).length(); // the same for CartVect
142 [ - + ][ # # ]: 9 : if( areaDouble < tolerance_segment * tolerance_segment ) { MBERRORR( MB_FAILURE, "area of triangle too small" ); }
[ # # ]
143 [ + - ][ + - ]: 9 : area_coord[0] = ( r1 * r2 ).length() / areaDouble;
144 [ + - ][ + - ]: 9 : area_coord[1] = ( r2 * r0 ).length() / areaDouble;
145 [ + - ][ + - ]: 9 : area_coord[2] = ( r0 * r1 ).length() / areaDouble;
146 : :
147 [ - + ]: 9 : if( fabs( area_coord[0] + area_coord[1] + area_coord[2] - 1 ) > tolerance_segment )
148 [ # # ][ # # ]: 0 : { MBERRORR( MB_FAILURE, "point outside triangle" ); }
149 : : // the tolerance is used here for area coordinates (0 to 1), and in other
150 : : // parts it is used as an absolute distance; pretty inconsistent.
151 : 9 : bool side0 = ( area_coord[0] < tolerance_segment );
152 : 9 : bool side1 = ( area_coord[1] < tolerance_segment );
153 : 9 : bool side2 = ( area_coord[2] < tolerance_segment );
154 [ + - ][ + - ]: 9 : if( !side0 && !side1 && !side2 ) return MB_SUCCESS; // interior point
[ + + ]
155 : : // now, find out what boundary is in question
156 : : // first, get all edges, in order
157 [ + - ]: 2 : std::vector< EntityHandle > edges;
158 : : EntityHandle nn2[2];
159 [ + + ]: 8 : for( int i = 0; i < 3; i++ )
160 : : {
161 : 6 : nn2[0] = conn3[( i + 1 ) % 3];
162 : 6 : nn2[1] = conn3[( i + 2 ) % 3];
163 [ + - ]: 6 : std::vector< EntityHandle > adjacent;
164 [ + - ]: 6 : rval = mbi->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT );
165 [ - + ][ # # ]: 6 : MBERRORR( rval, "Failed to get edges" );
[ # # ]
166 [ - + ][ # # ]: 6 : if( adjacent.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges" );
[ # # ]
167 : : // should be only one edge here
168 [ + - ][ + - ]: 6 : edges.push_back( adjacent[0] );
[ + - ]
169 : 6 : }
170 : :
171 [ - + ][ # # ]: 2 : if( side0 ) boundary_handle = edges[0];
172 [ - + ][ # # ]: 2 : if( side1 ) boundary_handle = edges[1];
173 [ + - ][ + - ]: 2 : if( side2 ) boundary_handle = edges[2];
174 : :
175 : 15 : return MB_SUCCESS;
176 : : }
177 : :
178 : 12 : FBEngine::FBEngine( Interface* impl, GeomTopoTool* topoTool, const bool smooth )
179 : : : _mbImpl( impl ), _my_geomTopoTool( topoTool ), _t_created( false ), _smooth( smooth ), _initialized( false ),
180 [ + - ][ + + ]: 72 : _smthFace( NULL ), _smthCurve( NULL )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + -
# # # # #
# ]
181 : : {
182 [ + + ]: 12 : if( !_my_geomTopoTool )
183 : : {
184 [ + - ][ + - ]: 11 : _my_geomTopoTool = new GeomTopoTool( _mbImpl );
185 : 11 : _t_created = true;
186 : : }
187 : : // should this be part of the constructor or not?
188 : : // Init();
189 [ # # ]: 12 : }
190 : 84 : FBEngine::~FBEngine()
191 : : {
192 : 12 : clean();
193 : 12 : _smooth = false;
194 [ + - ][ + + ]: 72 : }
195 : :
196 : 12 : void FBEngine::clean()
197 : : {
198 [ + + ]: 12 : if( _smooth )
199 : : {
200 : 8 : _faces.clear();
201 : 8 : _edges.clear();
202 : 8 : int size1 = _my_gsets[1].size();
203 : 8 : int i = 0;
204 [ + + ]: 38 : for( i = 0; i < size1; i++ )
205 [ + - ]: 30 : delete _smthCurve[i];
206 [ + - ]: 8 : delete[] _smthCurve;
207 : 8 : _smthCurve = NULL;
208 : 8 : size1 = _my_gsets[2].size();
209 [ + + ]: 20 : for( i = 0; i < size1; i++ )
210 [ + - ]: 12 : delete _smthFace[i];
211 [ + - ]: 8 : delete[] _smthFace;
212 : 8 : _smthFace = NULL;
213 : : //_smooth = false;
214 : : }
215 : :
216 [ + + ]: 72 : for( int j = 0; j < 5; j++ )
217 : 60 : _my_gsets[j].clear();
218 [ + + ][ + - ]: 12 : if( _t_created ) delete _my_geomTopoTool;
219 : 12 : _my_geomTopoTool = NULL;
220 : 12 : _t_created = false;
221 : 12 : }
222 : :
223 : 11 : ErrorCode FBEngine::Init()
224 : : {
225 [ + + ]: 11 : if( !_initialized )
226 : : {
227 [ - + ]: 10 : if( !_my_geomTopoTool ) return MB_FAILURE;
228 : :
229 : 10 : ErrorCode rval = _my_geomTopoTool->find_geomsets( _my_gsets );
230 [ - + ]: 10 : assert( rval == MB_SUCCESS );
231 [ - + ]: 10 : if( MB_SUCCESS != rval ) { return rval; }
232 : :
233 : 10 : rval = split_quads();
234 [ - + ]: 10 : assert( rval == MB_SUCCESS );
235 : :
236 : 10 : rval = _my_geomTopoTool->construct_obb_trees();
237 [ - + ]: 10 : assert( rval == MB_SUCCESS );
238 : :
239 [ + + ]: 10 : if( _smooth ) rval = initializeSmoothing();
240 [ - + ]: 10 : assert( rval == MB_SUCCESS );
241 : :
242 : 10 : _initialized = true;
243 : : }
244 : 11 : return MB_SUCCESS;
245 : : }
246 : 8 : ErrorCode FBEngine::initializeSmoothing()
247 : : {
248 : : //
249 : : /*ErrorCode rval = Init();
250 : : MBERRORR(rval, "failed initialize");*/
251 : : // first of all, we need to retrieve all the surfaces from the (root) set
252 : : // in icesheet_test we use iGeom, but maybe that is a stretch
253 : : // get directly the sets with geom dim 2, and from there create the SmoothFace
254 : : Tag geom_tag;
255 [ + - ]: 8 : ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
256 [ - + ][ # # ]: 8 : MBERRORR( rval, "can't get geom tag" );
[ # # ]
257 : :
258 [ + - ]: 8 : int numSurfaces = _my_gsets[2].size();
259 : : // SmoothFace ** smthFace = new SmoothFace *[numSurfaces];
260 [ + - ][ + - ]: 8 : _smthFace = new SmoothFace*[numSurfaces];
261 : :
262 : : // there should also be a map from surfaces to evaluators
263 : : // std::map<MBEntityHandle, SmoothFace*> mapSurfaces;
264 : :
265 : 8 : int i = 0;
266 [ + - ]: 8 : Range::iterator it;
267 [ + - ][ + - ]: 20 : for( it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it, i++ )
[ + - ][ + - ]
[ + + ]
268 : : {
269 [ + - ]: 12 : EntityHandle face = *it;
270 [ + - ][ + - ]: 12 : _smthFace[i] = new SmoothFace( MBI, face, _my_geomTopoTool ); // geom topo tool will be used for searching,
271 : : // among other things; also for senses in edge sets...
272 [ + - ]: 12 : _faces[face] = _smthFace[i];
273 : : }
274 : :
275 [ + - ]: 8 : int numCurves = _my_gsets[1].size(); // csets.size();
276 : : // SmoothCurve ** smthCurve = new SmoothCurve *[numCurves];
277 [ + - ][ + - ]: 8 : _smthCurve = new SmoothCurve*[numCurves];
278 : : // there should also be a map from surfaces to evaluators
279 : : // std::map<MBEntityHandle, SmoothCurve*> mapCurves;
280 : :
281 : 8 : i = 0;
282 [ + - ][ + - ]: 38 : for( it = _my_gsets[1].begin(); it != _my_gsets[1].end(); ++it, i++ )
[ + - ][ + - ]
[ + + ]
283 : : {
284 [ + - ]: 30 : EntityHandle curve = *it;
285 [ + - ][ + - ]: 30 : _smthCurve[i] = new SmoothCurve( MBI, curve, _my_geomTopoTool );
286 [ + - ]: 30 : _edges[curve] = _smthCurve[i];
287 : : }
288 : :
289 [ + + ]: 20 : for( i = 0; i < numSurfaces; i++ )
290 : : {
291 [ + - ]: 12 : _smthFace[i]->init_gradient(); // this will also retrieve the triangles in each surface
292 [ + - ]: 12 : _smthFace[i]->compute_tangents_for_each_edge(); // this one will consider all edges
293 : : // internal, so the
294 : : // tangents are all in the direction of the edge; a little bit of waste, as we store
295 : : // one tangent for each edge node , even though they are equal here...
296 : : // no loops are considered
297 : : }
298 : :
299 : : // this will be used to mark boundary edges, so for them the control points are computed earlier
300 : 8 : unsigned char value = 0; // default value is "not used"=0 for the tag
301 : : // unsigned char def_data_bit = 1;// valid by default
302 : : // rval = mb->tag_create("valid", 1, MB_TAG_BIT, validTag, &def_data_bit);
303 : : Tag markTag;
304 : : rval = MBI->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, markTag, MB_TAG_EXCL | MB_TAG_BIT,
305 [ + - ]: 8 : &value ); // default value : 0 = not computed yet
306 : : // each feature edge will need to have a way to retrieve at every moment the surfaces it belongs
307 : : // to from edge sets, using the sense tag, we can get faces, and from each face, using the map,
308 : : // we can get the SmoothFace (surface evaluator), that has everything, including the normals!!!
309 [ - + ]: 8 : assert( rval == MB_SUCCESS );
310 : :
311 : : // create the tag also for control points on the edges
312 : 8 : double defCtrlPoints[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
313 : : Tag edgeCtrlTag;
314 : : rval = MBI->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, edgeCtrlTag, MB_TAG_DENSE | MB_TAG_CREAT,
315 [ + - ]: 8 : &defCtrlPoints );
316 [ - + ]: 8 : assert( rval == MB_SUCCESS );
317 : :
318 : : Tag facetCtrlTag;
319 : 8 : double defControls[18] = { 0. };
320 : : rval = MBI->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, facetCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
321 [ + - ]: 8 : &defControls );
322 [ - + ]: 8 : assert( rval == MB_SUCCESS );
323 : :
324 : : Tag facetEdgeCtrlTag;
325 : 8 : double defControls2[27] = { 0. }; // corresponding to 9 control points on edges, in order
326 : : // from edge 0, 1, 2 ( 1-2, 2-0, 0-1 )
327 : : rval = MBI->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, facetEdgeCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
328 [ + - ]: 8 : &defControls2 );
329 [ - + ]: 8 : assert( rval == MB_SUCCESS );
330 : : // if the
331 : 8 : double min_dot = -1.0; // depends on _angle, but now we ignore it, for the time being
332 [ + + ]: 38 : for( i = 0; i < numCurves; i++ )
333 : : {
334 [ + - ]: 30 : _smthCurve[i]->compute_tangents_for_each_edge(); // do we need surfaces now? or just the chains?
335 : : // the computed edges will be marked accordingly; later one, only internal edges to surfaces
336 : : // are left
337 [ + - ]: 30 : _smthCurve[i]->compute_control_points_on_boundary_edges( min_dot, _faces, edgeCtrlTag, markTag );
338 : : }
339 : :
340 : : // when done with boundary edges, compute the control points on all edges in the surfaces
341 : :
342 [ + + ]: 20 : for( i = 0; i < numSurfaces; i++ )
343 : : {
344 : : // also pass the tags for
345 [ + - ]: 12 : _smthFace[i]->compute_control_points_on_edges( min_dot, edgeCtrlTag, markTag );
346 : : }
347 : :
348 : : // now we should be able to compute the control points for the facets
349 : :
350 [ + + ]: 20 : for( i = 0; i < numSurfaces; i++ )
351 : : {
352 : : // also pass the tags for edge and facet control points
353 [ + - ]: 12 : _smthFace[i]->compute_internal_control_points_on_facets( min_dot, facetCtrlTag, facetEdgeCtrlTag );
354 : : }
355 : : // we will need to compute the tangents for all edges in the model
356 : : // they will be needed for control points for each edge
357 : : // the boundary edges and the feature edges are more complicated
358 : : // the boundary edges need to consider their loops, but feature edges need to consider loops and
359 : : // the normals on each connected surface
360 : :
361 : : // some control points
362 : : if( Debug_surf_eval )
363 : : for( i = 0; i < numSurfaces; i++ )
364 : : _smthFace[i]->DumpModelControlPoints();
365 : :
366 : 8 : return MB_SUCCESS;
367 : : }
368 : :
369 : : // clean up the smooth tags data if created, so the files will be smaller
370 : : // if saved
371 : : // also, recompute the tags if topology is modified
372 : 4 : void FBEngine::delete_smooth_tags()
373 : : {
374 : : // get all tags from database that are created for smooth data, and
375 : : // delete them; it will delete all data associated with them
376 : : // first tags from faces, edges:
377 [ + - ]: 4 : std::vector< Tag > smoothTags;
378 [ + - ]: 4 : int size1 = (int)_my_gsets[2].size();
379 : :
380 [ + + ]: 10 : for( int i = 0; i < size1; i++ )
381 : : {
382 : : // these 2 will append gradient tag and plane tag
383 [ + - ]: 6 : _smthFace[i]->append_smooth_tags( smoothTags );
384 : : }
385 : : // then , get other tags:
386 : : // "TANGENTS", "MARKER", "CONTROLEDGE", "CONTROLFACE", "CONTROLEDGEFACE"
387 : : Tag tag_handle;
388 [ + - ]: 4 : ErrorCode rval = _mbImpl->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tag_handle );
389 [ + - ][ + - ]: 4 : if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
390 : :
391 [ + - ]: 4 : rval = _mbImpl->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, tag_handle );
392 [ + - ][ + - ]: 4 : if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
393 : :
394 [ + - ]: 4 : rval = _mbImpl->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, tag_handle );
395 [ + - ][ + - ]: 4 : if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
396 : :
397 [ + - ]: 4 : rval = _mbImpl->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, tag_handle );
398 [ + - ][ + - ]: 4 : if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
399 : :
400 [ + - ]: 4 : rval = _mbImpl->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, tag_handle );
401 [ + - ][ + - ]: 4 : if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
402 : :
403 : : // a lot of tags, delete them
404 [ + + ]: 36 : for( unsigned int k = 0; k < smoothTags.size(); k++ )
405 : : {
406 : : // could be a lot of data
407 [ + - ][ + - ]: 32 : _mbImpl->tag_delete( smoothTags[k] );
408 : 4 : }
409 : 4 : }
410 : : /*
411 : : #define COPY_RANGE(r, vec) { \
412 : : EntityHandle *tmp_ptr = reinterpret_cast<EntityHandle*>(vec); \
413 : : std::copy(r.begin(), r.end(), tmp_ptr);}
414 : : */
415 : :
416 : : /*static inline void
417 : : ProcessError(const char* desc);*/
418 : :
419 : 37 : ErrorCode FBEngine::getRootSet( EntityHandle* root_set )
420 : : {
421 : 37 : *root_set = _my_geomTopoTool->get_root_model_set();
422 : 37 : return MB_SUCCESS;
423 : : }
424 : :
425 : 3 : ErrorCode FBEngine::getNumEntSets( EntityHandle set, int num_hops, int* all_sets )
426 : : {
427 : 3 : ErrorCode rval = MBI->num_contained_meshsets( set, all_sets, num_hops + 1 );
428 : 3 : return rval;
429 : : }
430 : :
431 : 18 : ErrorCode FBEngine::createEntSet( int isList, EntityHandle* pSet )
432 : : {
433 : : ErrorCode rval;
434 : :
435 [ + - ]: 18 : if( isList )
436 : 18 : rval = MBI->create_meshset( MESHSET_ORDERED, *pSet );
437 : : else
438 : 0 : rval = MBI->create_meshset( MESHSET_SET, *pSet );
439 : :
440 : 18 : return rval;
441 : : }
442 : :
443 : 118 : ErrorCode FBEngine::getEntities( EntityHandle set_handle, int entity_type, Range& gentities )
444 : : {
445 : : int i;
446 [ + - ][ - + ]: 118 : if( 0 > entity_type || 4 < entity_type ) { return MB_FAILURE; }
447 [ + + ]: 118 : else if( entity_type < 4 )
448 : : { // 4 means all entities
449 [ + - ][ + - ]: 111 : gentities = _my_geomTopoTool->geoRanges()[entity_type]; // all from root set!
450 : : }
451 : : else
452 : : {
453 [ + - ]: 7 : gentities.clear();
454 [ + + ]: 35 : for( i = 0; i < 4; i++ )
455 : : {
456 [ + - ][ + - ]: 28 : gentities.merge( _my_geomTopoTool->geoRanges()[i] );
457 : : }
458 : : }
459 [ + - ]: 118 : Range sets;
460 : : // see now if they are in the set passed as input or not
461 [ + - ]: 118 : ErrorCode rval = MBI->get_entities_by_type( set_handle, MBENTITYSET, sets );
462 [ - + ][ # # ]: 118 : MBERRORR( rval, "can't get sets in the initial set" );
[ # # ]
463 [ + - ][ + - ]: 118 : gentities = intersect( gentities, sets );
464 : :
465 : 118 : return MB_SUCCESS;
466 : : }
467 : :
468 : 18 : ErrorCode FBEngine::addEntArrToSet( Range entities, EntityHandle set )
469 : : {
470 : 18 : return MBI->add_entities( set, entities );
471 : : }
472 : :
473 : 12 : ErrorCode FBEngine::addEntSet( EntityHandle entity_set_to_add, EntityHandle entity_set_handle )
474 : : {
475 : 12 : return MBI->add_entities( entity_set_handle, &entity_set_to_add, 1 );
476 : : }
477 : :
478 : 39 : ErrorCode FBEngine::getNumOfType( EntityHandle set, int ent_type, int* pNum )
479 : : {
480 [ + - ][ - + ]: 39 : if( 0 > ent_type || 4 < ent_type )
481 : : {
482 [ # # ]: 0 : std::cout << "Invalid type\n";
483 : 0 : return MB_FAILURE;
484 : : }
485 : : // get sets of geom dimension tag from here, and intersect with the gentities from geo
486 : : // ranges
487 : :
488 : : // get the geom dimensions sets in the set (AKA gentities)
489 [ + - ]: 39 : Range geom_sets;
490 : : Tag geom_tag;
491 : : ErrorCode rval =
492 [ + - ]: 39 : _mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
493 [ - + ][ # # ]: 39 : MBERRORR( rval, "Failed to get geom tag." );
[ # # ]
494 [ + - ]: 39 : rval = _mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, &geom_tag, NULL, 1, geom_sets, Interface::UNION );
495 [ - + ][ # # ]: 39 : MBERRORR( rval, "Failed to get gentities from set" );
[ # # ]
496 : :
497 [ - + ]: 39 : if( ent_type == 4 )
498 : : {
499 : 0 : *pNum = 0;
500 [ # # ]: 0 : for( int k = 0; k <= 3; k++ )
501 : : {
502 [ # # ][ # # ]: 0 : Range gEntsOfTypeK = intersect( geom_sets, _my_geomTopoTool->geoRanges()[k] );
503 [ # # ]: 0 : *pNum += (int)gEntsOfTypeK.size();
504 : 0 : }
505 : : }
506 : : else
507 : : {
508 [ + - ][ + - ]: 39 : Range gEntsOfType = intersect( geom_sets, _my_geomTopoTool->geoRanges()[ent_type] );
509 [ + - ]: 39 : *pNum = (int)gEntsOfType.size();
510 : : }
511 : : // we do not really check if it is in the set or not;
512 : : // _my_gsets[i].find(gent) != _my_gsets[i].end()
513 : 39 : return MB_SUCCESS;
514 : : }
515 : :
516 : 223 : ErrorCode FBEngine::getEntType( EntityHandle gent, int* type )
517 : : {
518 [ + - ]: 343 : for( int i = 0; i < 4; i++ )
519 : : {
520 [ + - ][ + - ]: 343 : if( _my_geomTopoTool->geoRanges()[i].find( gent ) != _my_geomTopoTool->geoRanges()[i].end() )
[ + - ][ + + ]
521 : : {
522 : 223 : *type = i;
523 : 223 : return MB_SUCCESS;
524 : : }
525 : : }
526 : 0 : *type = -1; // failure
527 : 223 : return MB_FAILURE;
528 : : }
529 : 22 : ErrorCode FBEngine::getEntBoundBox( EntityHandle gent, double* min_x, double* min_y, double* min_z, double* max_x,
530 : : double* max_y, double* max_z )
531 : : {
532 : : ErrorCode rval;
533 : : int type;
534 [ + - ]: 22 : rval = getEntType( gent, &type );
535 [ - + ][ # # ]: 22 : MBERRORR( rval, "Failed to get entity type." );
[ # # ]
536 : :
537 [ + + ]: 22 : if( type == 0 )
538 : : {
539 [ + - ]: 11 : rval = getVtxCoord( gent, min_x, min_y, min_z );
540 [ - + ][ # # ]: 11 : MBERRORR( rval, "Failed to get vertex coordinates." );
[ # # ]
541 : 11 : max_x = min_x;
542 : 11 : max_y = min_y;
543 : 11 : max_z = min_z;
544 : : }
545 [ - + ]: 11 : else if( type == 1 )
546 : : {
547 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, "iGeom_getEntBoundBox is not supported for Edge entity type." );
548 : : }
549 [ - + ][ # # ]: 11 : else if( type == 2 || type == 3 )
550 : : {
551 : :
552 : : EntityHandle root;
553 [ + - ][ + - ]: 44 : CartVect center, axis[3];
[ + + ]
554 [ + - ]: 11 : rval = _my_geomTopoTool->get_root( gent, root );
555 [ - + ][ # # ]: 11 : MBERRORR( rval, "Failed to get tree root in iGeom_getEntBoundBox." );
[ # # ]
556 : : rval = _my_geomTopoTool->obb_tree()->box( root, center.array(), axis[0].array(), axis[1].array(),
557 [ + - ][ + - ]: 11 : axis[2].array() );
[ + - ][ + - ]
[ + - ][ + - ]
558 [ - + ][ # # ]: 11 : MBERRORR( rval, "Failed to get closest point in iGeom_getEntBoundBox." );
[ # # ]
559 : :
560 [ + - ][ + + ]: 44 : CartVect absv[3];
561 [ + + ]: 44 : for( int i = 0; i < 3; i++ )
562 : : {
563 [ + - ][ + - ]: 33 : absv[i] = CartVect( fabs( axis[i][0] ), fabs( axis[i][1] ), fabs( axis[i][2] ) );
[ + - ][ + - ]
564 : : }
565 [ + - ][ + - ]: 11 : CartVect min, max;
566 [ + - ][ + - ]: 11 : min = center - absv[0] - absv[1] - absv[2];
[ + - ]
567 [ + - ][ + - ]: 11 : max = center + absv[0] + absv[1] + absv[2];
[ + - ]
568 [ + - ]: 11 : *min_x = min[0];
569 [ + - ]: 11 : *min_y = min[1];
570 [ + - ]: 11 : *min_z = min[2];
571 [ + - ]: 11 : *max_x = max[0];
572 [ + - ]: 11 : *max_y = max[1];
573 [ + - ]: 11 : *max_z = max[2];
574 : : }
575 : : else
576 : 0 : return MB_FAILURE;
577 : :
578 : 22 : return MB_SUCCESS;
579 : : }
580 : 51 : ErrorCode FBEngine::getEntClosestPt( EntityHandle this_gent, double near_x, double near_y, double near_z, double* on_x,
581 : : double* on_y, double* on_z )
582 : : {
583 : : ErrorCode rval;
584 : : int type;
585 [ + - ]: 51 : rval = getEntType( this_gent, &type );
586 [ - + ][ # # ]: 51 : MBERRORR( rval, "Failed to get entity type." );
[ # # ]
587 : :
588 [ + + ]: 51 : if( type == 0 )
589 : : {
590 [ + - ]: 29 : rval = getVtxCoord( this_gent, on_x, on_y, on_z );
591 [ - + ][ # # ]: 29 : MBERRORR( rval, "Failed to get vertex coordinates." );
[ # # ]
592 : : }
593 [ + + ][ + + ]: 22 : else if( _smooth && type == 1 )
594 : : {
595 : 12 : *on_x = near_x;
596 : 12 : *on_y = near_y;
597 : 12 : *on_z = near_z;
598 [ + - ]: 12 : SmoothCurve* smthcurve = _edges[this_gent];
599 : : // call the new method from smooth edge
600 [ + - ]: 12 : smthcurve->move_to_curve( *on_x, *on_y, *on_z );
601 : : }
602 [ - + ][ # # ]: 10 : else if( type == 2 || type == 3 )
603 : : {
604 : 10 : double point[3] = { near_x, near_y, near_z };
605 : : double point_out[3];
606 : : EntityHandle root, facet_out;
607 [ + + ][ + - ]: 10 : if( _smooth && 2 == type )
608 : : {
609 [ + - ]: 8 : SmoothFace* smthFace = _faces[this_gent];
610 : 8 : *on_x = near_x;
611 : 8 : *on_y = near_y;
612 : 8 : *on_z = near_z;
613 [ + - ]: 8 : smthFace->move_to_surface( *on_x, *on_y, *on_z );
614 : : }
615 : : else
616 : : {
617 [ + - ]: 2 : rval = _my_geomTopoTool->get_root( this_gent, root );
618 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to get tree root in iGeom_getEntClosestPt." );
[ # # ]
619 [ + - ][ + - ]: 2 : rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
620 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to get closest point in iGeom_getEntClosestPt." );
[ # # ]
621 : :
622 : 2 : *on_x = point_out[0];
623 : 2 : *on_y = point_out[1];
624 : 2 : *on_z = point_out[2];
625 : 10 : }
626 : : }
627 : : else
628 : 0 : return MB_TYPE_OUT_OF_RANGE;
629 : :
630 : 51 : return MB_SUCCESS;
631 : : }
632 : :
633 : 58 : ErrorCode FBEngine::getVtxCoord( EntityHandle vertex_handle, double* x0, double* y0, double* z0 )
634 : : {
635 : : int type;
636 [ + - ]: 58 : ErrorCode rval = getEntType( vertex_handle, &type );
637 [ - + ][ # # ]: 58 : MBERRORR( rval, "Failed to get entity type in getVtxCoord." );
[ # # ]
638 : :
639 [ - + ][ # # ]: 58 : if( type != 0 ) { MBERRORR( MB_FAILURE, "Entity is not a vertex type." ); }
[ # # ]
640 : :
641 [ + - ]: 58 : Range entities;
642 [ + - ]: 58 : rval = MBI->get_entities_by_type( vertex_handle, MBVERTEX, entities );
643 [ - + ][ # # ]: 58 : MBERRORR( rval, "can't get nodes in vertex set." );
[ # # ]
644 : :
645 [ + - ][ - + ]: 58 : if( entities.size() != 1 ) { MBERRORR( MB_FAILURE, "Vertex has multiple points." ); }
[ # # ][ # # ]
646 : : double coords[3];
647 [ + - ]: 58 : EntityHandle node = entities[0];
648 [ + - ]: 58 : rval = MBI->get_coords( &node, 1, coords );
649 [ - + ][ # # ]: 58 : MBERRORR( rval, "can't get coordinates." );
[ # # ]
650 : 58 : *x0 = coords[0];
651 : 58 : *y0 = coords[1];
652 : 58 : *z0 = coords[2];
653 : :
654 : 58 : return MB_SUCCESS;
655 : : }
656 : :
657 : 3 : ErrorCode FBEngine::gsubtract( EntityHandle entity_set_1, EntityHandle entity_set_2, EntityHandle result_entity_set )
658 : : {
659 : : /*result_entity_set = subtract(entity_set_1, entity_set_2);*/
660 [ + - ][ + - ]: 6 : Range ents1, ents2;
661 [ + - ]: 3 : ErrorCode rval = MBI->get_entities_by_type( entity_set_1, MBENTITYSET, ents1 );
662 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't get entities from set 1." );
[ # # ]
663 : :
664 [ + - ]: 3 : rval = MBI->get_entities_by_type( entity_set_2, MBENTITYSET, ents2 );
665 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't get entities from set 2." );
[ # # ]
666 : :
667 [ + - ][ + - ]: 3 : ents1 = subtract( ents1, ents2 );
668 [ + - ]: 3 : rval = MBI->clear_meshset( &result_entity_set, 1 );
669 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't empty set." );
[ # # ]
670 : :
671 [ + - ]: 3 : rval = MBI->add_entities( result_entity_set, ents1 );
672 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't add result to set." );
[ # # ]
673 : :
674 : 6 : return rval;
675 : : }
676 : :
677 : 8 : ErrorCode FBEngine::getEntNrmlXYZ( EntityHandle entity_handle, double x, double y, double z, double* nrml_i,
678 : : double* nrml_j, double* nrml_k )
679 : : {
680 : : // just do for surface and volume
681 : : int type;
682 [ + - ]: 8 : ErrorCode rval = getEntType( entity_handle, &type );
683 [ - + ][ # # ]: 8 : MBERRORR( rval, "Failed to get entity type in iGeom_getEntNrmlXYZ." );
[ # # ]
684 : :
685 [ - + ][ # # ]: 8 : if( type != 2 && type != 3 )
686 [ # # ][ # # ]: 0 : { MBERRORR( MB_FAILURE, "Entities passed into gentityNormal must be face or volume." ); }
687 : :
688 [ + - ][ + - ]: 8 : if( _smooth && 2 == type )
689 : : {
690 [ + - ]: 8 : SmoothFace* smthFace = _faces[entity_handle];
691 : : //*on_x = near_x; *on_y = near_y; *on_z = near_z;
692 [ + - ]: 8 : smthFace->normal_at( x, y, z, *nrml_i, *nrml_j, *nrml_k );
693 : : }
694 : : else
695 : : {
696 : : // get closest location and facet
697 : 0 : double point[3] = { x, y, z };
698 : : double point_out[3];
699 : : EntityHandle root, facet_out;
700 [ # # ]: 0 : _my_geomTopoTool->get_root( entity_handle, root );
701 [ # # ][ # # ]: 0 : rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
702 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to get closest location in iGeom_getEntNrmlXYZ." );
[ # # ]
703 : :
704 : : // get facet normal
705 : : const EntityHandle* conn;
706 : : int len;
707 [ # # ][ # # ]: 0 : CartVect coords[3], normal;
[ # # ]
708 [ # # ]: 0 : rval = MBI->get_connectivity( facet_out, conn, len );
709 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to get triangle connectivity in iGeom_getEntNrmlXYZ." );
[ # # ]
710 [ # # ][ # # ]: 0 : if( len != 3 ) MBERRORR( MB_FAILURE, " not a triangle, error " );
[ # # ]
711 : :
712 [ # # ][ # # ]: 0 : rval = MBI->get_coords( conn, len, coords[0].array() );
713 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to get triangle coordinates in iGeom_getEntNrmlXYZ." );
[ # # ]
714 : :
715 [ # # ]: 0 : coords[1] -= coords[0];
716 [ # # ]: 0 : coords[2] -= coords[0];
717 [ # # ]: 0 : normal = coords[1] * coords[2];
718 [ # # ]: 0 : normal.normalize();
719 [ # # ]: 0 : *nrml_i = normal[0];
720 [ # # ]: 0 : *nrml_j = normal[1];
721 [ # # ]: 0 : *nrml_k = normal[2];
722 : : }
723 : 8 : return MB_SUCCESS;
724 : : }
725 : :
726 : 5 : ErrorCode FBEngine::getPntRayIntsct( double x, double y, double z, double dir_x, double dir_y, double dir_z,
727 : : std::vector< EntityHandle >& intersect_entity_handles,
728 : : /* int storage_order,*/
729 : : std::vector< double >& intersect_coords, std::vector< double >& param_coords )
730 : : {
731 : : // this is pretty cool
732 : : // we will return only surfaces (gentities )
733 : : //
734 : : ErrorCode rval;
735 : :
736 [ + - ]: 5 : unsigned int numfaces = _my_gsets[2].size();
737 : : // do ray fire
738 : 5 : const double point[] = { x, y, z };
739 : 5 : const double dir[] = { dir_x, dir_y, dir_z };
740 [ + - ]: 5 : CartVect P( point );
741 [ + - ]: 5 : CartVect V( dir );
742 : :
743 : : // std::vector<double> distances;
744 [ + - ]: 5 : std::vector< EntityHandle > facets;
745 : : // std::vector<EntityHandle> sets;
746 : : unsigned int i;
747 [ + + ]: 13 : for( i = 0; i < numfaces; i++ )
748 : : {
749 [ + - ]: 8 : EntityHandle face = _my_gsets[2][i];
750 : : EntityHandle rootForFace;
751 [ + - ]: 8 : rval = _my_geomTopoTool->get_root( face, rootForFace );
752 [ - + ][ # # ]: 8 : MBERRORR( rval, "Failed to get root of face." );
[ # # ]
753 [ + - ]: 8 : std::vector< double > distances_out;
754 [ + - ][ + - ]: 16 : std::vector< EntityHandle > sets_out;
755 [ + - ][ + - ]: 16 : std::vector< EntityHandle > facets_out;
756 : : rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
757 [ + - ][ + - ]: 8 : tolerance, point, dir );
758 : : unsigned int j;
759 [ + + ]: 12 : for( j = 0; j < distances_out.size(); j++ )
760 [ + - ][ + - ]: 4 : param_coords.push_back( distances_out[j] );
761 [ + + ]: 12 : for( j = 0; j < sets_out.size(); j++ )
762 [ + - ][ + - ]: 4 : intersect_entity_handles.push_back( sets_out[j] );
763 [ + + ]: 12 : for( j = 0; j < facets_out.size(); j++ )
764 [ + - ][ + - ]: 4 : facets.push_back( facets_out[j] );
765 : :
766 [ - + ][ # # ]: 8 : MBERRORR( rval, "Failed to get ray intersections." );
[ # # ][ + - ]
767 : 8 : }
768 : : // facets.size == distances.size()!!
769 [ + + ]: 9 : for( i = 0; i < param_coords.size(); i++ )
770 : : {
771 [ + - ][ + - ]: 4 : CartVect intx = P + param_coords[i] * V;
[ + - ]
772 [ + + ]: 16 : for( int j = 0; j < 3; j++ )
773 [ + - ][ + - ]: 12 : intersect_coords.push_back( intx[j] );
774 : : }
775 [ + - ]: 5 : if( _smooth )
776 : : {
777 : : // correct the intersection point and the distance for smooth surfaces
778 [ + + ]: 9 : for( i = 0; i < intersect_entity_handles.size(); i++ )
779 : : {
780 : : // EntityHandle geoSet = MBH_cast(sets[i]);
781 [ + - ][ + - ]: 4 : SmoothFace* sFace = _faces[intersect_entity_handles[i]];
782 : : // correct coordinates and distance from point
783 : : /*moab::ErrorCode ray_intersection_correct(moab::EntityHandle facet, // (IN) the facet
784 : : where the patch is defined moab::CartVect &pt, // (IN) shoot from moab::CartVect &ray,
785 : : // (IN) ray direction moab::CartVect &eval_pt, // (INOUT) The intersection point double
786 : : & distance, // (IN OUT) the new distance bool &outside);*/
787 [ + - ][ + - ]: 4 : CartVect pos( &( intersect_coords[3 * i] ) );
788 [ + - ]: 4 : double dist = param_coords[i];
789 : 4 : bool outside = false;
790 [ + - ][ + - ]: 4 : rval = sFace->ray_intersection_correct( facets[i], P, V, pos, dist, outside );
791 [ - + ][ # # ]: 4 : MBERRORR( rval, "Failed to get better point on ray." );
[ # # ]
792 [ + - ]: 4 : param_coords[i] = dist;
793 : :
794 [ + + ]: 16 : for( int j = 0; j < 3; j++ )
795 [ + - ][ + - ]: 12 : intersect_coords[3 * i + j] = pos[j];
796 : : }
797 : : }
798 : 5 : return MB_SUCCESS;
799 : : }
800 : :
801 : 135 : ErrorCode FBEngine::getAdjacentEntities( const EntityHandle from, const int to_dim, Range& adjs )
802 : : {
803 : 135 : int this_dim = -1;
804 [ + - ]: 224 : for( int i = 0; i < 4; i++ )
805 : : {
806 [ + - ][ + - ]: 224 : if( _my_geomTopoTool->geoRanges()[i].find( from ) != _my_geomTopoTool->geoRanges()[i].end() )
[ + - ][ + + ]
807 : : {
808 : 135 : this_dim = i;
809 : 135 : break;
810 : : }
811 : : }
812 : :
813 : : // check target dimension
814 [ - + ]: 135 : if( -1 == this_dim )
815 : : {
816 : : // ProcessError(iBase_FAILURE, "Entity not a geometry entity.");
817 : 0 : return MB_FAILURE;
818 : : }
819 [ + - ][ - + ]: 135 : else if( 0 > to_dim || 3 < to_dim )
820 : : {
821 : : // ProcessError(iBase_FAILURE, "To dimension must be between 0 and 3.");
822 : 0 : return MB_FAILURE;
823 : : }
824 [ - + ]: 135 : else if( to_dim == this_dim )
825 : : {
826 : : // ProcessError(iBase_FAILURE,
827 : : // "To dimension must be different from entity dimension.");
828 : 0 : return MB_FAILURE;
829 : : }
830 : :
831 : 135 : ErrorCode rval = MB_SUCCESS;
832 : 135 : adjs.clear();
833 [ + + ]: 135 : if( to_dim > this_dim )
834 : : {
835 : 92 : int diffDim = to_dim - this_dim;
836 : 92 : rval = MBI->get_parent_meshsets( from, adjs, diffDim );
837 [ - + ]: 92 : if( MB_SUCCESS != rval ) return rval;
838 [ + + ]: 92 : if( diffDim > 1 )
839 : : {
840 : : // subtract the parents that come with diffDim-1 hops
841 [ + - ]: 24 : Range extra;
842 [ + - ]: 24 : rval = MBI->get_parent_meshsets( from, extra, diffDim - 1 );
843 [ - + ]: 24 : if( MB_SUCCESS != rval ) return rval;
844 [ + - ][ + - ]: 92 : adjs = subtract( adjs, extra );
[ + - ]
845 : : }
846 : : }
847 : : else
848 : : {
849 : 43 : int diffDim = this_dim - to_dim;
850 : 43 : rval = MBI->get_child_meshsets( from, adjs, diffDim );
851 [ - + ]: 43 : if( MB_SUCCESS != rval ) return rval;
852 [ + + ]: 43 : if( diffDim > 1 )
853 : : {
854 : : // subtract the children that come with diffDim-1 hops
855 [ + - ]: 6 : Range extra;
856 [ + - ]: 6 : rval = MBI->get_child_meshsets( from, extra, diffDim - 1 );
857 [ - + ]: 6 : if( MB_SUCCESS != rval ) return rval;
858 [ + - ][ + - ]: 6 : adjs = subtract( adjs, extra );
[ + - ]
859 : : }
860 : : }
861 : :
862 : 135 : return rval;
863 : : }
864 : :
865 : : // so far, this one is
866 : : // used only for __MKModelEntityGeo tag
867 : :
868 : 0 : ErrorCode FBEngine::createTag( const char* tag_name, int tag_size, int tag_type, Tag& tag_handle_out )
869 : : {
870 : : // this is copied from iMesh_MOAB.cpp; different name to not have trouble
871 : : // with it
872 : : // also, we do not want to depend on iMesh.h...
873 : : // iMesh is more complicated, because of the options passed
874 : :
875 : : DataType mb_data_type_table2[] = { MB_TYPE_OPAQUE, MB_TYPE_INTEGER, MB_TYPE_DOUBLE, MB_TYPE_HANDLE,
876 : 0 : MB_TYPE_HANDLE };
877 : 0 : moab::TagType storage = MB_TAG_SPARSE;
878 : : ErrorCode result;
879 : :
880 : : result =
881 [ # # ]: 0 : MBI->tag_get_handle( tag_name, tag_size, mb_data_type_table2[tag_type], tag_handle_out, storage | MB_TAG_EXCL );
882 : :
883 [ # # ]: 0 : if( MB_SUCCESS != result )
884 : : {
885 [ # # ]: 0 : std::string msg( "iMesh_createTag: " );
886 [ # # ]: 0 : if( MB_ALREADY_ALLOCATED == result )
887 : : {
888 [ # # ]: 0 : msg += "Tag already exists with name: \"";
889 [ # # ]: 0 : msg += tag_name;
890 [ # # ][ # # ]: 0 : std::cout << msg << "\n";
891 : : }
892 : : else
893 : : {
894 [ # # ][ # # ]: 0 : std::cout << "Failed to create tag with name: " << tag_name << "\n";
[ # # ]
895 [ # # ]: 0 : return MB_FAILURE;
896 : 0 : }
897 : : }
898 : :
899 : : // end copy
900 : 0 : return MB_SUCCESS;
901 : : }
902 : :
903 : 0 : ErrorCode FBEngine::getArrData( const EntityHandle* entity_handles, int entity_handles_size, Tag tag_handle,
904 : : void* tag_values_out )
905 : : {
906 : : // responsibility of the user to have tag_values_out properly allocated
907 : : // only some types of Tags are possible (double, int, etc)
908 : 0 : return MBI->tag_get_data( tag_handle, entity_handles, entity_handles_size, tag_values_out );
909 : : }
910 : :
911 : 0 : ErrorCode FBEngine::setArrData( const EntityHandle* entity_handles, int entity_handles_size, Tag tag_handle,
912 : : const void* tag_values )
913 : : {
914 : : // responsibility of the user to have tag_values_out properly allocated
915 : : // only some types of Tags are possible (double, int, etc)
916 : 0 : return MBI->tag_set_data( tag_handle, entity_handles, entity_handles_size, tag_values );
917 : : }
918 : :
919 : 123 : ErrorCode FBEngine::getEntAdj( EntityHandle handle, int type_requested, Range& adjEnts )
920 : : {
921 : 123 : return getAdjacentEntities( handle, type_requested, adjEnts );
922 : : }
923 : :
924 : 0 : ErrorCode FBEngine::getEgFcSense( EntityHandle mbedge, EntityHandle mbface, int& sense_out )
925 : : {
926 : :
927 : : // this one is important, for establishing the orientation of the edges in faces
928 : : // use senses
929 [ # # ]: 0 : std::vector< EntityHandle > faces;
930 [ # # ]: 0 : std::vector< int > senses; // 0 is forward and 1 is backward
931 [ # # ]: 0 : ErrorCode rval = _my_geomTopoTool->get_senses( mbedge, faces, senses );
932 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
933 : :
934 [ # # ]: 0 : for( unsigned int i = 0; i < faces.size(); i++ )
935 : : {
936 [ # # ][ # # ]: 0 : if( faces[i] == mbface )
937 : : {
938 [ # # ]: 0 : sense_out = senses[i];
939 : 0 : return MB_SUCCESS;
940 : : }
941 : : }
942 : 0 : return MB_FAILURE;
943 : : }
944 : : // we assume the measures array was allocated correctly
945 : 0 : ErrorCode FBEngine::measure( const EntityHandle* moab_entities, int entities_size, double* measures )
946 : : {
947 : : ErrorCode rval;
948 [ # # ]: 0 : for( int i = 0; i < entities_size; i++ )
949 : : {
950 : 0 : measures[i] = 0.;
951 : :
952 : : int type;
953 : 0 : EntityHandle gset = moab_entities[i];
954 [ # # ]: 0 : rval = getEntType( gset, &type );
955 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
956 [ # # ]: 0 : if( type == 1 )
957 : : { // edge: get all edges part of the edge set
958 [ # # ]: 0 : Range entities;
959 [ # # ]: 0 : rval = MBI->get_entities_by_type( gset, MBEDGE, entities );
960 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
961 : :
962 [ # # ][ # # ]: 0 : for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
963 : : {
964 [ # # ]: 0 : EntityHandle edge = *it;
965 [ # # ][ # # ]: 0 : CartVect vv[2];
966 : 0 : const EntityHandle* conn2 = NULL;
967 : : int num_nodes;
968 [ # # ]: 0 : rval = MBI->get_connectivity( edge, conn2, num_nodes );
969 [ # # ][ # # ]: 0 : if( MB_SUCCESS != rval || num_nodes != 2 ) return MB_FAILURE;
970 [ # # ][ # # ]: 0 : rval = MBI->get_coords( conn2, 2, (double*)&( vv[0][0] ) );
971 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
972 : :
973 [ # # ]: 0 : vv[0] = vv[1] - vv[0];
974 [ # # ]: 0 : measures[i] += vv[0].length();
975 : 0 : }
976 : : }
977 [ # # ]: 0 : if( type == 2 )
978 : : { // surface
979 : : // get triangles in surface; TODO: quads!
980 [ # # ]: 0 : Range entities;
981 [ # # ]: 0 : rval = MBI->get_entities_by_type( gset, MBTRI, entities );
982 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
983 : :
984 [ # # ][ # # ]: 0 : for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
985 : : {
986 [ # # ]: 0 : EntityHandle tri = *it;
987 [ # # ][ # # ]: 0 : CartVect vv[3];
988 : 0 : const EntityHandle* conn3 = NULL;
989 : : int num_nodes;
990 [ # # ]: 0 : rval = MBI->get_connectivity( tri, conn3, num_nodes );
991 [ # # ][ # # ]: 0 : if( MB_SUCCESS != rval || num_nodes != 3 ) return MB_FAILURE;
992 [ # # ][ # # ]: 0 : rval = MBI->get_coords( conn3, 3, (double*)&( vv[0][0] ) );
993 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
994 : :
995 [ # # ]: 0 : vv[1] = vv[1] - vv[0];
996 [ # # ]: 0 : vv[2] = vv[2] - vv[0];
997 [ # # ]: 0 : vv[0] = vv[1] * vv[2];
998 [ # # ]: 0 : measures[i] += vv[0].length() / 2; // area of triangle
999 : 0 : }
1000 : : }
1001 : : }
1002 : 0 : return MB_SUCCESS;
1003 : : }
1004 : :
1005 : 0 : ErrorCode FBEngine::getEntNrmlSense( EntityHandle /*face*/, EntityHandle /*region*/, int& /*sense*/ )
1006 : : {
1007 : 0 : return MB_NOT_IMPLEMENTED; // not implemented
1008 : : }
1009 : :
1010 : 0 : ErrorCode FBEngine::getEgEvalXYZ( EntityHandle /*edge*/, double /*x*/, double /*y*/, double /*z*/, double& /*on_x*/,
1011 : : double& /*on_y*/, double& /*on_z*/, double& /*tngt_i*/, double& /*tngt_j*/,
1012 : : double& /*tngt_k*/, double& /*cvtr_i*/, double& /*cvtr_j*/, double& /*cvtr_k*/ )
1013 : : {
1014 : 0 : return MB_NOT_IMPLEMENTED; // not implemented
1015 : : }
1016 : 0 : ErrorCode FBEngine::getFcEvalXYZ( EntityHandle /*face*/, double /*x*/, double /*y*/, double /*z*/, double& /*on_x*/,
1017 : : double& /*on_y*/, double& /*on_z*/, double& /*nrml_i*/, double& /*nrml_j*/,
1018 : : double& /*nrml_k*/, double& /*cvtr1_i*/, double& /*cvtr1_j*/, double& /*cvtr1_k*/,
1019 : : double& /*cvtr2_i*/, double& /*cvtr2_j*/, double& /*cvtr2_k*/ )
1020 : : {
1021 : 0 : return MB_NOT_IMPLEMENTED; // not implemented
1022 : : }
1023 : :
1024 : 8 : ErrorCode FBEngine::getEgVtxSense( EntityHandle edge, EntityHandle vtx1, EntityHandle vtx2, int& sense )
1025 : : {
1026 : : // need to decide first or second vertex
1027 : : // important for moab
1028 : : int type;
1029 : :
1030 : : EntityHandle v1, v2;
1031 [ + - ]: 8 : ErrorCode rval = getEntType( vtx1, &type );
1032 [ + - ][ - + ]: 8 : if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
1033 : : // edge: get one vertex as part of the vertex set
1034 [ + - ]: 8 : Range entities;
1035 [ + - ]: 8 : rval = MBI->get_entities_by_type( vtx1, MBVERTEX, entities );
1036 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
1037 [ + - ][ - + ]: 8 : if( entities.size() < 1 ) return MB_FAILURE;
1038 [ + - ]: 8 : v1 = entities[0]; // the first vertex
1039 [ + - ]: 8 : entities.clear();
1040 [ + - ]: 8 : rval = getEntType( vtx2, &type );
1041 [ + - ][ - + ]: 8 : if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
1042 [ + - ]: 8 : rval = MBI->get_entities_by_type( vtx2, MBVERTEX, entities );
1043 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
1044 [ + - ][ - + ]: 8 : if( entities.size() < 1 ) return MB_FAILURE;
1045 [ + - ]: 8 : v2 = entities[0]; // the first vertex
1046 [ + - ]: 8 : entities.clear();
1047 : : // now get the edges, and get the first node and the last node in sequence of edges
1048 : : // the order is important...
1049 : : // these are ordered sets !!
1050 [ + - ]: 16 : std::vector< EntityHandle > ents;
1051 [ + - ]: 8 : rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
1052 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
1053 [ - + ]: 8 : if( ents.size() < 1 ) return MB_FAILURE;
1054 : :
1055 : 8 : const EntityHandle* conn = NULL;
1056 : : int len;
1057 : : EntityHandle startNode, endNode;
1058 [ + - ][ + - ]: 8 : rval = MBI->get_connectivity( ents[0], conn, len );
1059 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
1060 : 8 : startNode = conn[0];
1061 [ + - ][ + - ]: 8 : rval = MBI->get_connectivity( ents[ents.size() - 1], conn, len );
1062 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
1063 : :
1064 : 8 : endNode = conn[1];
1065 : 8 : sense = 1; //
1066 [ - + ][ # # ]: 8 : if( ( startNode == endNode ) && ( v1 == startNode ) )
1067 : : {
1068 : 0 : sense = 0; // periodic
1069 : : }
1070 [ + + ][ + - ]: 8 : if( ( startNode == v1 ) && ( endNode == v2 ) )
1071 : : {
1072 : 3 : sense = 1; // forward
1073 : : }
1074 [ + + ][ + - ]: 8 : if( ( startNode == v2 ) && ( endNode == v1 ) )
1075 : : {
1076 : 5 : sense = -1; // reverse
1077 : : }
1078 : 16 : return MB_SUCCESS;
1079 : : }
1080 : :
1081 : 0 : ErrorCode FBEngine::getEntURange( EntityHandle edge, double& u_min, double& u_max )
1082 : : {
1083 : 0 : SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1084 : : // now, call smoothCurve methods
1085 : 0 : smoothCurve->get_param_range( u_min, u_max );
1086 : 0 : return MB_SUCCESS;
1087 : : }
1088 : :
1089 : 12 : ErrorCode FBEngine::getEntUtoXYZ( EntityHandle edge, double u, double& x, double& y, double& z )
1090 : : {
1091 : 12 : SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1092 : : // now, call smoothCurve methods
1093 : 12 : smoothCurve->position_from_u( u, x, y, z );
1094 : 12 : return MB_SUCCESS;
1095 : : }
1096 : :
1097 : 0 : ErrorCode FBEngine::getEntTgntU( EntityHandle edge, double u, double& i, double& j, double& k )
1098 : : {
1099 [ # # ]: 0 : SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1100 : : // now, call smoothCurve methods
1101 : : double tg[3];
1102 : : double x, y, z;
1103 [ # # ]: 0 : smoothCurve->position_from_u( u, x, y, z, tg );
1104 : 0 : i = tg[0];
1105 : 0 : j = tg[1];
1106 : 0 : k = tg[2];
1107 : 0 : return MB_SUCCESS;
1108 : : }
1109 : 0 : ErrorCode FBEngine::isEntAdj( EntityHandle entity1, EntityHandle entity2, bool& adjacent_out )
1110 : : {
1111 : : int type1, type2;
1112 [ # # ]: 0 : ErrorCode rval = getEntType( entity1, &type1 );
1113 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1114 [ # # ]: 0 : rval = getEntType( entity2, &type2 );
1115 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1116 : :
1117 [ # # ]: 0 : Range adjs;
1118 [ # # ]: 0 : if( type1 < type2 )
1119 : : {
1120 [ # # ]: 0 : rval = MBI->get_parent_meshsets( entity1, adjs, type2 - type1 );
1121 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get parent meshsets in iGeom_isEntAdj.");
1122 : : }
1123 : : else
1124 : : {
1125 : : // note: if they ave the same type, they will not be adjacent, in our definition
1126 [ # # ]: 0 : rval = MBI->get_child_meshsets( entity1, adjs, type1 - type2 );
1127 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get child meshsets in iGeom_isEntAdj.");
1128 : : }
1129 : :
1130 : : // adjacent_out = adjs.find(entity2) != _my_gsets[type2].end();
1131 : : // hmmm, possible bug here; is this called?
1132 [ # # ][ # # ]: 0 : adjacent_out = adjs.find( entity2 ) != adjs.end();
[ # # ]
1133 : :
1134 : 0 : return MB_SUCCESS;
1135 : : }
1136 : :
1137 : 4 : ErrorCode FBEngine::split_surface_with_direction( EntityHandle face, std::vector< double >& xyz, double* direction,
1138 : : int closed, double min_dot, EntityHandle& oNewFace )
1139 : : {
1140 : :
1141 : : // first of all, find all intersection points (piercing in the face along the direction)
1142 : : // assume it is robust; what if it is not sufficiently robust?
1143 : : // if the polyline is open, find the intersection with the boundary edges, of the
1144 : : // polyline extruded at ends
1145 : :
1146 : : ErrorCode rval;
1147 : :
1148 : : // then find the position
1149 : 4 : int numIniPoints = (int)xyz.size() / 3;
1150 [ + + ][ + - ]: 4 : if( ( closed && numIniPoints < 3 ) || ( !closed && numIniPoints < 2 ) )
[ + + ][ - + ]
1151 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, "not enough polyline points " );
1152 : : EntityHandle rootForFace;
1153 : :
1154 [ + - ]: 4 : rval = _my_geomTopoTool->get_root( face, rootForFace );
1155 [ - + ][ # # ]: 4 : MBERRORR( rval, "Failed to get root of face." );
[ # # ]
1156 : :
1157 : 4 : const double dir[] = { direction[0], direction[1], direction[2] };
1158 [ + - ]: 4 : std::vector< EntityHandle > nodes; // get the nodes closest to the ray traces of interest
1159 : :
1160 : : // these are nodes on the boundary of original face;
1161 : : // if the cutting line is not closed, the starting - ending vertices of the
1162 : : // polygonal line must come from this list
1163 : :
1164 [ + - ]: 8 : std::vector< CartVect > b_pos;
1165 [ + - ]: 8 : std::vector< EntityHandle > boundary_nodes;
1166 [ + - ]: 8 : std::vector< EntityHandle > splittingNodes;
1167 [ + - ]: 8 : Range boundary_mesh_edges;
1168 [ + + ]: 4 : if( !closed )
1169 : : {
1170 [ + - ]: 1 : rval = boundary_nodes_on_face( face, boundary_nodes );
1171 [ - + ][ # # ]: 1 : MBERRORR( rval, "Failed to get boundary nodes." );
[ # # ]
1172 [ + - ]: 1 : b_pos.resize( boundary_nodes.size() );
1173 [ + - ][ + - ]: 1 : rval = _mbImpl->get_coords( &( boundary_nodes[0] ), boundary_nodes.size(), (double*)( &b_pos[0][0] ) );
[ + - ][ + - ]
1174 [ - + ][ # # ]: 1 : MBERRORR( rval, "Failed to get coordinates for boundary nodes." );
[ # # ]
1175 [ + - ]: 1 : rval = boundary_mesh_edges_on_face( face, boundary_mesh_edges );
1176 [ - + ][ # # ]: 1 : MBERRORR( rval, "Failed to get mesh boundary edges for face." );
[ # # ]
1177 : : }
1178 : : //
1179 : 4 : int i = 0;
1180 [ + - ]: 4 : CartVect dirct( direction );
1181 [ + - ]: 4 : dirct.normalize(); // maybe an overkill?
1182 [ + + ]: 21 : for( ; i < numIniPoints; i++ )
1183 : : {
1184 : :
1185 [ + - ][ + - ]: 17 : const double point[] = { xyz[3 * i], xyz[3 * i + 1], xyz[3 * i + 2] }; // or even point( &(xyz[3*i]) ); //
[ + - ]
1186 [ + - ]: 17 : CartVect p1( point );
1187 [ + + ][ + + ]: 17 : if( !closed && ( ( 0 == i ) || ( numIniPoints - 1 == i ) ) )
[ + + ]
1188 : : {
1189 : :
1190 : : // find the intersection point between a plane and boundary mesh edges
1191 : : // this will be the closest point on the boundary of face
1192 : : /// the segment is the first or last segment in the polyline
1193 : 2 : int i1 = i + 1;
1194 [ + + ]: 2 : if( i == numIniPoints - 1 ) i1 = i - 1; // previous point if the last
1195 : : // the direction is from point to point1
1196 [ + - ][ + - ]: 2 : const double point1[] = { xyz[3 * i1], xyz[3 * i1 + 1], xyz[3 * i1 + 2] };
[ + - ]
1197 [ + - ]: 2 : CartVect p2( point1 );
1198 [ + - ][ + - ]: 2 : CartVect normPlane = ( p2 - p1 ) * dirct;
1199 [ + - ]: 2 : normPlane.normalize();
1200 : : //(roughly, from p1 to p2, perpendicular to dirct, in the "xy" plane
1201 : : // if the intx point is "outside" p1 - p2, skip if the intx point is closer to p2
1202 [ + - ]: 2 : CartVect perpDir = dirct * normPlane;
1203 [ + - ]: 2 : Range::iterator ite = boundary_mesh_edges.begin();
1204 : : // do a linear search for the best intersection point position (on a boundary edge)
1205 [ - + ]: 2 : if( debug_splits )
1206 : : {
1207 [ # # ][ # # ]: 0 : std::cout << " p1:" << p1 << "\n";
[ # # ]
1208 [ # # ][ # # ]: 0 : std::cout << " p2:" << p2 << "\n";
[ # # ]
1209 [ # # ][ # # ]: 0 : std::cout << " perpDir:" << perpDir << "\n";
[ # # ]
1210 [ # # ][ # # ]: 0 : std::cout << " boundary edges size:" << boundary_mesh_edges.size() << "\n";
[ # # ][ # # ]
1211 : : }
1212 [ + - ][ + - ]: 153 : for( ; ite != boundary_mesh_edges.end(); ++ite )
[ + - ][ + - ]
1213 : : {
1214 [ + - ]: 153 : EntityHandle candidateEdge = *ite;
1215 : : const EntityHandle* conn2;
1216 : : int nno;
1217 [ + - ]: 153 : rval = _mbImpl->get_connectivity( candidateEdge, conn2, nno );
1218 [ - + ][ # # ]: 153 : MBERRORR( rval, "Failed to get conn for boundary edge" );
[ # # ]
1219 [ + - ][ + + ]: 459 : CartVect pts[2];
1220 [ + - ][ + - ]: 153 : rval = _mbImpl->get_coords( conn2, 2, &( pts[0][0] ) );
1221 [ - + ][ # # ]: 153 : MBERRORR( rval, "Failed to get coords of nodes for boundary edge" );
[ # # ]
1222 [ + - ]: 153 : CartVect intx_point;
1223 : : double parPos;
1224 : : bool intersect =
1225 [ + - ]: 153 : intersect_segment_and_plane_slice( pts[0], pts[1], p1, p2, dirct, normPlane, intx_point, parPos );
1226 [ - + ]: 153 : if( debug_splits )
1227 : : {
1228 [ # # ][ # # ]: 0 : std::cout << " Edge:" << _mbImpl->id_from_handle( candidateEdge ) << "\n";
[ # # ][ # # ]
1229 [ # # ][ # # ]: 0 : std::cout << " Node 1:" << _mbImpl->id_from_handle( conn2[0] ) << pts[0] << "\n";
[ # # ][ # # ]
[ # # ]
1230 [ # # ][ # # ]: 0 : std::cout << " Node 2:" << _mbImpl->id_from_handle( conn2[1] ) << pts[1] << "\n";
[ # # ][ # # ]
[ # # ]
1231 [ # # ][ # # ]: 0 : std::cout << " Intersect bool:" << intersect << "\n";
[ # # ]
1232 : : }
1233 [ + + ]: 153 : if( intersect )
1234 : : {
1235 [ + - ][ + - ]: 3 : double proj1 = ( intx_point - p1 ) % perpDir;
1236 [ + - ][ + - ]: 3 : double proj2 = ( intx_point - p2 ) % perpDir;
1237 [ + + ]: 3 : if( ( fabs( proj1 ) > fabs( proj2 ) ) // this means it is closer to p2 than p1
1238 : : )
1239 : 1 : continue; // basically, this means the intersection point is with a
1240 : : // boundary edge on the other side, closer to p2 than p1, so we
1241 : : // skip it
1242 [ - + ]: 2 : if( parPos == 0 )
1243 : : {
1244 : : // close to vertex 1, nothing to do
1245 [ # # ]: 0 : nodes.push_back( conn2[0] );
1246 [ # # ]: 0 : splittingNodes.push_back( conn2[0] );
1247 : : }
1248 [ - + ]: 2 : else if( parPos == 1. )
1249 : : {
1250 : : // close to vertex 2, nothing to do
1251 [ # # ]: 0 : nodes.push_back( conn2[1] );
1252 [ # # ]: 0 : splittingNodes.push_back( conn2[1] );
1253 : : }
1254 : : else
1255 : : {
1256 : : // break the edge, create a new node at intersection point (will be smoothed
1257 : : // out)
1258 : : EntityHandle newVertex;
1259 [ + - ][ + - ]: 2 : rval = _mbImpl->create_vertex( &( intx_point[0] ), newVertex );
1260 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't create vertex" );
[ # # ]
1261 [ + - ]: 2 : nodes.push_back( newVertex );
1262 [ + - ]: 2 : split_internal_edge( candidateEdge, newVertex );
1263 [ + - ]: 2 : splittingNodes.push_back( newVertex );
1264 [ + - ]: 2 : _brokenEdges[newVertex] = candidateEdge;
1265 [ + - ]: 2 : _piercedEdges.insert( candidateEdge );
1266 : : }
1267 : 152 : break; // break from the loop over boundary edges, we are interested in the
1268 : : // first
1269 : : // split (hopefully, the only split)
1270 : : }
1271 : : }
1272 [ + - ][ + - ]: 2 : if( ite == boundary_mesh_edges.end() )
[ - + ]
1273 [ # # ][ # # ]: 2 : MBERRORR( MB_FAILURE, "Failed to find boundary intersection edge. Bail out" );
1274 : : }
1275 : : else
1276 : : {
1277 [ + - ]: 15 : std::vector< double > distances_out;
1278 [ + - ][ + - ]: 30 : std::vector< EntityHandle > sets_out;
1279 [ + - ][ + - ]: 30 : std::vector< EntityHandle > facets_out;
1280 : : rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
1281 [ + - ][ + - ]: 15 : tolerance, point, dir );
1282 [ - + ][ # # ]: 15 : MBERRORR( rval, "Failed to get ray intersections." );
[ # # ]
1283 [ - + ]: 15 : if( distances_out.size() < 1 )
1284 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, "Failed to get one intersection point, bad direction." );
1285 : :
1286 [ + + ]: 15 : if( distances_out.size() > 1 )
1287 [ + - ]: 7 : { std::cout << " too many intersection points. Only the first one considered\n"; }
1288 [ + - ]: 15 : std::vector< EntityHandle >::iterator pFace = std::find( sets_out.begin(), sets_out.end(), face );
1289 : :
1290 [ + - ][ - + ]: 15 : if( pFace == sets_out.end() ) MBERRORR( MB_FAILURE, "Failed to intersect given face, bad direction." );
[ # # ][ # # ]
1291 [ + - ]: 15 : unsigned int index = pFace - sets_out.begin();
1292 : : // get the closest node of the triangle, and modify locally the triangle(s), so the
1293 : : // intersection point is at a new vertex, if needed
1294 [ + - ]: 15 : CartVect P( point );
1295 [ + - ]: 15 : CartVect Dir( dir );
1296 [ + - ][ + - ]: 15 : CartVect newPoint = P + distances_out[index] * Dir;
[ + - ]
1297 : : // get the triangle coordinates
1298 : :
1299 : : double area_coord[3];
1300 : 15 : EntityHandle boundary_handle = 0; // if 0, not on a boundary
1301 [ + - ][ + - ]: 15 : rval = area_coordinates( _mbImpl, facets_out[index], newPoint, area_coord, boundary_handle );
1302 [ - + ][ # # ]: 15 : MBERRORR( rval, "Failed to get area coordinates" );
[ # # ]
1303 : :
1304 [ - + ]: 15 : if( debug_splits )
1305 : : {
1306 [ # # ][ # # ]: 0 : std::cout << " int point:" << newPoint << " area coord " << area_coord[0] << " " << area_coord[1] << " "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1307 [ # # ][ # # ]: 0 : << area_coord[2] << "\n";
1308 [ # # ][ # # ]: 0 : std::cout << " triangle: " << _mbImpl->id_from_handle( facets_out[index] )
[ # # ][ # # ]
1309 [ # # ][ # # ]: 0 : << " boundary:" << boundary_handle << "\n";
[ # # ]
1310 : : }
1311 : : EntityType type;
1312 [ + + ][ + - ]: 15 : if( boundary_handle ) type = _mbImpl->type_from_handle( boundary_handle );
1313 [ + + ][ + + ]: 15 : if( boundary_handle && ( type == MBVERTEX ) )
1314 : : {
1315 : : // nothing to do, we are happy
1316 [ + - ]: 6 : nodes.push_back( boundary_handle );
1317 : : }
1318 : : else
1319 : : {
1320 : : // for an edge, we will split 2 triangles
1321 : : // for interior point, we will create 3 triangles out of one
1322 : : // create a new vertex
1323 : : EntityHandle newVertex;
1324 [ + - ][ + - ]: 9 : rval = _mbImpl->create_vertex( &( newPoint[0] ), newVertex );
1325 [ + + ]: 9 : if( boundary_handle ) // this is edge
1326 : : {
1327 [ + - ]: 2 : split_internal_edge( boundary_handle, newVertex );
1328 [ + - ]: 2 : _piercedEdges.insert( boundary_handle ); // to be removed at the end
1329 : : }
1330 : : else
1331 [ + - ][ + - ]: 7 : divide_triangle( facets_out[index], newVertex );
1332 : :
1333 [ + - ][ + - ]: 15 : nodes.push_back( newVertex );
1334 : 15 : }
1335 : : }
1336 : : }
1337 : : // now, we have to find more intersection points, either interior to triangles, or on edges, or
1338 : : // on vertices use the same tolerance as before starting from 2 points on 2 triangles, and
1339 : : // having the direction, get more intersection points between the plane formed by direction and
1340 : : // those 2 points, and edges from triangulation (the triangles involved will be part of the same
1341 : : // gentity , original face ( moab set)
1342 : : // int closed = 1;// closed = 0 if the polyline is not closed
1343 : :
1344 [ + - ]: 4 : CartVect Dir( direction );
1345 [ + - ]: 8 : std::vector< EntityHandle > chainedEdges;
1346 : :
1347 [ + + ]: 20 : for( i = 0; i < numIniPoints - 1 + closed; i++ )
1348 : : {
1349 : 16 : int nextIndex = ( i + 1 ) % numIniPoints;
1350 [ + - ]: 16 : std::vector< EntityHandle > trianglesAlong;
1351 [ + - ][ + - ]: 32 : std::vector< CartVect > points;
1352 : : // otherwise to edges or even nodes
1353 [ + - ][ + - ]: 32 : std::vector< EntityHandle > entities;
1354 : : // start with initial points, intersect along the direction, find the facets
1355 [ + - ][ + - ]: 16 : rval = compute_intersection_points( face, nodes[i], nodes[nextIndex], Dir, points, entities, trianglesAlong );
[ + - ]
1356 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't get intersection points along a line" );
[ # # ]
1357 [ + - ][ + - ]: 32 : std::vector< EntityHandle > nodesAlongPolyline;
1358 : : // refactor code; move some edge creation for each 2 intersection points
1359 [ + - ][ + - ]: 16 : nodesAlongPolyline.push_back( entities[0] ); // it is for sure a node
1360 : 16 : int num_points = (int)points.size(); // it should be num_triangles + 1
1361 [ + + ]: 1791 : for( int j = 0; j < num_points - 1; j++ )
1362 : : {
1363 [ + - ]: 1775 : EntityHandle tri = trianglesAlong[j]; // this is happening in trianglesAlong i
1364 [ + - ]: 1775 : EntityHandle e1 = entities[j];
1365 [ + - ]: 1775 : EntityHandle e2 = entities[j + 1];
1366 [ + - ]: 1775 : EntityType et1 = _mbImpl->type_from_handle( e1 );
1367 : : // EntityHandle vertex1 = nodesAlongPolyline[i];// irrespective of the entity type i,
1368 : : // we already have the vertex there
1369 [ + - ]: 1775 : EntityType et2 = _mbImpl->type_from_handle( e2 );
1370 [ + + ][ + - ]: 1775 : if( et2 == MBVERTEX ) { nodesAlongPolyline.push_back( e2 ); }
1371 : : else // if (et2==MBEDGE)
1372 : : {
1373 [ + - ]: 1744 : CartVect coord_vert = points[j + 1];
1374 : : EntityHandle newVertex;
1375 [ + - ]: 1744 : rval = _mbImpl->create_vertex( (double*)&coord_vert, newVertex );
1376 [ - + ][ # # ]: 1744 : MBERRORR( rval, "can't create vertex" );
[ # # ]
1377 [ + - ]: 1744 : nodesAlongPolyline.push_back( newVertex );
1378 : : }
1379 : : // if vertices, do not split anything, just get the edge for polyline
1380 [ + + ][ - + ]: 1775 : if( et2 == MBVERTEX && et1 == MBVERTEX )
1381 : : {
1382 : : // nothing to do, just continue;
1383 : 0 : continue; // continue the for loop
1384 : : }
1385 : :
1386 [ - + ]: 1775 : if( debug_splits )
1387 : : {
1388 [ # # ][ # # ]: 0 : std::cout << "tri: type: " << _mbImpl->type_from_handle( tri )
[ # # ]
1389 [ # # ][ # # ]: 0 : << " id:" << _mbImpl->id_from_handle( tri ) << "\n e1:" << e1
[ # # ][ # # ]
[ # # ]
1390 [ # # ][ # # ]: 0 : << " id:" << _mbImpl->id_from_handle( e1 ) << " e2:" << e2
[ # # ][ # # ]
[ # # ]
1391 [ # # ][ # # ]: 0 : << " id:" << _mbImpl->id_from_handle( e2 ) << "\n";
[ # # ][ # # ]
1392 : : }
1393 : : // here, at least one is an edge
1394 [ + - ][ + - ]: 1775 : rval = BreakTriangle2( tri, e1, e2, nodesAlongPolyline[j], nodesAlongPolyline[j + 1] );
[ + - ]
1395 [ - + ][ # # ]: 1775 : MBERRORR( rval, "can't break triangle 2" );
[ # # ]
1396 [ + + ][ + - ]: 1775 : if( et2 == MBEDGE ) _piercedEdges.insert( e2 );
1397 [ + - ]: 1775 : _piercedTriangles.insert( tri );
1398 : : }
1399 : : // nodesAlongPolyline will define a new geometric edge
1400 [ - + ]: 16 : if( debug_splits )
1401 : : {
1402 [ # # ][ # # ]: 0 : std::cout << "nodesAlongPolyline: " << nodesAlongPolyline.size() << "\n";
[ # # ]
1403 [ # # ][ # # ]: 0 : std::cout << "points: " << num_points << "\n";
[ # # ]
1404 : : }
1405 : : // if needed, create edges along polyline, or revert the existing ones, to
1406 : : // put them in a new edge set
1407 : : EntityHandle new_geo_edge;
1408 [ + - ]: 16 : rval = create_new_gedge( nodesAlongPolyline, new_geo_edge );
1409 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't create a new edge" );
[ # # ]
1410 [ + - ][ + - ]: 16 : chainedEdges.push_back( new_geo_edge );
1411 : : // end copy
1412 : 16 : }
1413 : : // the segment between point_i and point_i+1 is in trianglesAlong_i
1414 : : // points_i is on entities_i
1415 : : // all these edges are oriented correctly
1416 [ + - ]: 4 : rval = split_surface( face, chainedEdges, splittingNodes, oNewFace );
1417 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't split surface" );
[ # # ]
1418 : : //
1419 [ + - ]: 4 : rval = chain_edges( min_dot ); // acos(0.8)~= 36 degrees
1420 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't chain edges" );
[ # # ]
1421 : 8 : return MB_SUCCESS;
1422 : : }
1423 : : /**
1424 : : * this method splits along the polyline defined by points and entities
1425 : : * the polyline will be defined with
1426 : : * // the entities are now only nodes and edges, no triangles!!!
1427 : : * the first and last ones are also nodes for sure
1428 : : */
1429 : 4 : ErrorCode FBEngine::split_surface( EntityHandle face, std::vector< EntityHandle >& chainedEdges,
1430 : : std::vector< EntityHandle >& splittingNodes, EntityHandle& newFace )
1431 : : {
1432 : : // use now the chained edges to create a new face (loop or clean split)
1433 : : // use a fill to determine the new sets, up to the polyline
1434 : : // points on the polyline will be moved to the closest point location, with some constraints
1435 : : // then the sets will be reset, geometry recomputed. new vertices, new edges, etc.
1436 : :
1437 [ + - ]: 4 : Range iniTris;
1438 : : ErrorCode rval;
1439 [ + - ]: 4 : rval = _mbImpl->get_entities_by_type( face, MBTRI, iniTris );
1440 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get initial triangles" );
[ # # ]
1441 : :
1442 : : // start from a triangle that is not in the triangles to delete
1443 : : // flood fill
1444 : :
1445 : 4 : bool closed = splittingNodes.size() == 0;
1446 [ + + ]: 4 : if( !closed )
1447 : : {
1448 : : //
1449 [ - + ][ # # ]: 1 : if( splittingNodes.size() != 2 ) MBERRORR( MB_FAILURE, "need to have exactly 2 nodes for splitting" );
[ # # ]
1450 : : // we will have to split the boundary edges
1451 : : // first, find the actual boundary, and try to split with the 2 end points (nodes)
1452 : : // get the adjacent edges, and see which one has the end nodes
1453 [ + - ][ + - ]: 1 : rval = split_boundary( face, splittingNodes[0] );
1454 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't split with first node" );
[ # # ]
1455 [ + - ][ + - ]: 1 : rval = split_boundary( face, splittingNodes[1] );
1456 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't split with second node)" );
[ # # ]
1457 : : }
1458 : : // we will separate triangles to delete, unaffected, new_triangles,
1459 : : // nodesAlongPolyline,
1460 [ + - ][ + - ]: 8 : Range first, second;
1461 [ + - ]: 4 : rval = separate( face, chainedEdges, first, second );
1462 : :
1463 : : // now, we are done with the computations;
1464 : : // we need to put the new nodes on the smooth surface
1465 [ + - ]: 4 : if( this->_smooth )
1466 : : {
1467 [ + - ]: 4 : rval = smooth_new_intx_points( face, chainedEdges );
1468 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't smooth new points" );
[ # # ]
1469 : : }
1470 : :
1471 : : // create the new set
1472 [ + - ]: 4 : rval = _mbImpl->create_meshset( MESHSET_SET, newFace );
1473 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't create a new face" );
[ # # ]
1474 : :
1475 [ + - ]: 4 : _my_geomTopoTool->add_geo_set( newFace, 2 );
1476 : :
1477 : : // the new face will have the first set (positive sense triangles, to the left)
1478 [ + - ]: 4 : rval = _mbImpl->add_entities( newFace, first );
1479 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add first range triangles to new face" );
[ # # ]
1480 : :
1481 [ + + ]: 20 : for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1482 : : {
1483 [ + - ]: 16 : EntityHandle new_geo_edge = chainedEdges[j];
1484 : : // both faces will have the edge now
1485 [ + - ]: 16 : rval = _mbImpl->add_parent_child( face, new_geo_edge );
1486 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't add parent child relations for new edge" );
[ # # ]
1487 : :
1488 [ + - ]: 16 : rval = _mbImpl->add_parent_child( newFace, new_geo_edge );
1489 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't add parent child relations for new edge" );
[ # # ]
1490 : : // add senses
1491 : : // sense newFace is 1, old face is -1
1492 [ + - ]: 16 : rval = _my_geomTopoTool->set_sense( new_geo_edge, newFace, 1 );
1493 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't set sense for new edge" );
[ # # ]
1494 : :
1495 [ + - ]: 16 : rval = _my_geomTopoTool->set_sense( new_geo_edge, face, -1 );
1496 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't set sense for new edge in original face" );
[ # # ]
1497 : : }
1498 : :
1499 [ + - ]: 4 : rval = set_neumann_tags( face, newFace );
1500 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set NEUMANN set tags" );
[ # # ]
1501 : :
1502 : : // now, we should remove from the original set all tris, and put the "second" range
1503 [ + - ]: 4 : rval = _mbImpl->remove_entities( face, iniTris );
1504 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't remove original tris from initial face set" );
[ # # ]
1505 : :
1506 [ + - ]: 4 : rval = _mbImpl->add_entities( face, second );
1507 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add second range to the original set" );
[ # # ]
1508 : :
1509 [ + + ]: 4 : if( !closed )
1510 : : {
1511 [ + - ]: 1 : rval = redistribute_boundary_edges_to_faces( face, newFace, chainedEdges );
1512 [ - + ][ # # ]: 1 : MBERRORR( rval, "fail to reset the proper boundary faces" );
[ # # ]
1513 : : }
1514 : :
1515 : : /*if (_smooth)
1516 : : delete_smooth_tags();// they need to be recomputed, anyway
1517 : : // this will remove the extra smooth faces and edges
1518 : : clean();*/
1519 : : // also, these nodes need to be moved to the smooth surface, sometimes before deleting the old
1520 : : // triangles
1521 : : // remove the triangles from the set, then delete triangles (also some edges need to be
1522 : : // deleted!)
1523 [ + - ]: 4 : rval = _mbImpl->delete_entities( _piercedTriangles );
1524 : :
1525 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't delete triangles" );
[ # # ]
1526 [ + - ]: 4 : _piercedTriangles.clear();
1527 : : // delete edges that are broke up in 2
1528 [ + - ]: 4 : rval = _mbImpl->delete_entities( _piercedEdges );
1529 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't delete edges" );
[ # # ]
1530 [ + - ]: 4 : _piercedEdges.clear();
1531 : :
1532 [ - + ]: 4 : if( debug_splits )
1533 : : {
1534 [ # # ]: 0 : _mbImpl->write_file( "newFace.vtk", "vtk", 0, &newFace, 1 );
1535 [ # # ]: 0 : _mbImpl->write_file( "leftoverFace.vtk", "vtk", 0, &face, 1 );
1536 : : }
1537 : 8 : return MB_SUCCESS;
1538 : : }
1539 : :
1540 : 4 : ErrorCode FBEngine::smooth_new_intx_points( EntityHandle face, std::vector< EntityHandle >& chainedEdges )
1541 : : {
1542 : :
1543 : : // do not move nodes from the original face
1544 : : // first get all triangles, and then all nodes from those triangles
1545 : :
1546 [ + - ]: 4 : Range tris;
1547 [ + - ]: 4 : ErrorCode rval = _mbImpl->get_entities_by_type( face, MBTRI, tris );
1548 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get triangles" );
[ # # ]
1549 : :
1550 [ + - ]: 8 : Range ini_nodes;
1551 [ + - ]: 4 : rval = _mbImpl->get_connectivity( tris, ini_nodes );
1552 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get connectivities" );
[ # # ]
1553 : :
1554 [ + - ]: 4 : SmoothFace* smthFace = _faces[face];
1555 : :
1556 : : // get all nodes from chained edges
1557 [ + - ]: 8 : Range mesh_edges;
1558 [ + + ]: 20 : for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1559 : : {
1560 : : // keep adding to the range of mesh edges
1561 [ + - ][ + - ]: 16 : rval = _mbImpl->get_entities_by_dimension( chainedEdges[j], 1, mesh_edges );
1562 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't get mesh edges" );
[ # # ]
1563 : : }
1564 : : // nodes along polyline
1565 [ + - ]: 8 : Range nodes_on_polyline;
1566 [ + - ]: 4 : rval = _mbImpl->get_connectivity( mesh_edges, nodes_on_polyline, true ); // corners only
1567 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get nodes on the polyline" );
[ # # ]
1568 : :
1569 [ + - ]: 8 : Range new_intx_nodes = subtract( nodes_on_polyline, ini_nodes );
1570 : :
1571 [ + - ]: 8 : std::vector< double > ini_coords;
1572 [ + - ]: 4 : int num_points = (int)new_intx_nodes.size();
1573 [ + - ]: 4 : ini_coords.resize( 3 * num_points );
1574 [ + - ][ + - ]: 4 : rval = _mbImpl->get_coords( new_intx_nodes, &( ini_coords[0] ) );
1575 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get coordinates" );
[ # # ]
1576 : :
1577 : 4 : int i = 0;
1578 [ + - ][ + - ]: 1759 : for( Range::iterator it = new_intx_nodes.begin(); it != new_intx_nodes.end(); ++it )
[ + - ][ + - ]
[ + + ]
1579 : : {
1580 : : /*EntityHandle node = *it;*/
1581 : 1755 : int i3 = 3 * i;
1582 [ + - ][ + - ]: 1755 : smthFace->move_to_surface( ini_coords[i3], ini_coords[i3 + 1], ini_coords[i3 + 2] );
[ + - ][ + - ]
1583 : : // reset the coordinates of this node
1584 : 1755 : ++i;
1585 : : }
1586 [ + - ][ + - ]: 4 : rval = _mbImpl->set_coords( new_intx_nodes, &( ini_coords[0] ) );
[ + - ]
1587 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set smoothed coordinates for the new nodes" );
[ # # ]
1588 : :
1589 : 8 : return MB_SUCCESS;
1590 : : }
1591 : : // we will use the fact that the splitting edge is oriented right now
1592 : : // to the left will be new face, to the right, old face
1593 : : // (to the left, positively oriented triangles)
1594 : 4 : ErrorCode FBEngine::separate( EntityHandle face, std::vector< EntityHandle >& chainedEdges, Range& first,
1595 : : Range& second )
1596 : : {
1597 : : // Range unaffectedTriangles = subtract(iniTriangles, _piercedTriangles);
1598 : : // insert in each
1599 : : // start with a new triangle, and flood to get the first range; what is left is the
1600 : : // second range
1601 : : // flood fill is considering edges adjacent to seed triangles; if there is
1602 : : // an edge in the new_geo_edge, it is skipped; triangles in the
1603 : : // triangles to delete are not added
1604 : : // first, create all edges of the new triangles
1605 : :
1606 : : //
1607 : : // new face will have the new edge oriented positively
1608 : : // get mesh edges from geo edge (splitting gedge);
1609 : :
1610 [ + - ]: 4 : Range mesh_edges;
1611 : : ErrorCode rval;
1612 : : // mesh_edges
1613 [ + + ]: 20 : for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1614 : : {
1615 : : // this will keep adding edges to the mesh_edges range
1616 : : // when we get out, the mesh_edges will be in this range, but not ordered
1617 [ + - ][ + - ]: 16 : rval = _mbImpl->get_entities_by_type( chainedEdges[j], MBEDGE, mesh_edges );
1618 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't get new polyline edges" );
[ # # ]
1619 [ - + ]: 16 : if( debug_splits )
1620 : : {
1621 [ # # ][ # # ]: 0 : std::cout << " At chained edge " << j << " " << _mbImpl->id_from_handle( chainedEdges[j] )
[ # # ][ # # ]
[ # # ][ # # ]
1622 [ # # ][ # # ]: 0 : << " mesh_edges Range size:" << mesh_edges.size() << "\n";
[ # # ][ # # ]
1623 : : }
1624 : : }
1625 : :
1626 : : // get a positive triangle adjacent to mesh_edge[0]
1627 : : // add to first triangles to the left, second triangles to the right of the mesh_edges ;
1628 : :
1629 : : // create a temp tag, and when done, delete it
1630 : : // default value: 0
1631 : : // 3 to be deleted, pierced
1632 : : // 1 first set
1633 : : // 2 second set
1634 : : // create the tag also for control points on the edges
1635 : 4 : int defVal = 0;
1636 : : Tag separateTag;
1637 [ + - ]: 4 : rval = MBI->tag_get_handle( "SEPARATE_TAG", 1, MB_TYPE_INTEGER, separateTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
1638 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't create temp tag for separation" );
[ # # ]
1639 : : // the deleted triangles will get a value 3, from start
1640 : 4 : int delVal = 3;
1641 [ + - ][ + - ]: 1792 : for( Range::iterator it1 = this->_piercedTriangles.begin(); it1 != _piercedTriangles.end(); ++it1 )
[ + - ][ + - ]
[ + + ]
1642 : : {
1643 [ + - ]: 1788 : EntityHandle trToDelete = *it1;
1644 [ + - ]: 1788 : rval = _mbImpl->tag_set_data( separateTag, &trToDelete, 1, &delVal );
1645 [ - + ][ # # ]: 1788 : MBERRORR( rval, "can't set delete tag value" );
[ # # ]
1646 : : }
1647 : :
1648 : : // find a triangle that will be in the first range, positively oriented about the splitting edge
1649 : 4 : EntityHandle seed1 = 0;
1650 [ + - ][ + - ]: 8 : for( Range::iterator it = mesh_edges.begin(); it != mesh_edges.end() && !seed1; ++it )
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ]
[ + + # # ]
1651 : : {
1652 [ + - ]: 4 : EntityHandle meshEdge = *it;
1653 [ + - ]: 4 : Range adj_tri;
1654 [ + - ]: 4 : rval = _mbImpl->get_adjacencies( &meshEdge, 1, 2, false, adj_tri );
1655 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get adj_tris to mesh edge" );
[ # # ]
1656 : :
1657 [ + - ][ + - ]: 12 : for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
[ + - ][ + - ]
[ + - ][ + - ]
1658 : : {
1659 [ + - ]: 8 : EntityHandle tr = *it2;
1660 [ + - ][ + - ]: 8 : if( _piercedTriangles.find( tr ) != _piercedTriangles.end() )
[ + - ][ - + ]
1661 : 0 : continue; // do not attach pierced triangles, they are not good
1662 : : int num1, sense, offset;
1663 [ + - ]: 8 : rval = _mbImpl->side_number( tr, meshEdge, num1, sense, offset );
1664 [ - + ][ # # ]: 8 : MBERRORR( rval, "edge not adjacent" );
[ # # ]
1665 [ + + ]: 8 : if( sense == 1 )
1666 : : {
1667 : : // firstSet.insert(tr);
1668 [ + - ]: 4 : if( !seed1 )
1669 : : {
1670 : 4 : seed1 = tr;
1671 : 8 : break;
1672 : : }
1673 : : }
1674 : : }
1675 : 4 : }
1676 : :
1677 : : // flood fill first set, the rest will be in second set
1678 : : // the edges from new_geo_edge will not be crossed
1679 : :
1680 : : // get edges of face (adjacencies)
1681 : : // also get the old boundary edges, from face; they will be edges to not cross, too
1682 [ + - ]: 8 : Range bound_edges;
1683 [ + - ]: 4 : rval = getAdjacentEntities( face, 1, bound_edges );
1684 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get boundary edges" );
[ # # ]
1685 : :
1686 : : // add to the do not cross edges range, all edges from initial boundary
1687 [ + - ]: 8 : Range initialBoundaryEdges;
1688 [ + - ][ + - ]: 10 : for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
[ + - ][ + - ]
[ + + ]
1689 : : {
1690 [ + - ]: 6 : EntityHandle bound_edge = *it;
1691 [ + - ]: 6 : rval = _mbImpl->get_entities_by_dimension( bound_edge, 1, initialBoundaryEdges );
1692 : : }
1693 : :
1694 [ + - ]: 8 : Range doNotCrossEdges = unite( initialBoundaryEdges, mesh_edges ); // add the splitting edges !
1695 : :
1696 : : // use a second method, with tags
1697 : : //
1698 [ + - ][ + - ]: 8 : std::queue< EntityHandle > queue1;
1699 [ + - ]: 4 : queue1.push( seed1 );
1700 [ + - ]: 8 : std::vector< EntityHandle > arr1;
1701 [ + - ][ + + ]: 17575 : while( !queue1.empty() )
1702 : : {
1703 : : // start copy
1704 [ + - ]: 17571 : EntityHandle currentTriangle = queue1.front();
1705 [ + - ]: 17571 : queue1.pop();
1706 [ + - ]: 17571 : arr1.push_back( currentTriangle );
1707 : : // add new triangles that share an edge
1708 [ + - ]: 17571 : Range currentEdges;
1709 [ + - ]: 17571 : rval = _mbImpl->get_adjacencies( ¤tTriangle, 1, 1, true, currentEdges, Interface::UNION );
1710 [ - + ][ # # ]: 17571 : MBERRORR( rval, "can't get adjacencies" );
[ # # ]
1711 [ + - ][ + - ]: 70284 : for( Range::iterator it = currentEdges.begin(); it != currentEdges.end(); ++it )
[ + - ][ + - ]
[ + + ][ + - ]
1712 : : {
1713 [ + - ]: 52713 : EntityHandle frontEdge = *it;
1714 [ + - ][ + - ]: 52713 : if( doNotCrossEdges.find( frontEdge ) == doNotCrossEdges.end() )
[ + - ][ + + ]
1715 : : {
1716 : : // this is an edge that can be crossed
1717 [ + - ]: 50828 : Range adj_tri;
1718 [ + - ]: 50828 : rval = _mbImpl->get_adjacencies( &frontEdge, 1, 2, false, adj_tri, Interface::UNION );
1719 [ - + ][ # # ]: 50828 : MBERRORR( rval, "can't get adj_tris" );
[ # # ]
1720 : : // if the triangle is not in first range, add it to the queue
1721 [ + - ][ + - ]: 154336 : for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
[ + - ][ + - ]
[ + + ][ + - ]
1722 : : {
1723 [ + - ]: 103508 : EntityHandle tri2 = *it2;
1724 : 103508 : int val = 0;
1725 [ + - ]: 103508 : rval = _mbImpl->tag_get_data( separateTag, &tri2, 1, &val );
1726 [ - + ][ # # ]: 103508 : MBERRORR( rval, "can't get tag value" );
[ # # ]
1727 [ + + ]: 103508 : if( val ) continue;
1728 : : // else, set it to 1
1729 : 17567 : val = 1;
1730 [ + - ]: 17567 : rval = _mbImpl->tag_set_data( separateTag, &tri2, 1, &val );
1731 [ - + ][ # # ]: 17567 : MBERRORR( rval, "can't get tag value" );
[ # # ]
1732 : :
1733 [ + - ]: 17567 : queue1.push( tri2 );
1734 : 50828 : }
1735 : : } // end edge do not cross
1736 : : } // end while
1737 : 17571 : }
1738 : :
1739 [ + - ]: 4 : std::sort( arr1.begin(), arr1.end() );
1740 : : // Range first1;
1741 [ + - ][ + - ]: 4 : std::copy( arr1.rbegin(), arr1.rend(), range_inserter( first ) );
1742 : :
1743 : : // std::cout<< "\n first1.size() " << first1.size() << " first.size(): " << first.size() <<
1744 : : // "\n";
1745 [ - + ]: 4 : if( debug_splits )
1746 : : {
1747 : : EntityHandle tmpSet;
1748 [ # # ]: 0 : _mbImpl->create_meshset( MESHSET_SET, tmpSet );
1749 [ # # ]: 0 : _mbImpl->add_entities( tmpSet, first );
1750 [ # # ]: 0 : _mbImpl->write_file( "dbg1.vtk", "vtk", 0, &tmpSet, 1 );
1751 : : }
1752 : : // now, decide the set 2:
1753 : : // first, get all ini tris
1754 [ + - ]: 8 : Range initr;
1755 [ + - ]: 4 : rval = _mbImpl->get_entities_by_type( face, MBTRI, initr );
1756 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get tris " );
[ # # ]
1757 [ + - ][ + - ]: 4 : second = unite( initr, _newTriangles );
1758 [ + - ]: 8 : Range second2 = subtract( second, _piercedTriangles );
1759 [ + - ][ + - ]: 4 : second = subtract( second2, first );
1760 [ + - ]: 4 : _newTriangles.clear();
1761 [ - + ]: 4 : if( debug_splits )
1762 : : {
1763 [ # # ][ # # ]: 0 : std::cout << "\n second.size() " << second.size() << " first.size(): " << first.size() << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1764 : : // debugging code
1765 : : EntityHandle tmpSet2;
1766 [ # # ]: 0 : _mbImpl->create_meshset( MESHSET_SET, tmpSet2 );
1767 [ # # ]: 0 : _mbImpl->add_entities( tmpSet2, second );
1768 [ # # ]: 0 : _mbImpl->write_file( "dbg2.vtk", "vtk", 0, &tmpSet2, 1 );
1769 : : }
1770 : : /*Range intex = intersect(first, second);
1771 : : if (!intex.empty() && debug_splits)
1772 : : {
1773 : : std::cout << "error, the sets should be disjoint\n";
1774 : : for (Range::iterator it1=intex.begin(); it1!=intex.end(); ++it1)
1775 : : {
1776 : : std::cout<<_mbImpl->id_from_handle(*it1) << "\n";
1777 : : }
1778 : : }*/
1779 [ + - ]: 4 : rval = _mbImpl->tag_delete( separateTag );
1780 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't delete tag " );
[ # # ]
1781 : 8 : return MB_SUCCESS;
1782 : : }
1783 : : // if there is an edge between 2 nodes, then check it's orientation, and revert it if needed
1784 : 16 : ErrorCode FBEngine::create_new_gedge( std::vector< EntityHandle >& nodesAlongPolyline, EntityHandle& new_geo_edge )
1785 : : {
1786 : :
1787 [ + - ]: 16 : ErrorCode rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_geo_edge );
1788 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't create geo edge" );
[ # # ]
1789 : :
1790 : : // now, get the edges, or create if not existing
1791 [ + - ]: 16 : std::vector< EntityHandle > mesh_edges;
1792 [ + + ]: 1791 : for( unsigned int i = 0; i < nodesAlongPolyline.size() - 1; i++ )
1793 : : {
1794 [ + - ][ + - ]: 1775 : EntityHandle n1 = nodesAlongPolyline[i], n2 = nodesAlongPolyline[i + 1];
1795 : :
1796 : : EntityHandle nn2[2];
1797 : 1775 : nn2[0] = n1;
1798 : 1775 : nn2[1] = n2;
1799 : :
1800 [ + - ]: 1775 : std::vector< EntityHandle > adjacent;
1801 [ + - ]: 1775 : rval = _mbImpl->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT );
1802 : : // see if there is an edge between those 2 already, and if it is oriented as we like
1803 : 1775 : bool new_edge = true;
1804 [ - + ]: 1775 : if( adjacent.size() >= 1 )
1805 : : {
1806 : : // check the orientation
1807 : 0 : const EntityHandle* conn2 = NULL;
1808 : 0 : int nnodes = 0;
1809 [ # # ][ # # ]: 0 : rval = _mbImpl->get_connectivity( adjacent[0], conn2, nnodes );
1810 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't get connectivity" );
[ # # ]
1811 [ # # ][ # # ]: 0 : if( conn2[0] == nn2[0] && conn2[1] == nn2[1] )
1812 : : {
1813 : : // everything is fine
1814 [ # # ][ # # ]: 0 : mesh_edges.push_back( adjacent[0] );
1815 : 0 : new_edge = false; // we found one that's good, no need to create a new one
1816 : : }
1817 : : else
1818 : : {
1819 [ # # ][ # # ]: 0 : _piercedEdges.insert( adjacent[0] ); // we want to remove this one, it will be not needed
1820 : : }
1821 : : }
1822 [ + - ]: 1775 : if( new_edge )
1823 : : {
1824 : : // there is no edge between n1 and n2, create one
1825 : : EntityHandle mesh_edge;
1826 [ + - ]: 1775 : rval = _mbImpl->create_element( MBEDGE, nn2, 2, mesh_edge );
1827 [ - + ][ # # ]: 1775 : MBERRORR( rval, "Failed to create a new edge" );
[ # # ]
1828 [ + - ][ + - ]: 1775 : mesh_edges.push_back( mesh_edge );
1829 : : }
1830 : 1775 : }
1831 : :
1832 : : // add loops edges to the edge set
1833 [ + - ]: 16 : rval = _mbImpl->add_entities( new_geo_edge, &mesh_edges[0],
1834 [ + - ]: 32 : mesh_edges.size() ); // only one edge
1835 [ - + ][ # # ]: 16 : MBERRORR( rval, "can't add edges to new_geo_set" );
[ # # ]
1836 : : // check vertex sets for vertex 1 and vertex 2?
1837 : : // get all sets of dimension 0 from database, and see if our ends are here or not
1838 : :
1839 [ + - ]: 32 : Range ends_geo_edge;
1840 [ + - ][ + - ]: 16 : ends_geo_edge.insert( nodesAlongPolyline[0] );
1841 [ + - ][ + - ]: 16 : ends_geo_edge.insert( nodesAlongPolyline[nodesAlongPolyline.size() - 1] );
1842 : :
1843 [ + - ][ + + ]: 48 : for( unsigned int k = 0; k < ends_geo_edge.size(); k++ )
1844 : : {
1845 [ + - ]: 32 : EntityHandle node = ends_geo_edge[k];
1846 : : EntityHandle nodeSet;
1847 [ + - ]: 32 : bool found = find_vertex_set_for_node( node, nodeSet );
1848 : :
1849 [ + + ]: 32 : if( !found )
1850 : : {
1851 : : // create a node set and add the node
1852 : :
1853 [ + - ]: 17 : rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
1854 [ - + ][ # # ]: 17 : MBERRORR( rval, "Failed to create a new vertex set" );
[ # # ]
1855 : :
1856 [ + - ]: 17 : rval = _mbImpl->add_entities( nodeSet, &node, 1 );
1857 [ - + ][ # # ]: 17 : MBERRORR( rval, "Failed to add the node to the set" );
[ # # ]
1858 : :
1859 [ + - ]: 17 : rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
1860 [ - + ][ # # ]: 17 : MBERRORR( rval, "Failed to commit the node set" );
[ # # ]
1861 : :
1862 [ - + ]: 17 : if( debug_splits )
1863 : : {
1864 [ # # ][ # # ]: 0 : std::cout << " create a vertex set " << _mbImpl->id_from_handle( nodeSet )
[ # # ]
1865 [ # # ][ # # ]: 0 : << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << " for node " << node
[ # # ][ # # ]
[ # # ]
1866 [ # # ]: 0 : << "\n";
1867 : : }
1868 : : }
1869 : :
1870 [ + - ]: 32 : rval = _mbImpl->add_parent_child( new_geo_edge, nodeSet );
1871 [ - + ][ # # ]: 32 : MBERRORR( rval, "Failed to add parent child relation" );
[ # # ]
1872 : : }
1873 : : // finally, put the edge in the range of edges
1874 [ + - ][ - + ]: 16 : rval = _my_geomTopoTool->add_geo_set( new_geo_edge, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
1875 : :
1876 : 32 : return rval;
1877 : : }
1878 : :
1879 : 0 : void FBEngine::print_debug_triangle( EntityHandle t )
1880 : : {
1881 [ # # ][ # # ]: 0 : std::cout << " triangle id:" << _mbImpl->id_from_handle( t ) << "\n";
[ # # ][ # # ]
1882 : 0 : const EntityHandle* conn3 = NULL;
1883 : 0 : int nnodes = 0;
1884 [ # # ]: 0 : _mbImpl->get_connectivity( t, conn3, nnodes );
1885 : : // get coords
1886 [ # # ][ # # ]: 0 : CartVect P[3];
1887 [ # # ]: 0 : _mbImpl->get_coords( conn3, 3, (double*)&P[0] );
1888 [ # # ][ # # ]: 0 : std::cout << " nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1889 [ # # ][ # # ]: 0 : CartVect PP[3];
1890 [ # # ]: 0 : PP[0] = P[1] - P[0];
1891 [ # # ]: 0 : PP[1] = P[2] - P[1];
1892 [ # # ]: 0 : PP[2] = P[0] - P[2];
1893 : :
1894 [ # # ][ # # ]: 0 : std::cout << " pos:" << P[0] << " " << P[1] << " " << P[2] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1895 [ # # ][ # # ]: 0 : std::cout << " x,y diffs " << PP[0][0] << " " << PP[0][1] << ", " << PP[1][0] << " " << PP[1][1] << ", "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1896 [ # # ][ # # ]: 0 : << PP[2][0] << " " << PP[2][1] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
1897 : 0 : return;
1898 : : }
1899 : : // actual breaking of triangles
1900 : : // case 1: n2 interior to triangle
1901 : 0 : ErrorCode FBEngine::BreakTriangle( EntityHandle, EntityHandle, EntityHandle, EntityHandle, EntityHandle, EntityHandle )
1902 : : {
1903 : 0 : std::cout << "FBEngine::BreakTriangle not implemented yet\n";
1904 : 0 : return MB_FAILURE;
1905 : : }
1906 : : // case 2, n1 and n2 on boundary
1907 : 1775 : ErrorCode FBEngine::BreakTriangle2( EntityHandle tri, EntityHandle e1, EntityHandle e2, EntityHandle n1,
1908 : : EntityHandle n2 ) // nodesAlongPolyline are on entities!
1909 : : {
1910 : : // we have the nodes, we just need to reconnect to form new triangles
1911 : : ErrorCode rval;
1912 : 1775 : const EntityHandle* conn3 = NULL;
1913 : 1775 : int nnodes = 0;
1914 [ + - ]: 1775 : rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
1915 [ - + ][ # # ]: 1775 : MBERRORR( rval, "Failed to get connectivity" );
[ # # ]
1916 [ - + ]: 1775 : assert( 3 == nnodes );
1917 : :
1918 [ + - ]: 1775 : EntityType et1 = _mbImpl->type_from_handle( e1 );
1919 [ + - ]: 1775 : EntityType et2 = _mbImpl->type_from_handle( e2 );
1920 : :
1921 [ + + ]: 1775 : if( MBVERTEX == et1 )
1922 : : {
1923 : : // find the vertex in conn3, and form 2 other triangles
1924 : 31 : int index = -1;
1925 [ + - ]: 64 : for( index = 0; index < 3; index++ )
1926 : : {
1927 [ + + ]: 64 : if( conn3[index] == e1 ) // also n1
1928 : 31 : break;
1929 : : }
1930 [ - + ]: 31 : if( index == 3 ) return MB_FAILURE;
1931 : : // 2 triangles: n1, index+1, n2, and n1, n2, index+2
1932 : 31 : EntityHandle conn[6] = { n1, conn3[( index + 1 ) % 3], n2, n1, n2, conn3[( index + 2 ) % 3] };
1933 : : EntityHandle newTriangle;
1934 [ + - ]: 31 : rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
1935 [ - + ][ # # ]: 31 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
1936 [ + - ]: 31 : _newTriangles.insert( newTriangle );
1937 [ - + ][ # # ]: 31 : if( debug_splits ) print_debug_triangle( newTriangle );
1938 [ + - ]: 31 : rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
1939 [ - + ][ # # ]: 31 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
1940 [ + - ]: 31 : _newTriangles.insert( newTriangle );
1941 [ - + ][ # # ]: 31 : if( debug_splits ) print_debug_triangle( newTriangle );
1942 : : }
1943 [ + + ]: 1744 : else if( MBVERTEX == et2 )
1944 : : {
1945 : 31 : int index = -1;
1946 [ + - ]: 72 : for( index = 0; index < 3; index++ )
1947 : : {
1948 [ + + ]: 72 : if( conn3[index] == e2 ) // also n2
1949 : 31 : break;
1950 : : }
1951 [ - + ]: 31 : if( index == 3 ) return MB_FAILURE;
1952 : : // 2 triangles: n1, index+1, n2, and n1, n2, index+2
1953 : 31 : EntityHandle conn[6] = { n2, conn3[( index + 1 ) % 3], n1, n2, n1, conn3[( index + 2 ) % 3] };
1954 : : EntityHandle newTriangle;
1955 [ + - ]: 31 : rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
1956 [ - + ][ # # ]: 31 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
1957 [ + - ]: 31 : _newTriangles.insert( newTriangle );
1958 [ - + ][ # # ]: 31 : if( debug_splits ) print_debug_triangle( newTriangle );
1959 [ + - ]: 31 : rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
1960 [ - + ][ # # ]: 31 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
1961 [ + - ]: 31 : _newTriangles.insert( newTriangle );
1962 [ - + ][ # # ]: 31 : if( debug_splits ) print_debug_triangle( newTriangle );
1963 : : }
1964 : : else
1965 : : {
1966 : : // both are edges adjacent to triangle tri
1967 : : // there are several configurations possible for n1, n2, between conn3 nodes.
1968 : : int num1, num2, sense, offset;
1969 [ + - ]: 1713 : rval = _mbImpl->side_number( tri, e1, num1, sense, offset );
1970 [ - + ][ # # ]: 1713 : MBERRORR( rval, "edge not adjacent" );
[ # # ]
1971 : :
1972 [ + - ]: 1713 : rval = _mbImpl->side_number( tri, e2, num2, sense, offset );
1973 [ - + ][ # # ]: 1713 : MBERRORR( rval, "edge not adjacent" );
[ # # ]
1974 : :
1975 : : const EntityHandle* conn12; // connectivity for edge 1
1976 : : const EntityHandle* conn22; // connectivity for edge 2
1977 : : // int nnodes;
1978 [ + - ]: 1713 : rval = _mbImpl->get_connectivity( e1, conn12, nnodes );
1979 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to get connectivity of edge 1" );
[ # # ]
1980 [ - + ]: 1713 : assert( 2 == nnodes );
1981 [ + - ]: 1713 : rval = _mbImpl->get_connectivity( e2, conn22, nnodes );
1982 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to get connectivity of edge 2" );
[ # # ]
1983 [ - + ]: 1713 : assert( 2 == nnodes );
1984 : : // now, having the side number, work through
1985 [ - + ]: 1713 : if( debug_splits )
1986 : : {
1987 [ # # ][ # # ]: 0 : std::cout << "tri conn3:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1988 [ # # ][ # # ]: 0 : std::cout << " edge1: conn12:" << conn12[0] << " " << conn12[1] << " side: " << num1 << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1989 [ # # ][ # # ]: 0 : std::cout << " edge2: conn22:" << conn22[0] << " " << conn22[1] << " side: " << num2 << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1990 : : }
1991 : 1713 : int unaffectedSide = 3 - num1 - num2;
1992 : 1713 : int i3 = ( unaffectedSide + 2 ) % 3; // to 0 is 2, to 1 is 0, to 2 is 1
1993 : : // triangles will be formed with triVertexIndex , n1, n2 (in what order?)
1994 : : EntityHandle v1, v2; // to hold the 2 nodes on edges
1995 [ + + ]: 1713 : if( num1 == i3 )
1996 : : {
1997 : 848 : v1 = n1;
1998 : 848 : v2 = n2;
1999 : : }
2000 : : else // if (num2==i3)
2001 : : {
2002 : 865 : v1 = n2;
2003 : 865 : v2 = n1;
2004 : : }
2005 : : // three triangles are formed
2006 : 1713 : int i1 = ( i3 + 1 ) % 3;
2007 : 1713 : int i2 = ( i3 + 2 ) % 3;
2008 : : // we could break the surface differently
2009 : 1713 : EntityHandle conn[9] = { conn3[i3], v1, v2, v1, conn3[i1], conn3[i2], v2, v1, conn3[i2] };
2010 : : EntityHandle newTriangle;
2011 [ - + ][ # # ]: 1713 : if( debug_splits ) std::cout << "Split 2 edges :\n";
2012 [ + - ]: 1713 : rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
2013 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
2014 [ + - ]: 1713 : _newTriangles.insert( newTriangle );
2015 [ - + ][ # # ]: 1713 : if( debug_splits ) print_debug_triangle( newTriangle );
2016 [ + - ]: 1713 : rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
2017 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
2018 [ + - ]: 1713 : _newTriangles.insert( newTriangle );
2019 [ - + ][ # # ]: 1713 : if( debug_splits ) print_debug_triangle( newTriangle );
2020 [ + - ]: 1713 : rval = _mbImpl->create_element( MBTRI, conn + 6, 3, newTriangle ); // the second triangle
2021 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to create a new triangle" );
[ # # ]
2022 [ + - ]: 1713 : _newTriangles.insert( newTriangle );
2023 [ - + ][ # # ]: 1713 : if( debug_splits ) print_debug_triangle( newTriangle );
2024 : : }
2025 : :
2026 : 1775 : return MB_SUCCESS;
2027 : : }
2028 : :
2029 : : // build the list of intx points and entities from database involved
2030 : : // vertices, edges, triangles
2031 : : // it could be just a list of vertices (easiest case to handle after)
2032 : :
2033 : 16 : ErrorCode FBEngine::compute_intersection_points( EntityHandle&, EntityHandle from, EntityHandle to, CartVect& Dir,
2034 : : std::vector< CartVect >& points, std::vector< EntityHandle >& entities,
2035 : : std::vector< EntityHandle >& triangles )
2036 : : {
2037 : : // keep a stack of triangles to process, and do not add those already processed
2038 : : // either mark them, or maybe keep them in a local set?
2039 : : // from and to are now nodes, start from them
2040 [ + - ][ + - ]: 16 : CartVect p1, p2; // the position of from and to
2041 [ + - ]: 16 : ErrorCode rval = _mbImpl->get_coords( &from, 1, (double*)&p1 );
2042 [ - + ][ # # ]: 16 : MBERRORR( rval, "failed to get 'from' coordinates" );
[ # # ]
2043 [ + - ]: 16 : rval = _mbImpl->get_coords( &to, 1, (double*)&p2 );
2044 [ - + ][ # # ]: 16 : MBERRORR( rval, "failed to get 'from' coordinates" );
[ # # ]
2045 : :
2046 [ + - ]: 16 : CartVect vect( p2 - p1 );
2047 [ + - ]: 16 : double dist2 = vect.length();
2048 [ - + ]: 16 : if( dist2 < tolerance_segment )
2049 : : {
2050 : : // we are done, return
2051 : 0 : return MB_SUCCESS;
2052 : : }
2053 [ + - ]: 16 : CartVect normPlane = Dir * vect;
2054 [ + - ]: 16 : normPlane.normalize();
2055 [ + - ]: 16 : std::set< EntityHandle > visitedTriangles;
2056 : 16 : CartVect currentPoint = p1;
2057 : : // push the current point if it is empty
2058 [ + - ]: 16 : if( points.size() == 0 )
2059 : : {
2060 [ + - ]: 16 : points.push_back( p1 );
2061 [ + - ]: 16 : entities.push_back( from ); // this is a node now
2062 : : }
2063 : :
2064 : : // these will be used in many loops
2065 : 16 : CartVect intx = p1; // somewhere to start
2066 : 16 : double param = -1.;
2067 : :
2068 : : // first intersection
2069 : 16 : EntityHandle currentBoundary = from; // it is a node, in the beginning
2070 : :
2071 [ + - ]: 16 : vect = p2 - currentPoint;
2072 [ + - ][ + + ]: 1791 : while( vect.length() > 0. )
2073 : : {
2074 : : // advance towards "to" node, from boundary handle
2075 [ + - ]: 1775 : EntityType etype = _mbImpl->type_from_handle( currentBoundary );
2076 : : // if vertex, look for other triangles connected which intersect our plane (defined by p1,
2077 : : // p2, dir)
2078 [ + - ]: 1775 : std::vector< EntityHandle > adj_tri;
2079 [ + - ]: 1775 : rval = _mbImpl->get_adjacencies( ¤tBoundary, 1, 2, false, adj_tri );
2080 : 1775 : unsigned int j = 0;
2081 : : EntityHandle tri;
2082 [ + + ]: 5330 : for( ; j < adj_tri.size(); j++ )
2083 : : {
2084 [ + - ]: 3602 : tri = adj_tri[j];
2085 [ + - ][ + - ]: 3602 : if( visitedTriangles.find( tri ) != visitedTriangles.end() )
[ + + ]
2086 : 1770 : continue; // get another triangle, this one was already visited
2087 : : // check if it is one of the triangles that was pierced already
2088 [ + - ][ + - ]: 1856 : if( _piercedTriangles.find( tri ) != _piercedTriangles.end() ) continue;
[ + - ][ + + ]
2089 : : // if vertex, look for opposite edge
2090 : : // if edge, look for 2 opposite edges
2091 : : // get vertices
2092 : : int nnodes;
2093 : : const EntityHandle* conn3;
2094 [ + - ]: 1832 : rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
2095 [ - + ][ # # ]: 1832 : MBERRORR( rval, "Failed to get connectivity" );
[ # # ]
2096 : : // if one of the nodes is to, stop right there
2097 : : {
2098 [ + + ][ + + ]: 1832 : if( conn3[0] == to || conn3[1] == to || conn3[2] == to )
[ + + ]
2099 : : {
2100 [ + - ]: 16 : visitedTriangles.insert( tri );
2101 [ + - ]: 16 : triangles.push_back( tri );
2102 : 16 : currentPoint = p2;
2103 [ + - ]: 16 : points.push_back( p2 );
2104 [ + - ]: 16 : entities.push_back( to ); // we are done
2105 : 47 : break; // this is break from for loop, we still need to get out of while
2106 : : // we will get out, because vect will become 0, (p2-p2)
2107 : : }
2108 : : }
2109 : : EntityHandle nn2[2];
2110 [ + + ]: 1816 : if( MBVERTEX == etype )
2111 : : {
2112 : 88 : nn2[0] = conn3[0];
2113 : 88 : nn2[1] = conn3[1];
2114 [ + + ]: 88 : if( nn2[0] == currentBoundary ) nn2[0] = conn3[2];
2115 [ + + ]: 88 : if( nn2[1] == currentBoundary ) nn2[1] = conn3[2];
2116 : : // get coordinates
2117 [ + - ][ + + ]: 264 : CartVect Pt[2];
2118 : :
2119 [ + - ]: 88 : rval = _mbImpl->get_coords( nn2, 2, (double*)&Pt[0] );
2120 [ - + ][ # # ]: 88 : MBERRORR( rval, "Failed to get coordinates" );
[ # # ]
2121 : : // see the intersection
2122 [ + - ][ + + ]: 88 : if( intersect_segment_and_plane_slice( Pt[0], Pt[1], currentPoint, p2, Dir, normPlane, intx, param ) )
2123 : : {
2124 : : // we should stop for loop, and decide if it is edge or vertex
2125 [ - + ]: 31 : if( param == 0. )
2126 : 0 : currentBoundary = nn2[0];
2127 : : else
2128 : : {
2129 [ - + ]: 31 : if( param == 1. )
2130 : 0 : currentBoundary = nn2[1];
2131 : : else // param between 0 and 1, so edge
2132 : : {
2133 : : // find the edge between vertices
2134 [ + - ]: 31 : std::vector< EntityHandle > edges1;
2135 : : // change the create flag to true, because that edge must exist in
2136 : : // current triangle if we want to advance; nn2 are 2 nodes in current
2137 : : // triangle!!
2138 [ + - ]: 31 : rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
2139 [ - + ][ # # ]: 31 : MBERRORR( rval, "Failed to get edges" );
[ # # ]
2140 [ - + ][ # # ]: 31 : if( edges1.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
[ # # ]
2141 [ + - ][ + - ]: 31 : currentBoundary = edges1[0];
2142 : : }
2143 : : }
2144 [ + - ]: 31 : visitedTriangles.insert( tri );
2145 : 31 : currentPoint = intx;
2146 [ + - ]: 31 : points.push_back( intx );
2147 [ + - ]: 31 : entities.push_back( currentBoundary );
2148 [ + - ]: 31 : triangles.push_back( tri );
2149 [ - + ]: 31 : if( debug_splits )
2150 [ # # ][ # # ]: 0 : std::cout << "vtx new tri : " << _mbImpl->id_from_handle( tri )
[ # # ]
2151 [ # # ][ # # ]: 0 : << " type bdy:" << _mbImpl->type_from_handle( currentBoundary ) << "\n";
[ # # ][ # # ]
2152 : 88 : break; // out of for loop over triangles
2153 : : }
2154 : : }
2155 : : else // this is MBEDGE, we have the other triangle to try out
2156 : : {
2157 : : // first find the nodes from existing boundary
2158 : : int nnodes2;
2159 : : const EntityHandle* conn2;
2160 [ + - ]: 1728 : rval = _mbImpl->get_connectivity( currentBoundary, conn2, nnodes2 );
2161 [ - + ][ # # ]: 3513 : MBERRORR( rval, "Failed to get connectivity" );
[ # # ]
2162 : 1728 : int thirdIndex = -1;
2163 [ + - ]: 3493 : for( int tj = 0; tj < 3; tj++ )
2164 : : {
2165 [ + + ][ + + ]: 3493 : if( ( conn3[tj] != conn2[0] ) && ( conn3[tj] != conn2[1] ) )
2166 : : {
2167 : 1728 : thirdIndex = tj;
2168 : 1728 : break;
2169 : : }
2170 : : }
2171 [ - + ][ # # ]: 1728 : if( -1 == thirdIndex ) MBERRORR( MB_FAILURE, " can't get third node" );
[ # # ]
2172 [ + - ][ + + ]: 6912 : CartVect Pt[3];
2173 [ + - ]: 1728 : rval = _mbImpl->get_coords( conn3, 3, (double*)&Pt[0] );
2174 [ - + ][ # # ]: 1728 : MBERRORR( rval, "Failed to get coords" );
[ # # ]
2175 : 1728 : int indexFirst = ( thirdIndex + 1 ) % 3;
2176 : 1728 : int indexSecond = ( thirdIndex + 2 ) % 3;
2177 : 1728 : int index[2] = { indexFirst, indexSecond };
2178 [ + - ]: 4321 : for( int k = 0; k < 2; k++ )
2179 : : {
2180 : 2593 : nn2[0] = conn3[index[k]], nn2[1] = conn3[thirdIndex];
2181 [ + - ][ + + ]: 2593 : if( intersect_segment_and_plane_slice( Pt[index[k]], Pt[thirdIndex], currentPoint, p2, Dir,
2182 : 5186 : normPlane, intx, param ) )
2183 : : {
2184 : : // we should stop for loop, and decide if it is edge or vertex
2185 [ - + ]: 1728 : if( param == 0. )
2186 : 0 : currentBoundary = conn3[index[k]]; // it is not really possible
2187 : : else
2188 : : {
2189 [ + + ]: 1728 : if( param == 1. )
2190 : 15 : currentBoundary = conn3[thirdIndex];
2191 : : else // param between 0 and 1, so edge is fine
2192 : : {
2193 : : // find the edge between vertices
2194 [ + - ]: 1713 : std::vector< EntityHandle > edges1;
2195 : : // change the create flag to true, because that edge must exist in
2196 : : // current triangle if we want to advance; nn2 are 2 nodes in
2197 : : // current triangle!!
2198 [ + - ]: 1713 : rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
2199 [ - + ][ # # ]: 1713 : MBERRORR( rval, "Failed to get edges" );
[ # # ]
2200 [ - + ]: 1713 : if( edges1.size() != 1 )
2201 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
2202 [ + - ][ + - ]: 1713 : currentBoundary = edges1[0];
2203 : : }
2204 : : }
2205 [ + - ]: 1728 : visitedTriangles.insert( tri );
2206 : 1728 : currentPoint = intx;
2207 [ + - ]: 1728 : points.push_back( intx );
2208 [ + - ]: 1728 : entities.push_back( currentBoundary );
2209 [ + - ]: 1728 : triangles.push_back( tri );
2210 [ - + ]: 1728 : if( debug_splits )
2211 : : {
2212 [ # # ][ # # ]: 0 : std::cout << "edge new tri : " << _mbImpl->id_from_handle( tri )
[ # # ]
2213 [ # # ][ # # ]: 0 : << " type bdy: " << _mbImpl->type_from_handle( currentBoundary )
[ # # ]
2214 [ # # ][ # # ]: 0 : << " id: " << _mbImpl->id_from_handle( currentBoundary ) << "\n";
[ # # ][ # # ]
2215 [ # # ]: 0 : _mbImpl->list_entity( currentBoundary );
2216 : : }
2217 : 1728 : break; // out of for loop over triangles
2218 : : }
2219 : : // we should not reach here
2220 : : }
2221 : : }
2222 : : }
2223 : : /*if (j==adj_tri.size())
2224 : : {
2225 : : MBERRORR(MB_FAILURE, "did not advance");
2226 : : }*/
2227 [ + - ][ + - ]: 1775 : vect = p2 - currentPoint;
2228 : 1775 : }
2229 : :
2230 [ - + ]: 16 : if( debug_splits )
2231 [ # # ][ # # ]: 0 : std::cout << "nb entities: " << entities.size() << " triangles:" << triangles.size()
[ # # ][ # # ]
2232 [ # # ][ # # ]: 0 : << " points.size(): " << points.size() << "\n";
[ # # ]
2233 : :
2234 : 16 : return MB_SUCCESS;
2235 : : }
2236 : :
2237 : 0 : ErrorCode FBEngine::split_edge_at_point( EntityHandle edge, CartVect& point, EntityHandle& new_edge )
2238 : : {
2239 : : // return MB_NOT_IMPLEMENTED;
2240 : : // first, we need to find the closest point on the smooth edge, then
2241 : : // split the smooth edge, then call the split_edge_at_mesh_node
2242 : : // or maybe just find the closest node??
2243 : : // first of all, we need to find a point on the smooth edge, close by
2244 : : // then split the mesh edge (if needed)
2245 : : // this would be quite a drama, as the new edge has to be inserted in
2246 : : // order for proper geo edge definition
2247 : :
2248 : : // first of all, find a closest point
2249 : : // usually called for
2250 [ # # ]: 0 : if( debug_splits )
2251 [ # # ][ # # ]: 0 : { std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) << " at point:" << point << "\n"; }
[ # # ][ # # ]
[ # # ][ # # ]
2252 [ # # ]: 0 : int dim = _my_geomTopoTool->dimension( edge );
2253 [ # # ]: 0 : if( dim != 1 ) return MB_FAILURE;
2254 [ # # ]: 0 : if( !_smooth ) return MB_FAILURE; // call it only for smooth option...
2255 : : // maybe we should do something for linear too
2256 : :
2257 [ # # ]: 0 : SmoothCurve* curve = this->_edges[edge];
2258 : : EntityHandle closeNode;
2259 : : int edgeIndex;
2260 [ # # ][ # # ]: 0 : double u = curve->u_from_position( point[0], point[1], point[2], closeNode, edgeIndex );
[ # # ][ # # ]
2261 [ # # ]: 0 : if( 0 == closeNode )
2262 : : {
2263 : : // we really need to split an existing edge
2264 : : // do not do that yet
2265 : : // evaluate from u:
2266 : : /*double pos[3];
2267 : : curve->position_from_u(u, pos[0], pos[1], pos[2] );*/
2268 : : // create a new node here, and split one edge
2269 : : // change also connectivity, create new triangles on both sides ...
2270 [ # # ][ # # ]: 0 : std::cout << "not found a close node, u is: " << u << " edge index: " << edgeIndex << "\n";
[ # # ][ # # ]
[ # # ]
2271 : 0 : return MB_FAILURE; // not implemented yet
2272 : : }
2273 [ # # ]: 0 : return split_edge_at_mesh_node( edge, closeNode, new_edge );
2274 : : }
2275 : :
2276 : 0 : ErrorCode FBEngine::split_edge_at_mesh_node( EntityHandle edge, EntityHandle node, EntityHandle& new_edge )
2277 : : {
2278 : : // the node should be in the list of nodes
2279 : :
2280 [ # # ]: 0 : int dim = _my_geomTopoTool->dimension( edge );
2281 [ # # ]: 0 : if( dim != 1 ) return MB_FAILURE; // not an edge
2282 : :
2283 [ # # ]: 0 : if( debug_splits )
2284 : : {
2285 [ # # ][ # # ]: 0 : std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
[ # # ]
2286 [ # # ][ # # ]: 0 : << " with global id: " << _my_geomTopoTool->global_id( edge )
[ # # ]
2287 [ # # ][ # # ]: 0 : << " at node:" << _mbImpl->id_from_handle( node ) << "\n";
[ # # ][ # # ]
2288 : : }
2289 : :
2290 : : // now get the edges
2291 : : // the order is important...
2292 : : // these are ordered sets !!
2293 [ # # ]: 0 : std::vector< EntityHandle > ents;
2294 [ # # ]: 0 : ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
2295 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
2296 [ # # ]: 0 : if( ents.size() < 1 ) return MB_FAILURE; // no edges
2297 : :
2298 : 0 : const EntityHandle* conn = NULL;
2299 : : int len;
2300 : : // find the edge connected to the splitting node
2301 : 0 : int num_mesh_edges = (int)ents.size();
2302 : : int index_edge;
2303 : 0 : EntityHandle firstNode = 0;
2304 [ # # ]: 0 : for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
2305 : : {
2306 [ # # ][ # # ]: 0 : rval = MBI->get_connectivity( ents[index_edge], conn, len );
2307 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
2308 [ # # ]: 0 : if( index_edge == 0 ) firstNode = conn[0]; // will be used to decide vertex sets adjacencies
2309 [ # # ]: 0 : if( conn[0] == node )
2310 : : {
2311 [ # # ]: 0 : if( index_edge == 0 )
2312 : : {
2313 : 0 : new_edge = 0; // no need to split, it is the first node
2314 : 0 : return MB_SUCCESS; // no split
2315 : : }
2316 : : else
2317 : 0 : return MB_FAILURE; // we should have found the node already , wrong
2318 : : // connectivity
2319 : : }
2320 [ # # ]: 0 : else if( conn[1] == node )
2321 : : {
2322 : : // we found the index of interest
2323 : 0 : break;
2324 : : }
2325 : : }
2326 [ # # ]: 0 : if( index_edge == num_mesh_edges - 1 )
2327 : : {
2328 : 0 : new_edge = 0; // no need to split, it is the last node
2329 : 0 : return MB_SUCCESS; // no split
2330 : : }
2331 : :
2332 : : // here, we have 0 ... index_edge edges in the first set,
2333 : : // create a vertex set and add the node to it
2334 : :
2335 [ # # ]: 0 : if( debug_splits )
2336 [ # # ][ # # ]: 0 : { std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n"; }
[ # # ][ # # ]
[ # # ]
2337 : :
2338 : : // at this moment, the node set should have been already created in new_geo_edge
2339 : : EntityHandle nodeSet; // the node set that has the node (find it!)
2340 [ # # ]: 0 : bool found = find_vertex_set_for_node( node, nodeSet );
2341 : :
2342 [ # # ]: 0 : if( !found )
2343 : : {
2344 : : // create a node set and add the node
2345 : :
2346 : : // must be an error, but create one nevertheless
2347 [ # # ]: 0 : rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
2348 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to create a new vertex set" );
[ # # ]
2349 : :
2350 [ # # ]: 0 : rval = _mbImpl->add_entities( nodeSet, &node, 1 );
2351 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to add the node to the set" );
[ # # ]
2352 : :
2353 [ # # ]: 0 : rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
2354 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to commit the node set" );
[ # # ]
2355 : :
2356 [ # # ]: 0 : if( debug_splits )
2357 : : {
2358 [ # # ]: 0 : std::cout << " create a vertex set (this should have been created before!)"
2359 [ # # ][ # # ]: 0 : << _mbImpl->id_from_handle( nodeSet )
2360 [ # # ][ # # ]: 0 : << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
[ # # ][ # # ]
2361 : : }
2362 : : }
2363 : :
2364 : : // we need to remove the remaining mesh edges from first set, and add it
2365 : : // to the second set, in order
2366 : :
2367 [ # # ]: 0 : rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
2368 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't create geo edge" );
[ # # ]
2369 : :
2370 : 0 : int remaining = num_mesh_edges - 1 - index_edge;
2371 : :
2372 : : // add edges to the edge set
2373 [ # # ][ # # ]: 0 : rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
2374 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't add edges to the new edge" );
[ # # ]
2375 : :
2376 : : // also, remove the second node set from old edge
2377 : : // remove the edges index_edge+1 and up
2378 : :
2379 [ # # ][ # # ]: 0 : rval = _mbImpl->remove_entities( edge, &ents[index_edge + 1], remaining );
2380 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't remove edges from the old edge" );
[ # # ]
2381 : :
2382 : : // need to find the adjacent vertex sets
2383 [ # # ]: 0 : Range vertexRange;
2384 [ # # ]: 0 : rval = getAdjacentEntities( edge, 0, vertexRange );
2385 : :
2386 : : EntityHandle secondSet;
2387 [ # # ][ # # ]: 0 : if( vertexRange.size() == 1 )
2388 : : {
2389 : : // initially a periodic edge, OK to add the new set to both edges, and the
2390 : : // second set
2391 [ # # ]: 0 : secondSet = vertexRange[0];
2392 : : }
2393 : : else
2394 : : {
2395 [ # # ][ # # ]: 0 : if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices
2396 : : // find first node
2397 : : int k;
2398 [ # # ]: 0 : for( k = 0; k < 2; k++ )
2399 : : {
2400 [ # # ]: 0 : Range verts;
2401 [ # # ][ # # ]: 0 : rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );
2402 : :
2403 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't get vertices from vertex set" );
[ # # ]
2404 [ # # ][ # # ]: 0 : if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
[ # # ][ # # ]
2405 [ # # ][ # # ]: 0 : if( firstNode == verts[0] )
2406 : : {
2407 [ # # ]: 0 : secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0
2408 [ # # # ]: 0 : break;
2409 : : }
2410 : 0 : }
2411 [ # # ]: 0 : if( k >= 2 )
2412 : : {
2413 : : // it means we didn't find the right node set
2414 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, " can't find the right vertex" );
2415 : : }
2416 : : // remove the second set from the connectivity with the
2417 : : // edge (parent-child relation)
2418 : : // remove_parent_child( EntityHandle parent,
2419 : : // EntityHandle child )
2420 [ # # ]: 0 : rval = _mbImpl->remove_parent_child( edge, secondSet );
2421 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't remove the second vertex from edge" );
[ # # ]
2422 : : }
2423 : : // at this point, we just need to add to both edges the new node set (vertex)
2424 [ # # ]: 0 : rval = _mbImpl->add_parent_child( edge, nodeSet );
2425 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't add new vertex to old edge" );
[ # # ]
2426 : :
2427 [ # # ]: 0 : rval = _mbImpl->add_parent_child( new_edge, nodeSet );
2428 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't add new vertex to new edge" );
[ # # ]
2429 : :
2430 : : // one thing that I forgot: add the secondSet as a child to new edge!!!
2431 : : // (so far, the new edge has only one end vertex!)
2432 [ # # ]: 0 : rval = _mbImpl->add_parent_child( new_edge, secondSet );
2433 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't add second vertex to new edge" );
[ # # ]
2434 : :
2435 : : // now, add the edge and vertex to geo tool
2436 : :
2437 [ # # ]: 0 : rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
2438 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't add new edge" );
[ # # ]
2439 : :
2440 : : // next, get the adjacent faces to initial edge, and add them as parents
2441 : : // to the new edge
2442 : :
2443 : : // need to find the adjacent face sets
2444 [ # # ]: 0 : Range faceRange;
2445 [ # # ]: 0 : rval = getAdjacentEntities( edge, 2, faceRange );
2446 : :
2447 : : // these faces will be adjacent to the new edge too!
2448 : : // need to set the senses of new edge within faces
2449 : :
2450 [ # # ][ # # ]: 0 : for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
[ # # ][ # # ]
[ # # ]
2451 : : {
2452 [ # # ]: 0 : EntityHandle face = *it;
2453 [ # # ]: 0 : rval = _mbImpl->add_parent_child( face, new_edge );
2454 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't add new edge - face parent relation" );
[ # # ]
2455 : : int sense;
2456 [ # # ]: 0 : rval = _my_geomTopoTool->get_sense( edge, face, sense );
2457 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
[ # # ]
2458 : : // keep the same sense for the new edge within the faces
2459 [ # # ]: 0 : rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
2460 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
[ # # ]
2461 : : }
2462 : :
2463 : 0 : return MB_SUCCESS;
2464 : : }
2465 : :
2466 : 2 : ErrorCode FBEngine::split_bedge_at_new_mesh_node( EntityHandle edge, EntityHandle node, EntityHandle brokenEdge,
2467 : : EntityHandle& new_edge )
2468 : : {
2469 : : // the node should be in the list of nodes
2470 : :
2471 [ + - ]: 2 : int dim = _my_geomTopoTool->dimension( edge );
2472 [ - + ]: 2 : if( dim != 1 ) return MB_FAILURE; // not an edge
2473 : :
2474 [ - + ]: 2 : if( debug_splits )
2475 : : {
2476 [ # # ][ # # ]: 0 : std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
[ # # ]
2477 [ # # ][ # # ]: 0 : << " with global id: " << _my_geomTopoTool->global_id( edge )
[ # # ]
2478 [ # # ][ # # ]: 0 : << " at new node:" << _mbImpl->id_from_handle( node ) << "breaking mesh edge"
[ # # ][ # # ]
2479 [ # # ][ # # ]: 0 : << _mbImpl->id_from_handle( brokenEdge ) << "\n";
[ # # ]
2480 : : }
2481 : :
2482 : : // now get the edges
2483 : : // the order is important...
2484 : : // these are ordered sets !!
2485 [ + - ]: 2 : std::vector< EntityHandle > ents;
2486 [ + - ]: 2 : ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
2487 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
2488 [ - + ]: 2 : if( ents.size() < 1 ) return MB_FAILURE; // no edges
2489 : :
2490 : 2 : const EntityHandle* conn = NULL;
2491 : : int len;
2492 : : // the edge connected to the splitting node is brokenEdge
2493 : : // find the small edges it is broken into, which are connected to
2494 : : // new vertex (node) and its ends; also, if necessary, reorient these small edges
2495 : : // for proper orientation
2496 [ + - ]: 2 : rval = _mbImpl->get_connectivity( brokenEdge, conn, len );
2497 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to get connectivity of broken edge" );
[ # # ]
2498 : :
2499 : : // we already created the new edges, make sure they are oriented fine; if not, revert them
2500 : 2 : EntityHandle conn02[] = { conn[0], node }; // first node and new node
2501 : : // default, intersect
2502 [ + - ]: 4 : std::vector< EntityHandle > adj_edges;
2503 [ + - ]: 2 : rval = _mbImpl->get_adjacencies( conn02, 2, 1, false, adj_edges );
2504 [ + - ][ - + ]: 2 : if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find small split edge" );
[ - + ][ # # ]
[ # # ]
2505 : :
2506 : : // get this edge connectivity;
2507 [ + - ]: 2 : EntityHandle firstEdge = adj_edges[0];
2508 : 2 : const EntityHandle* connActual = NULL;
2509 [ + - ]: 2 : rval = _mbImpl->get_connectivity( firstEdge, connActual, len );
2510 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to get connectivity of first split edge" );
[ # # ]
2511 : : // if it is the same as conn02, we are happy
2512 [ - + ]: 2 : if( conn02[0] != connActual[0] )
2513 : : {
2514 : : // reset connectivity of edge
2515 [ # # ]: 0 : rval = _mbImpl->set_connectivity( firstEdge, conn02, 2 );
2516 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to reset connectivity of first split edge" );
[ # # ]
2517 : : }
2518 : : // now treat the second edge
2519 : 2 : adj_edges.clear(); //
2520 : 2 : EntityHandle conn21[] = { node, conn[1] }; // new node and second node
2521 [ + - ]: 2 : rval = _mbImpl->get_adjacencies( conn21, 2, 1, false, adj_edges );
2522 [ + - ][ - + ]: 2 : if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find second small split edge" );
[ - + ][ # # ]
[ # # ]
2523 : :
2524 : : // get this edge connectivity;
2525 [ + - ]: 2 : EntityHandle secondEdge = adj_edges[0];
2526 [ + - ]: 2 : rval = _mbImpl->get_connectivity( firstEdge, connActual, len );
2527 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to get connectivity of first split edge" );
[ # # ]
2528 : : // if it is the same as conn21, we are happy
2529 [ + - ]: 2 : if( conn21[0] != connActual[0] )
2530 : : {
2531 : : // reset connectivity of edge
2532 [ + - ]: 2 : rval = _mbImpl->set_connectivity( secondEdge, conn21, 2 );
2533 [ - + ][ # # ]: 2 : MBERRORR( rval, "Failed to reset connectivity of first split edge" );
[ # # ]
2534 : : }
2535 : :
2536 : 2 : int num_mesh_edges = (int)ents.size();
2537 : : int index_edge; // this is the index of the edge that will be removed (because it is split)
2538 : : // the rest of edges will be put in order in the (remaining) edge and new edge
2539 : :
2540 [ + - ]: 251 : for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
2541 [ + - ][ + + ]: 251 : if( brokenEdge == ents[index_edge] ) break;
2542 [ - + ][ # # ]: 2 : if( index_edge >= num_mesh_edges ) MBERRORR( MB_FAILURE, "can't find the broken edge" );
[ # # ]
2543 : :
2544 : : // so the edges up to index_edge and firstEdge, will form the "old" edge
2545 : : // the edges secondEdge and from index_edge+1 to end will form the new_edge
2546 : :
2547 : : // here, we have 0 ... index_edge edges in the first set,
2548 : : // create a vertex set and add the node to it
2549 : :
2550 [ - + ]: 2 : if( debug_splits )
2551 [ # # ][ # # ]: 0 : { std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n"; }
[ # # ][ # # ]
[ # # ]
2552 : :
2553 : : // at this moment, the node set should have been already created in new_geo_edge
2554 : : EntityHandle nodeSet; // the node set that has the node (find it!)
2555 [ + - ]: 2 : bool found = find_vertex_set_for_node( node, nodeSet );
2556 : :
2557 [ - + ]: 2 : if( !found )
2558 : : {
2559 : : // create a node set and add the node
2560 : :
2561 : : // must be an error, but create one nevertheless
2562 [ # # ]: 0 : rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
2563 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to create a new vertex set" );
[ # # ]
2564 : :
2565 [ # # ]: 0 : rval = _mbImpl->add_entities( nodeSet, &node, 1 );
2566 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to add the node to the set" );
[ # # ]
2567 : :
2568 [ # # ]: 0 : rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
2569 [ # # ][ # # ]: 0 : MBERRORR( rval, "Failed to commit the node set" );
[ # # ]
2570 : :
2571 [ # # ]: 0 : if( debug_splits )
2572 : : {
2573 [ # # ]: 0 : std::cout << " create a vertex set (this should have been created before!)"
2574 [ # # ][ # # ]: 0 : << _mbImpl->id_from_handle( nodeSet )
2575 [ # # ][ # # ]: 0 : << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
[ # # ][ # # ]
2576 : : }
2577 : : }
2578 : :
2579 : : // we need to remove the remaining mesh edges from first set, and add it
2580 : : // to the second set, in order
2581 : :
2582 [ + - ]: 2 : rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
2583 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't create geo edge" );
[ # # ]
2584 : :
2585 : 2 : int remaining = num_mesh_edges - 1 - index_edge;
2586 : :
2587 : : // add edges to the new edge set
2588 [ + - ]: 2 : rval = _mbImpl->add_entities( new_edge, &secondEdge, 1 ); // add the second split edge to new edge
2589 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't add second split edge to the new edge" );
[ # # ]
2590 : : // then add the rest
2591 [ + - ][ + - ]: 2 : rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
2592 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't add edges to the new edge" );
[ # # ]
2593 : :
2594 : : // also, remove the second node set from old edge
2595 : : // remove the edges index_edge and up
2596 : :
2597 [ + - ][ + - ]: 2 : rval = _mbImpl->remove_entities( edge, &ents[index_edge], remaining + 1 ); // include the
2598 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove edges from the old edge" );
[ # # ]
2599 : : // add the firstEdge too
2600 [ + - ]: 2 : rval = _mbImpl->add_entities( edge, &firstEdge, 1 ); // add the second split edge to new edge
2601 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't add first split edge to the old edge" );
[ # # ]
2602 : :
2603 : : // need to find the adjacent vertex sets
2604 [ + - ]: 4 : Range vertexRange;
2605 [ + - ]: 2 : rval = getAdjacentEntities( edge, 0, vertexRange );
2606 : :
2607 : : EntityHandle secondSet;
2608 [ + - ][ + + ]: 2 : if( vertexRange.size() == 1 )
2609 : : {
2610 : : // initially a periodic edge, OK to add the new set to both edges, and the
2611 : : // second set
2612 [ + - ]: 1 : secondSet = vertexRange[0];
2613 : : }
2614 : : else
2615 : : {
2616 [ + - ][ - + ]: 1 : if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices
2617 : : // find first node
2618 : : EntityHandle firstNode;
2619 : :
2620 [ + - ][ + - ]: 1 : rval = MBI->get_connectivity( ents[0], conn, len );
2621 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
2622 : 1 : firstNode = conn[0]; // this is the first node of the initial edge
2623 : : // we will use it to identify the vertex set associated with it
2624 : : int k;
2625 [ + - ]: 2 : for( k = 0; k < 2; k++ )
2626 : : {
2627 [ + - ]: 1 : Range verts;
2628 [ + - ][ + - ]: 1 : rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );
2629 : :
2630 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get vertices from vertex set" );
[ # # ]
2631 [ + - ][ - + ]: 1 : if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
[ # # ][ # # ]
2632 [ + - ][ + - ]: 1 : if( firstNode == verts[0] )
2633 : : {
2634 [ + - ]: 1 : secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0
2635 [ - - + ]: 1 : break;
2636 : : }
2637 : 0 : }
2638 [ - + ]: 1 : if( k >= 2 )
2639 : : {
2640 : : // it means we didn't find the right node set
2641 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, " can't find the right vertex" );
2642 : : }
2643 : : // remove the second set from the connectivity with the
2644 : : // edge (parent-child relation)
2645 : : // remove_parent_child( EntityHandle parent,
2646 : : // EntityHandle child )
2647 [ + - ]: 1 : rval = _mbImpl->remove_parent_child( edge, secondSet );
2648 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't remove the second vertex from edge" );
[ # # ]
2649 : : }
2650 : : // at this point, we just need to add to both edges the new node set (vertex)
2651 [ + - ]: 2 : rval = _mbImpl->add_parent_child( edge, nodeSet );
2652 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't add new vertex to old edge" );
[ # # ]
2653 : :
2654 [ + - ]: 2 : rval = _mbImpl->add_parent_child( new_edge, nodeSet );
2655 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't add new vertex to new edge" );
[ # # ]
2656 : :
2657 : : // one thing that I forgot: add the secondSet as a child to new edge!!!
2658 : : // (so far, the new edge has only one end vertex!)
2659 [ + - ]: 2 : rval = _mbImpl->add_parent_child( new_edge, secondSet );
2660 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't add second vertex to new edge" );
[ # # ]
2661 : :
2662 : : // now, add the edge and vertex to geo tool
2663 : :
2664 [ + - ]: 2 : rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
2665 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't add new edge" );
[ # # ]
2666 : :
2667 : : // next, get the adjacent faces to initial edge, and add them as parents
2668 : : // to the new edge
2669 : :
2670 : : // need to find the adjacent face sets
2671 [ + - ]: 4 : Range faceRange;
2672 [ + - ]: 2 : rval = getAdjacentEntities( edge, 2, faceRange );
2673 : :
2674 : : // these faces will be adjacent to the new edge too!
2675 : : // need to set the senses of new edge within faces
2676 : :
2677 [ + - ][ + - ]: 4 : for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
[ + - ][ + - ]
[ + + ]
2678 : : {
2679 [ + - ]: 2 : EntityHandle face = *it;
2680 [ + - ]: 2 : rval = _mbImpl->add_parent_child( face, new_edge );
2681 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't add new edge - face parent relation" );
[ # # ]
2682 : : int sense;
2683 [ + - ]: 2 : rval = _my_geomTopoTool->get_sense( edge, face, sense );
2684 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
[ # # ]
2685 : : // keep the same sense for the new edge within the faces
2686 [ + - ]: 2 : rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
2687 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
[ # # ]
2688 : : }
2689 : :
2690 : 4 : return MB_SUCCESS;
2691 : : }
2692 : :
2693 : 2 : ErrorCode FBEngine::split_boundary( EntityHandle face, EntityHandle atNode )
2694 : : {
2695 : : // find the boundary edges, and split the one that we find it is a part of
2696 [ - + ]: 2 : if( debug_splits )
2697 : : {
2698 [ # # ][ # # ]: 0 : std::cout << "Split face " << _mbImpl->id_from_handle( face )
[ # # ]
2699 [ # # ][ # # ]: 0 : << " at node:" << _mbImpl->id_from_handle( atNode ) << "\n";
[ # # ][ # # ]
2700 : : }
2701 [ + - ]: 2 : Range bound_edges;
2702 [ + - ]: 2 : ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
2703 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't get boundary edges" );
[ # # ]
2704 [ + - ][ + - ]: 2 : bool brokEdge = _brokenEdges.find( atNode ) != _brokenEdges.end();
2705 : :
2706 [ + - ][ # # ]: 4 : for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
[ + - ][ + - ]
[ + - ]
2707 : : {
2708 [ + - ]: 2 : EntityHandle b_edge = *it;
2709 : : // get all edges in range
2710 [ + - ]: 2 : Range mesh_edges;
2711 [ + - ]: 2 : rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
2712 [ - + ][ # # ]: 2 : MBERRORR( rval, " can't get mesh edges" );
[ # # ]
2713 [ + - ]: 2 : if( brokEdge )
2714 : : {
2715 [ + - ]: 2 : EntityHandle brokenEdge = _brokenEdges[atNode];
2716 [ + - ][ + - ]: 2 : if( mesh_edges.find( brokenEdge ) != mesh_edges.end() )
[ + - ][ + - ]
2717 : : {
2718 : : EntityHandle new_edge;
2719 [ + - ]: 2 : rval = split_bedge_at_new_mesh_node( b_edge, atNode, brokenEdge, new_edge );
2720 : 2 : return rval;
2721 : : }
2722 : : }
2723 : : else
2724 : : {
2725 [ # # ]: 0 : Range nodes;
2726 [ # # ]: 0 : rval = _mbImpl->get_connectivity( mesh_edges, nodes );
2727 [ # # ][ # # ]: 0 : MBERRORR( rval, " can't get nodes from mesh edges" );
[ # # ]
2728 : :
2729 [ # # ][ # # ]: 0 : if( nodes.find( atNode ) != nodes.end() )
[ # # ][ # # ]
2730 : : {
2731 : : // we found our boundary edge candidate
2732 : : EntityHandle new_edge;
2733 [ # # ]: 0 : rval = split_edge_at_mesh_node( b_edge, atNode, new_edge );
2734 [ # # ]: 0 : return rval;
2735 [ - + ]: 2 : }
2736 : : }
2737 : 0 : }
2738 : : // if the node was not found in any "current" boundary, it broke an existing
2739 : : // boundary edge
2740 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, " we did not find an appropriate boundary edge" );
2741 : 2 : ; //
2742 : : }
2743 : :
2744 : 34 : bool FBEngine::find_vertex_set_for_node( EntityHandle iNode, EntityHandle& oVertexSet )
2745 : : {
2746 : 34 : bool found = false;
2747 [ + - ]: 34 : Range vertex_sets;
2748 : :
2749 : 34 : const int zero = 0;
2750 : 34 : const void* const zero_val[] = { &zero };
2751 : : Tag geom_tag;
2752 [ + - ]: 34 : ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
2753 [ - + ]: 34 : if( MB_SUCCESS != rval ) return false;
2754 [ + - ]: 34 : rval = _mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag, zero_val, 1, vertex_sets );
2755 [ - + ]: 34 : if( MB_SUCCESS != rval ) return false;
2756 : : // local _gsets, as we might have not updated the local lists
2757 : : // see if ends of geo edge generated is in a geo set 0 or not
2758 : :
2759 [ + - ][ + - ]: 172 : for( Range::iterator vsit = vertex_sets.begin(); vsit != vertex_sets.end(); ++vsit )
[ + - ][ + - ]
[ + + ]
2760 : : {
2761 [ + - ]: 155 : EntityHandle vset = *vsit;
2762 : : // is the node part of this set?
2763 [ + - ][ + + ]: 155 : if( _mbImpl->contains_entities( vset, &iNode, 1 ) )
2764 : : {
2765 : :
2766 : 17 : found = true;
2767 : 17 : oVertexSet = vset;
2768 : 17 : break;
2769 : : }
2770 : : }
2771 : 34 : return found;
2772 : : }
2773 : 1 : ErrorCode FBEngine::redistribute_boundary_edges_to_faces( EntityHandle face, EntityHandle newFace,
2774 : : std::vector< EntityHandle >& chainedEdges )
2775 : : {
2776 : :
2777 : : // so far, original boundary edges are all parent/child relations for face
2778 : : // we should get them all, and see if they are truly adjacent to face or newFace
2779 : : // decide also on the orientation/sense with respect to the triangles
2780 [ + - ]: 1 : Range r1; // range in old face
2781 [ + - ]: 2 : Range r2; // range of tris in new face
2782 [ + - ]: 1 : ErrorCode rval = _mbImpl->get_entities_by_dimension( face, 2, r1 );
2783 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get triangles from old face" );
[ # # ]
2784 [ + - ]: 1 : rval = _mbImpl->get_entities_by_dimension( newFace, 2, r2 );
2785 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get triangles from new face" );
[ # # ]
2786 : : // right now, all boundary edges are children of face
2787 : : // we need to get them all, and verify if they indeed are adjacent to face
2788 [ + - ]: 2 : Range children;
2789 [ + - ]: 1 : rval = _mbImpl->get_child_meshsets( face, children ); // only direct children are of interest
2790 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get children sets from face" );
[ # # ]
2791 : :
2792 [ + - ][ + - ]: 7 : for( Range::iterator it = children.begin(); it != children.end(); ++it )
[ + - ][ + - ]
[ + + ]
2793 : : {
2794 [ + - ]: 6 : EntityHandle edge = *it;
2795 [ + - ][ + - ]: 6 : if( std::find( chainedEdges.begin(), chainedEdges.end(), edge ) != chainedEdges.end() )
[ + + ]
2796 : 5 : continue; // we already set this one fine
2797 : : // get one mesh edge from the edge; we have to get all of them!!
2798 [ + - ][ - + ]: 3 : if( _my_geomTopoTool->dimension( edge ) != 1 ) continue; // not an edge
2799 [ + - ]: 3 : Range mesh_edges;
2800 [ + - ]: 3 : rval = _mbImpl->get_entities_by_handle( edge, mesh_edges );
2801 [ - + ][ # # ]: 3 : MBERRORR( rval, " can't get mesh edges from edge" );
[ # # ]
2802 [ + - ][ - + ]: 3 : if( mesh_edges.empty() ) MBERRORR( MB_FAILURE, " no mesh edges" );
[ # # ][ # # ]
2803 [ + - ]: 3 : EntityHandle mesh_edge = mesh_edges[0]; // first one is enough
2804 : : // get its triangles; see which one is in r1 or r2; it should not be in both
2805 [ + - ]: 6 : Range triangles;
[ + + - ]
2806 [ + - ]: 3 : rval = _mbImpl->get_adjacencies( &mesh_edge, 1, 2, false, triangles );
2807 [ - + ][ # # ]: 3 : MBERRORR( rval, " can't get adjacent triangles" );
[ # # ]
2808 [ + - ]: 6 : Range intx1 = intersect( triangles, r1 );
[ + - + ]
2809 [ + - ]: 6 : Range intx2 = intersect( triangles, r2 );
[ + - + ]
2810 [ + - ][ + + ]: 3 : if( !intx1.empty() && !intx2.empty() ) MBERRORR( MB_FAILURE, " at least one should be empty" );
[ + - ][ - + ]
[ - + ][ # # ]
[ # # ]
2811 : :
2812 [ + - ][ + + ]: 3 : if( intx2.empty() )
2813 : : {
2814 : : // it means it is in the range r1; the sense should be fine, no need to reset
2815 : : // the sense should have been fine, also
2816 : 2 : continue;
2817 : : }
2818 : : // so it must be a triangle in r2;
2819 [ + - ]: 1 : EntityHandle triangle = intx2[0]; // one triangle only
2820 : : // remove the edge from boundary of face, and add it to the boundary of newFace
2821 : : // remove_parent_child( EntityHandle parent, EntityHandle child )
2822 [ + - ]: 1 : rval = _mbImpl->remove_parent_child( face, edge );
2823 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't remove parent child relation for edge" );
[ # # ]
2824 : : // add to the new face
2825 [ + - ]: 1 : rval = _mbImpl->add_parent_child( newFace, edge );
2826 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't add parent child relation for edge" );
[ # # ]
2827 : :
2828 : : // set some sense, based on the sense of the mesh_edge in triangle
2829 : : int num1, sense, offset;
2830 [ + - ]: 1 : rval = _mbImpl->side_number( triangle, mesh_edge, num1, sense, offset );
2831 [ - + ][ # # ]: 1 : MBERRORR( rval, "mesh edge not adjacent to triangle" );
[ # # ]
2832 : :
2833 [ + - ]: 1 : rval = this->_my_geomTopoTool->set_sense( edge, newFace, sense );
2834 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't set proper sense of edge in face" );
[ # # ]
[ + - + ]
2835 : 1 : }
2836 : :
2837 : 2 : return MB_SUCCESS;
2838 : : }
2839 : :
2840 : 4 : ErrorCode FBEngine::set_neumann_tags( EntityHandle face, EntityHandle newFace )
2841 : : {
2842 : : // these are for debugging purposes only
2843 : : // check the initial tag, then
2844 : : Tag ntag;
2845 [ + - ]: 4 : ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag );
2846 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get tag handle" );
[ # # ]
2847 : : // check the value for face
2848 : : int nval;
2849 [ + - ]: 4 : rval = _mbImpl->tag_get_data( ntag, &face, 1, &nval );
2850 [ + - ]: 4 : if( MB_SUCCESS == rval ) { nval++; }
2851 : : else
2852 : : {
2853 : 0 : nval = 1;
2854 [ # # ]: 0 : rval = _mbImpl->tag_set_data( ntag, &face, 1, &nval );
2855 [ # # ][ # # ]: 0 : MBERRORR( rval, "can't set tag" );
[ # # ]
2856 : 0 : nval = 2;
2857 : : }
2858 [ + - ]: 4 : rval = _mbImpl->tag_set_data( ntag, &newFace, 1, &nval );
2859 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set tag" );
[ # # ]
2860 : :
2861 : 4 : return MB_SUCCESS;
2862 : : }
2863 : :
2864 : : // split the quads if needed; it will create a new gtt, which will
2865 : : // contain triangles instead of quads
2866 : 10 : ErrorCode FBEngine::split_quads()
2867 : : {
2868 : : // first see if we have any quads in the 2d faces
2869 : : // _my_gsets[2] is a range of surfaces (moab sets of dimension 2)
2870 : 10 : int num_quads = 0;
2871 [ + - ][ + - ]: 26 : for( Range::iterator it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it )
[ + - ][ + - ]
[ + + ]
2872 : : {
2873 [ + - ]: 16 : EntityHandle surface = *it;
2874 : 16 : int num_q = 0;
2875 [ + - ]: 16 : _mbImpl->get_number_entities_by_type( surface, MBQUAD, num_q );
2876 : 16 : num_quads += num_q;
2877 : : }
2878 : :
2879 [ + + ]: 10 : if( num_quads == 0 ) return MB_SUCCESS; // nothing to do
2880 : :
2881 : 1 : GeomTopoTool* new_gtt = NULL;
2882 [ + - ]: 1 : ErrorCode rval = _my_geomTopoTool->duplicate_model( new_gtt );
2883 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't duplicate model" );
[ # # ]
2884 [ + - ][ + - ]: 1 : if( this->_t_created ) delete _my_geomTopoTool;
2885 : :
2886 : 1 : _t_created = true; // this one is for sure created here, even if the original gtt was not
2887 : :
2888 : : // if we were using smart pointers, we would decrease the reference to the _my_geomTopoTool, at
2889 : : // least
2890 : 1 : _my_geomTopoTool = new_gtt;
2891 : :
2892 : : // replace the _my_gsets with the new ones, from the new set
2893 [ + - ]: 1 : _my_geomTopoTool->find_geomsets( _my_gsets );
2894 : :
2895 : : // we have duplicated now the model, and the quads are part of the new _my_gsets[2]
2896 : : // we will split them now, by creating triangles along the smallest diagonal
2897 : : // maybe we will come up with a better criteria, but for the time being, it should be fine.
2898 : : // what we will do: we will remove all quads from the surface sets, and split them
2899 : :
2900 [ + - ][ + - ]: 3 : for( Range::iterator it2 = _my_gsets[2].begin(); it2 != _my_gsets[2].end(); ++it2 )
[ + - ][ + - ]
[ + + ]
2901 : : {
2902 [ + - ]: 2 : EntityHandle surface = *it2;
2903 [ + - ]: 2 : Range quads;
2904 [ + - ]: 2 : rval = _mbImpl->get_entities_by_type( surface, MBQUAD, quads );
2905 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't get quads from the surface set" );
[ # # ]
2906 [ + - ]: 2 : rval = _mbImpl->remove_entities( surface, quads );
2907 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove quads from the surface set" ); // they are not deleted, just
[ # # ]
2908 : : // removed from the set
2909 [ + - ][ + - ]: 129 : for( Range::iterator it = quads.begin(); it != quads.end(); ++it )
[ + - ][ + - ]
[ + + ][ + - ]
2910 : : {
2911 [ + - ]: 127 : EntityHandle quad = *it;
2912 : : int nnodes;
2913 : : const EntityHandle* conn;
2914 [ + - ]: 127 : rval = _mbImpl->get_connectivity( quad, conn, nnodes );
2915 [ - + ][ # # ]: 127 : MBERRORR( rval, "can't get quad connectivity" );
[ # # ]
2916 : : // get all vertices position, to see which one is the shortest diagonal
2917 [ + - ][ + + ]: 635 : CartVect pos[4];
2918 [ + - ]: 127 : rval = _mbImpl->get_coords( conn, 4, (double*)&pos[0] );
2919 [ - + ][ # # ]: 127 : MBERRORR( rval, "can't get node coordinates" );
[ # # ]
2920 [ + - ][ + - ]: 127 : bool diag1 = ( ( pos[2] - pos[0] ).length_squared() < ( pos[3] - pos[1] ).length_squared() );
[ + - ][ + - ]
2921 : : EntityHandle newTris[2];
2922 : 127 : EntityHandle tri1[3] = { conn[0], conn[1], conn[2] };
2923 : 127 : EntityHandle tri2[3] = { conn[0], conn[2], conn[3] };
2924 [ + + ]: 127 : if( !diag1 )
2925 : : {
2926 : 68 : tri1[2] = conn[3];
2927 : 68 : tri2[0] = conn[1];
2928 : : }
2929 [ + - ]: 127 : rval = _mbImpl->create_element( MBTRI, tri1, 3, newTris[0] );
2930 [ - + ][ # # ]: 127 : MBERRORR( rval, "can't create triangle 1" );
[ # # ]
2931 [ + - ]: 127 : rval = _mbImpl->create_element( MBTRI, tri2, 3, newTris[1] );
2932 [ - + ][ # # ]: 127 : MBERRORR( rval, "can't create triangle 2" );
[ # # ]
2933 [ + - ]: 127 : rval = _mbImpl->add_entities( surface, newTris, 2 );
2934 [ - + ][ # # ]: 127 : MBERRORR( rval, "can't add triangles to the set" );
[ # # ]
2935 : : }
2936 : : //
2937 : 2 : }
2938 : 10 : return MB_SUCCESS;
2939 : : }
2940 : 1 : ErrorCode FBEngine::boundary_mesh_edges_on_face( EntityHandle face, Range& boundary_mesh_edges )
2941 : : {
2942 : : // this list is used only for finding the intersecting mesh edge for starting the
2943 : : // polygonal cutting line at boundary (if !closed)
2944 [ + - ]: 1 : Range bound_edges;
2945 [ + - ]: 1 : ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
2946 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get boundary edges" );
[ # # ]
2947 [ + - ][ + - ]: 2 : for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
[ + - ][ + - ]
[ + + ]
2948 : : {
2949 [ + - ]: 1 : EntityHandle b_edge = *it;
2950 : : // get all edges in range
2951 : : // Range mesh_edges;
2952 [ + - ]: 1 : rval = _mbImpl->get_entities_by_dimension( b_edge, 1, boundary_mesh_edges );
2953 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get mesh edges" );
[ # # ]
2954 : : }
2955 : 1 : return MB_SUCCESS;
2956 : : }
2957 : 1 : ErrorCode FBEngine::boundary_nodes_on_face( EntityHandle face, std::vector< EntityHandle >& boundary_nodes )
2958 : : {
2959 : : // even if we repeat some nodes, it is OK
2960 : : // this list is used only for finding the closest boundary node for starting the
2961 : : // polygonal cutting line at boundary (if !closed)
2962 [ + - ]: 1 : Range bound_edges;
2963 [ + - ]: 1 : ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
2964 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get boundary edges" );
[ # # ]
2965 [ + - ]: 2 : Range b_nodes;
2966 [ + - ][ + - ]: 2 : for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
[ + - ][ + - ]
[ + + ]
2967 : : {
2968 [ + - ]: 1 : EntityHandle b_edge = *it;
2969 : : // get all edges in range
2970 [ + - ]: 1 : Range mesh_edges;
2971 [ + - ]: 1 : rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
2972 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get mesh edges" );
[ # # ]
2973 [ + - ]: 1 : rval = _mbImpl->get_connectivity( mesh_edges, b_nodes );
2974 [ - + ][ # # ]: 1 : MBERRORR( rval, " can't get nodes from mesh edges" );
[ # # ][ + - ]
2975 : 1 : }
2976 : : // create now a vector based on Range of nodes
2977 [ + - ][ + - ]: 1 : boundary_nodes.resize( b_nodes.size() );
2978 [ + - ][ + - ]: 1 : std::copy( b_nodes.begin(), b_nodes.end(), boundary_nodes.begin() );
[ + - ]
2979 : 2 : return MB_SUCCESS;
2980 : : }
2981 : : // used for splitting an edge
2982 : 4 : ErrorCode FBEngine::split_internal_edge( EntityHandle& edge, EntityHandle& newVertex )
2983 : : {
2984 : : // split the edge, and form 4 new triangles and 2 edges
2985 : : // get 2 triangles adjacent to edge
2986 [ + - ]: 4 : Range adj_tri;
2987 [ + - ]: 4 : ErrorCode rval = _mbImpl->get_adjacencies( &edge, 1, 2, false, adj_tri );
2988 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get adj_tris" );
[ # # ]
2989 [ + - ][ + - ]: 4 : adj_tri = subtract( adj_tri, _piercedTriangles );
2990 [ + - ][ - + ]: 4 : if( adj_tri.size() >= 3 ) { std::cout << "WARNING: non manifold geometry. Are you sure?"; }
[ # # ]
2991 [ + - ][ + - ]: 10 : for( Range::iterator it = adj_tri.begin(); it != adj_tri.end(); ++it )
[ + - ][ + - ]
[ + + ]
2992 : : {
2993 [ + - ]: 6 : EntityHandle tri = *it;
2994 [ + - ]: 6 : _piercedTriangles.insert( tri );
2995 : : const EntityHandle* conn3;
2996 : : int nnodes;
2997 [ + - ]: 6 : rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
2998 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't get nodes" );
[ # # ]
2999 : : int num1, sense, offset;
3000 [ + - ]: 6 : rval = _mbImpl->side_number( tri, edge, num1, sense, offset );
3001 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't get side number" );
[ # # ]
3002 : : // after we know the side number, we can split in 2 triangles
3003 : : // node i is opposite to edge i
3004 : 6 : int num2 = ( num1 + 1 ) % 3;
3005 : 6 : int num3 = ( num2 + 1 ) % 3;
3006 : : // the edge from num1 to num2 is split into 2 edges
3007 : 6 : EntityHandle t1[] = { conn3[num2], conn3[num3], newVertex };
3008 : 6 : EntityHandle t2[] = { conn3[num1], newVertex, conn3[num3] };
3009 : : EntityHandle newTriangle, newTriangle2;
3010 [ + - ]: 6 : rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
3011 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3012 [ + - ]: 6 : _newTriangles.insert( newTriangle );
3013 [ + - ]: 6 : rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle2 );
3014 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3015 [ + - ]: 6 : _newTriangles.insert( newTriangle2 );
3016 : : // create edges with this, indirectly
3017 [ + - ]: 6 : std::vector< EntityHandle > edges0;
3018 [ + - ]: 6 : rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
3019 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't get new edges" );
[ # # ]
3020 : 6 : edges0.clear();
3021 [ + - ]: 6 : rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
3022 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't get new edges" );
[ # # ]
3023 [ - + ]: 6 : if( debug_splits )
3024 : : {
3025 [ # # ]: 0 : std::cout << "2 (out of 4) triangles formed:\n";
3026 [ # # ]: 0 : _mbImpl->list_entity( newTriangle );
3027 [ # # ]: 0 : print_debug_triangle( newTriangle );
3028 [ # # ]: 0 : _mbImpl->list_entity( newTriangle2 );
3029 [ # # ][ + - ]: 6 : print_debug_triangle( newTriangle2 );
3030 : : }
3031 : 6 : }
3032 : 4 : return MB_SUCCESS;
3033 : : }
3034 : : // triangle split
3035 : 7 : ErrorCode FBEngine::divide_triangle( EntityHandle triangle, EntityHandle& newVertex )
3036 : : {
3037 : : //
3038 [ + - ]: 7 : _piercedTriangles.insert( triangle );
3039 : 7 : int nnodes = 0;
3040 : 7 : const EntityHandle* conn3 = NULL;
3041 [ + - ]: 7 : ErrorCode rval = _mbImpl->get_connectivity( triangle, conn3, nnodes );
3042 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't get nodes" );
[ # # ]
3043 : 7 : EntityHandle t1[] = { conn3[0], conn3[1], newVertex };
3044 : 7 : EntityHandle t2[] = { conn3[1], conn3[2], newVertex };
3045 : 7 : EntityHandle t3[] = { conn3[2], conn3[0], newVertex };
3046 : : EntityHandle newTriangle, newTriangle2, newTriangle3;
3047 [ + - ]: 7 : rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
3048 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3049 [ + - ]: 7 : _newTriangles.insert( newTriangle );
3050 [ + - ]: 7 : rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle3 );
3051 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3052 [ + - ]: 7 : _newTriangles.insert( newTriangle3 );
3053 [ + - ]: 7 : rval = _mbImpl->create_element( MBTRI, t3, 3, newTriangle2 );
3054 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3055 [ + - ]: 7 : _newTriangles.insert( newTriangle2 );
3056 : :
3057 : : // create all edges
3058 [ + - ]: 7 : std::vector< EntityHandle > edges0;
3059 [ + - ]: 7 : rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
3060 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't get new edges" );
[ # # ]
3061 : 7 : edges0.clear();
3062 [ + - ]: 7 : rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
3063 [ - + ][ # # ]: 7 : MBERRORR( rval, "can't get new edges" );
[ # # ]
3064 [ - + ]: 7 : if( debug_splits )
3065 : : {
3066 [ # # ]: 0 : std::cout << "3 triangles formed:\n";
3067 [ # # ]: 0 : _mbImpl->list_entity( newTriangle );
3068 [ # # ]: 0 : print_debug_triangle( newTriangle );
3069 [ # # ]: 0 : _mbImpl->list_entity( newTriangle3 );
3070 [ # # ]: 0 : print_debug_triangle( newTriangle3 );
3071 [ # # ]: 0 : _mbImpl->list_entity( newTriangle2 );
3072 [ # # ]: 0 : print_debug_triangle( newTriangle2 );
3073 [ # # ]: 0 : std::cout << "original nodes in tri:\n";
3074 [ # # ]: 0 : _mbImpl->list_entity( conn3[0] );
3075 [ # # ]: 0 : _mbImpl->list_entity( conn3[1] );
3076 [ # # ]: 0 : _mbImpl->list_entity( conn3[2] );
3077 : : }
3078 : 7 : return MB_SUCCESS;
3079 : : }
3080 : :
3081 : 1 : ErrorCode FBEngine::create_volume_with_direction( EntityHandle newFace1, EntityHandle newFace2, double* direction,
3082 : : EntityHandle& volume )
3083 : : {
3084 : :
3085 : : // MESHSET
3086 : : // ErrorCode rval = _mbImpl->create_meshset(MESHSET_ORDERED, new_geo_edge);
3087 [ + - ]: 1 : ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, volume );
3088 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't create volume" );
[ # # ]
3089 : :
3090 : 1 : int volumeMatId = 1; // just give a mat id, for debugging, mostly
3091 : : Tag matTag;
3092 [ + - ]: 1 : rval = _mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matTag );
3093 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get material tag" );
[ # # ]
3094 : :
3095 [ + - ]: 1 : rval = _mbImpl->tag_set_data( matTag, &volume, 1, &volumeMatId );
3096 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't set material tag value on volume" );
[ # # ]
3097 : :
3098 : : // get the edges of those 2 faces, and get the vertices of those edges
3099 : : // in order, they should be created in the same order (?); check for that, anyway
3100 [ + - ]: 1 : rval = _mbImpl->add_parent_child( volume, newFace1 );
3101 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add first face to volume" );
[ # # ]
3102 : :
3103 [ + - ]: 1 : rval = _mbImpl->add_parent_child( volume, newFace2 );
3104 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add second face to volume" );
[ # # ]
3105 : :
3106 : : // first is bottom, so it is negatively oriented
3107 [ + - ]: 1 : rval = _my_geomTopoTool->add_geo_set( volume, 3 );
3108 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add volume to the gtt" );
[ # # ]
3109 : :
3110 : : // set senses
3111 : : // bottom face is negatively oriented, its normal is toward interior of the volume
3112 [ + - ]: 1 : rval = _my_geomTopoTool->set_sense( newFace1, volume, -1 );
3113 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't set bottom face sense to the volume" );
[ # # ]
3114 : :
3115 : : // the top face is positively oriented
3116 [ + - ]: 1 : rval = _my_geomTopoTool->set_sense( newFace2, volume, 1 );
3117 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't set top face sense to the volume" );
[ # # ]
3118 : :
3119 : : // the children should be in the same direction
3120 : : // get the side edges of each face, and form lateral faces, along direction
3121 [ + - ]: 1 : std::vector< EntityHandle > edges1;
3122 [ + - ]: 2 : std::vector< EntityHandle > edges2;
3123 : :
3124 [ + - ]: 1 : rval = _mbImpl->get_child_meshsets( newFace1, edges1 ); // no hops
3125 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get children edges or first face, bottom" );
[ # # ]
3126 : :
3127 [ + - ]: 1 : rval = _mbImpl->get_child_meshsets( newFace2, edges2 ); // no hops
3128 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get children edges for second face, top" );
[ # # ]
3129 : :
3130 [ - + ][ # # ]: 1 : if( edges1.size() != edges2.size() ) MBERRORR( MB_FAILURE, "wrong correspondence " );
[ # # ]
3131 : :
3132 [ + + ]: 5 : for( unsigned int i = 0; i < edges1.size(); ++i )
3133 : : {
3134 : : EntityHandle newLatFace;
3135 [ + - ][ + - ]: 4 : rval = weave_lateral_face_from_edges( edges1[i], edges2[i], direction, newLatFace );
[ + - ]
3136 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't weave lateral face" );
[ # # ]
3137 [ + - ]: 4 : rval = _mbImpl->add_parent_child( volume, newLatFace );
3138 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add lateral face to volume" );
[ # # ]
3139 : :
3140 : : // set sense as positive
3141 [ + - ]: 4 : rval = _my_geomTopoTool->set_sense( newLatFace, volume, 1 );
3142 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set lateral face sense to the volume" );
[ # # ]
3143 : : }
3144 : :
3145 [ + - ]: 1 : rval = set_default_neumann_tags();
3146 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't set new neumann tags" );
[ # # ]
3147 : :
3148 : 2 : return MB_SUCCESS;
3149 : : }
3150 : :
3151 : 8 : ErrorCode FBEngine::get_nodes_from_edge( EntityHandle gedge, std::vector< EntityHandle >& nodes )
3152 : : {
3153 [ + - ]: 8 : std::vector< EntityHandle > ents;
3154 [ + - ]: 8 : ErrorCode rval = _mbImpl->get_entities_by_type( gedge, MBEDGE, ents );
3155 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
3156 [ - + ]: 8 : if( ents.size() < 1 ) return MB_FAILURE;
3157 : :
3158 [ + - ]: 8 : nodes.resize( ents.size() + 1 );
3159 : 8 : const EntityHandle* conn = NULL;
3160 : : int len;
3161 [ + + ]: 1414 : for( unsigned int i = 0; i < ents.size(); ++i )
3162 : : {
3163 [ + - ][ + - ]: 1406 : rval = _mbImpl->get_connectivity( ents[i], conn, len );
3164 [ - + ][ # # ]: 1406 : MBERRORR( rval, "can't get edge connectivity" );
[ # # ]
3165 [ + - ]: 1406 : nodes[i] = conn[0];
3166 : : }
3167 : : // the last one is conn[1]
3168 [ + - ]: 8 : nodes[ents.size()] = conn[1];
3169 : 8 : return MB_SUCCESS;
3170 : : }
3171 : 4 : ErrorCode FBEngine::weave_lateral_face_from_edges( EntityHandle bEdge, EntityHandle tEdge, double* direction,
3172 : : EntityHandle& newLatFace )
3173 : : {
3174 : : // in weird cases might need to create new vertices in the interior;
3175 : : // in most cases, it is OK
3176 : :
3177 [ + - ]: 4 : ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, newLatFace );
3178 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't create new lateral face" );
[ # # ]
3179 : :
3180 : : EntityHandle v[4]; // vertex sets
3181 : : // bot edge will be v1->v2
3182 : : // top edge will be v3->v4
3183 : : // we need to create edges from v1 to v3 and from v2 to v4
3184 [ + - ]: 4 : std::vector< EntityHandle > adj;
3185 [ + - ]: 4 : rval = _mbImpl->get_child_meshsets( bEdge, adj );
3186 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get children nodes" );
[ # # ]
3187 : 4 : bool periodic = false;
3188 [ - + ]: 4 : if( adj.size() == 1 )
3189 : : {
3190 [ # # ]: 0 : v[0] = v[1] = adj[0];
3191 : 0 : periodic = true;
3192 : : }
3193 : : else
3194 : : {
3195 [ + - ]: 4 : v[0] = adj[0];
3196 [ + - ]: 4 : v[1] = adj[1];
3197 : : }
3198 : : int senseB;
3199 [ + - ]: 4 : rval = getEgVtxSense( bEdge, v[0], v[1], senseB );
3200 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get bottom edge sense" );
[ # # ]
3201 [ + + ]: 4 : if( -1 == senseB )
3202 : : {
3203 [ + - ]: 1 : v[1] = adj[0]; // so , bEdge will be oriented from v[0] to v[1], and will start at
3204 : : // nodes1, coords1..
3205 [ + - ]: 1 : v[0] = adj[1];
3206 : : }
3207 : 4 : adj.clear();
3208 [ + - ]: 4 : rval = _mbImpl->get_child_meshsets( tEdge, adj );
3209 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get children nodes" );
[ # # ]
3210 [ - + ]: 4 : if( adj.size() == 1 )
3211 : : {
3212 [ # # ]: 0 : v[2] = v[3] = adj[0];
3213 [ # # ][ # # ]: 0 : if( !periodic ) MBERRORR( MB_FAILURE, "top edge is periodic, but bottom edge is not" );
[ # # ]
3214 : : }
3215 : : else
3216 : : {
3217 [ + - ]: 4 : v[2] = adj[0];
3218 [ + - ]: 4 : v[3] = adj[1];
3219 [ - + ][ # # ]: 4 : if( periodic ) MBERRORR( MB_FAILURE, "top edge is not periodic, but bottom edge is" );
[ # # ]
3220 : : }
3221 : :
3222 : : // now, using nodes on bottom edge and top edge, create triangles, oriented outwards the
3223 : : // volume (sense positive on bottom edge)
3224 [ + - ]: 8 : std::vector< EntityHandle > nodes1;
3225 [ + - ]: 4 : rval = get_nodes_from_edge( bEdge, nodes1 );
3226 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get nodes from bott edge" );
[ # # ]
3227 : :
3228 [ + - ]: 8 : std::vector< EntityHandle > nodes2;
3229 [ + - ]: 4 : rval = get_nodes_from_edge( tEdge, nodes2 );
3230 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get nodes from top edge" );
[ # # ]
3231 : :
3232 [ + - ][ + - ]: 8 : std::vector< CartVect > coords1, coords2;
3233 [ + - ]: 4 : coords1.resize( nodes1.size() );
3234 [ + - ]: 4 : coords2.resize( nodes2.size() );
3235 : :
3236 : 4 : int N1 = (int)nodes1.size();
3237 : 4 : int N2 = (int)nodes2.size();
3238 : :
3239 [ + - ][ + - ]: 4 : rval = _mbImpl->get_coords( &( nodes1[0] ), nodes1.size(), (double*)&( coords1[0] ) );
[ + - ]
3240 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get coords of nodes from bott edge" );
[ # # ]
3241 : :
3242 [ + - ][ + - ]: 4 : rval = _mbImpl->get_coords( &( nodes2[0] ), nodes2.size(), (double*)&( coords2[0] ) );
[ + - ]
3243 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get coords of nodes from top edge" );
[ # # ]
3244 [ + - ]: 4 : CartVect up( direction );
3245 : :
3246 : : // see if the start and end coordinates are matching, if not, reverse edge 2 nodes and
3247 : : // coordinates
3248 [ + - ][ + - ]: 4 : CartVect v1 = ( coords1[0] - coords2[0] ) * up;
[ + - ][ + - ]
3249 [ + - ][ + - ]: 4 : CartVect v2 = ( coords1[0] - coords2[N2 - 1] ) * up;
[ + - ][ + - ]
3250 [ + - ][ + - ]: 4 : if( v1.length_squared() > v2.length_squared() )
[ - + ]
3251 : : {
3252 : : // we need to reverse coords2 and node2, as nodes are not above each other
3253 : : // the edges need to be found between v0 and v3, v1 and v2!
3254 [ # # ]: 0 : for( unsigned int k = 0; k < nodes2.size() / 2; k++ )
3255 : : {
3256 [ # # ]: 0 : EntityHandle tmp = nodes2[k];
3257 [ # # ][ # # ]: 0 : nodes2[k] = nodes2[N2 - 1 - k];
3258 [ # # ]: 0 : nodes2[N2 - 1 - k] = tmp;
3259 [ # # ]: 0 : CartVect tv = coords2[k];
3260 [ # # ][ # # ]: 0 : coords2[k] = coords2[N2 - 1 - k];
3261 [ # # ]: 0 : coords2[N2 - 1 - k] = tv;
3262 : : }
3263 : : }
3264 : : // make sure v[2] has nodes2[0], if not, reverse v[2] and v[3]
3265 [ + - ][ + - ]: 4 : if( !_mbImpl->contains_entities( v[2], &( nodes2[0] ), 1 ) )
[ + + ]
3266 : : {
3267 : : // reverse v[2] and v[3], so v[2] will be above v[0]
3268 : 1 : EntityHandle tmp = v[2];
3269 : 1 : v[2] = v[3];
3270 : 1 : v[3] = tmp;
3271 : : }
3272 : : // now we know that v[0]--v[3] will be vertex sets in the order we want
3273 [ + - ][ + - ]: 4 : EntityHandle nd[4] = { nodes1[0], nodes1[N1 - 1], nodes2[0], nodes2[N2 - 1] };
[ + - ][ + - ]
3274 : :
3275 : 4 : adj.clear();
3276 : : EntityHandle e1, e2;
3277 : : // find edge 1 between v[0] and v[2], and e2 between v[1] and v[3]
3278 [ + - ]: 4 : rval = _mbImpl->get_parent_meshsets( v[0], adj );
3279 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get edges connected to vertex set 1" );
[ # # ]
3280 : 4 : bool found = false;
3281 [ + + ]: 15 : for( unsigned int j = 0; j < adj.size(); j++ )
3282 : : {
3283 [ + - ]: 11 : EntityHandle ed = adj[j];
3284 [ + - ]: 11 : Range vertices;
3285 [ + - ]: 11 : rval = _mbImpl->get_child_meshsets( ed, vertices );
3286 [ + - ][ + - ]: 11 : if( vertices.find( v[2] ) != vertices.end() )
[ + - ][ + + ]
3287 : : {
3288 : 3 : found = true;
3289 : 3 : e1 = ed;
3290 [ + + ]: 11 : break;
3291 : : }
3292 : 8 : }
3293 [ + + ]: 4 : if( !found )
3294 : : {
3295 : : // create an edge from v[0] to v[2]
3296 [ + - ]: 1 : rval = _mbImpl->create_meshset( MESHSET_SET, e1 );
3297 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't create edge 1" );
[ # # ]
3298 : :
3299 [ + - ]: 1 : rval = _mbImpl->add_parent_child( e1, v[0] );
3300 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3301 : :
3302 [ + - ]: 1 : rval = _mbImpl->add_parent_child( e1, v[2] );
3303 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3304 : :
3305 : 1 : EntityHandle nn2[2] = { nd[0], nd[2] };
3306 : : EntityHandle medge;
3307 [ + - ]: 1 : rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
3308 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't create mesh edge" );
[ # # ]
3309 : :
3310 [ + - ]: 1 : rval = _mbImpl->add_entities( e1, &medge, 1 );
3311 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add mesh edge to geo edge" );
[ # # ]
3312 : :
3313 [ + - ]: 1 : rval = this->_my_geomTopoTool->add_geo_set( e1, 1 );
3314 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't add edge to gtt" );
[ # # ]
3315 : : }
3316 : :
3317 : : // find the edge from v2 to v4 (if not, create one)
3318 [ + - ]: 4 : rval = _mbImpl->get_parent_meshsets( v[1], adj );
3319 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get edges connected to vertex set 2" );
[ # # ]
3320 : 4 : found = false;
3321 [ + + ]: 20 : for( unsigned int i = 0; i < adj.size(); i++ )
3322 : : {
3323 [ + - ]: 16 : EntityHandle ed = adj[i];
3324 [ + - ]: 16 : Range vertices;
3325 [ + - ]: 16 : rval = _mbImpl->get_child_meshsets( ed, vertices );
3326 [ + - ][ + - ]: 16 : if( vertices.find( v[3] ) != vertices.end() )
[ + - ][ + + ]
3327 : : {
3328 : 1 : found = true;
3329 : 1 : e2 = ed;
3330 [ + + ]: 16 : break;
3331 : : }
3332 : 15 : }
3333 [ + + ]: 4 : if( !found )
3334 : : {
3335 : : // create an edge from v2 to v4
3336 [ + - ]: 3 : rval = _mbImpl->create_meshset( MESHSET_SET, e2 );
3337 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't create edge 1" );
[ # # ]
3338 : :
3339 [ + - ]: 3 : rval = _mbImpl->add_parent_child( e2, v[1] );
3340 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3341 : :
3342 [ + - ]: 3 : rval = _mbImpl->add_parent_child( e2, v[3] );
3343 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3344 : :
3345 : 3 : EntityHandle nn2[2] = { nd[1], nd[3] };
3346 : : EntityHandle medge;
3347 [ + - ]: 3 : rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
3348 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't create mesh edge" );
[ # # ]
3349 : :
3350 [ + - ]: 3 : rval = _mbImpl->add_entities( e2, &medge, 1 );
3351 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't add mesh edge to geo edge" );
[ # # ]
3352 : :
3353 [ + - ]: 3 : rval = _my_geomTopoTool->add_geo_set( e2, 1 );
3354 [ - + ][ # # ]: 3 : MBERRORR( rval, "can't add edge to gtt" );
[ # # ]
3355 : : }
3356 : :
3357 : : // now we have the four edges, add them to the face, as children
3358 : :
3359 : : // add children to face
3360 [ + - ]: 4 : rval = _mbImpl->add_parent_child( newLatFace, bEdge );
3361 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3362 : :
3363 [ + - ]: 4 : rval = _mbImpl->add_parent_child( newLatFace, tEdge );
3364 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3365 : :
3366 [ + - ]: 4 : rval = _mbImpl->add_parent_child( newLatFace, e1 );
3367 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3368 : :
3369 [ + - ]: 4 : rval = _mbImpl->add_parent_child( newLatFace, e2 );
3370 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add parent - child relation" );
[ # # ]
3371 : :
3372 [ + - ]: 4 : rval = _my_geomTopoTool->add_geo_set( newLatFace, 2 );
3373 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't add face to gtt" );
[ # # ]
3374 : : // add senses
3375 : : //
3376 [ + - ]: 4 : rval = _my_geomTopoTool->set_sense( bEdge, newLatFace, 1 );
3377 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set bottom edge sense to the lateral face" );
[ # # ]
3378 : :
3379 : : int Tsense;
3380 [ + - ]: 4 : rval = getEgVtxSense( tEdge, v[3], v[2], Tsense );
3381 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't get top edge sense" );
[ # # ]
3382 : : // we need to see what sense has topEdge in face
3383 [ + - ]: 4 : rval = _my_geomTopoTool->set_sense( tEdge, newLatFace, Tsense );
3384 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set top edge sense to the lateral face" );
[ # # ]
3385 : :
3386 [ + - ]: 4 : rval = _my_geomTopoTool->set_sense( e1, newLatFace, -1 );
3387 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set first vert edge sense" );
[ # # ]
3388 : :
3389 [ + - ]: 4 : rval = _my_geomTopoTool->set_sense( e2, newLatFace, 1 );
3390 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't set second edge sense to the lateral face" );
[ # # ]
3391 : : // first, create edges along direction, for the
3392 : :
3393 : 4 : int indexB = 0, indexT = 0; // indices of the current nodes in the weaving process
3394 : : // weaving is either up or down; the triangles are oriented positively either way
3395 : : // up is nodes1[indexB], nodes2[indexT+1], nodes2[indexT]
3396 : : // down is nodes1[indexB], nodes1[indexB+1], nodes2[indexT]
3397 : : // the decision to weave up or down is based on which one is closer to the direction normal
3398 : : /*
3399 : : *
3400 : : * --------*------*-----------* ^
3401 : : * / . \ . . ------> dir1 | up
3402 : : * /. \ . . |
3403 : : * -----*------------*----*
3404 : : *
3405 : : */
3406 : : // we have to change the logic to account for cases when the curve in xy plane is not straight
3407 : : // (for example, when merging curves with a min_dot < -0.5, which means that the edges
3408 : : // could be even closed (periodic), with one vertex
3409 : : // the top and bottom curves should follow the same path in the "xy" plane (better to say
3410 : : // the plane perpendicular to the direction of weaving)
3411 : : // in this logic, the vector dir1 varies along the curve !!!
3412 [ + - ][ + - ]: 4 : CartVect dir1 = coords1[1] - coords1[0]; // we should have at least 2 nodes, N1>=2!!
[ + - ]
3413 : :
3414 [ + - ]: 4 : CartVect planeNormal = dir1 * up;
3415 [ + - ]: 4 : dir1 = up * planeNormal;
3416 [ + - ]: 4 : dir1.normalize();
3417 : : // this direction will be updated constantly with the direction of last edge added
3418 : 4 : bool weaveDown = true;
3419 : :
3420 [ + - ]: 4 : CartVect startP = coords1[0]; // used for origin of comparisons
3421 : : while( 1 )
3422 : : {
3423 [ + + ][ + + ]: 1410 : if( ( indexB == N1 - 1 ) && ( indexT == N2 - 1 ) ) break; // we cannot advance anymore
3424 [ + + ]: 1406 : if( indexB == N1 - 1 ) { weaveDown = false; }
3425 [ + + ]: 1403 : else if( indexT == N2 - 1 )
3426 : : {
3427 : 1 : weaveDown = true;
3428 : : }
3429 : : else
3430 : : {
3431 : : // none are at the end, we have to decide which way to go, based on which index + 1 is
3432 : : // closer
3433 [ + - ][ + - ]: 1402 : double proj1 = ( coords1[indexB + 1] - startP ) % dir1;
[ + - ]
3434 [ + - ][ + - ]: 1402 : double proj2 = ( coords2[indexT + 1] - startP ) % dir1;
[ + - ]
3435 [ + + ]: 1402 : if( proj1 < proj2 )
3436 : 748 : weaveDown = true;
3437 : : else
3438 : 654 : weaveDown = false;
3439 : : }
3440 [ + - ][ + - ]: 1406 : EntityHandle nTri[3] = { nodes1[indexB], 0, nodes2[indexT] };
3441 [ + + ]: 1406 : if( weaveDown )
3442 : : {
3443 [ + - ]: 749 : nTri[1] = nodes1[indexB + 1];
3444 [ + - ]: 749 : nTri[2] = nodes2[indexT];
3445 : 749 : indexB++;
3446 : : }
3447 : : else
3448 : : {
3449 [ + - ]: 657 : nTri[1] = nodes2[indexT + 1];
3450 : 657 : indexT++;
3451 : : }
3452 : : EntityHandle triangle;
3453 [ + - ]: 1406 : rval = _mbImpl->create_element( MBTRI, nTri, 3, triangle );
3454 [ - + ][ # # ]: 1406 : MBERRORR( rval, "can't create triangle" );
[ # # ]
3455 : :
3456 [ + - ]: 1406 : rval = _mbImpl->add_entities( newLatFace, &triangle, 1 );
3457 [ - + ][ # # ]: 1406 : MBERRORR( rval, "can't add triangle to face set" );
[ # # ]
3458 [ + + ]: 1406 : if( weaveDown )
3459 : : {
3460 : : // increase was from nodes1[indexB-1] to nodes1[indexb]
3461 [ + - ][ + - ]: 749 : dir1 = coords1[indexB] - coords1[indexB - 1]; // we should have at least 2 nodes, N1>=2!!
[ + - ]
3462 : : }
3463 : : else
3464 : : {
3465 [ + - ][ + - ]: 657 : dir1 = coords2[indexT] - coords2[indexT - 1];
[ + - ]
3466 : : }
3467 [ + - ][ + - ]: 1406 : dir1 = up * ( dir1 * up );
3468 [ + - ]: 1406 : dir1.normalize();
3469 : : }
3470 : : // we do not check yet if the triangles are inverted
3471 : : // if yes, we should correct them. HOW?
3472 : : // we probably need a better meshing strategy, one that can overcome really bad meshes.
3473 : : // again, these faces are not what we should use for geometry, maybe we should use the
3474 : : // extruded quads, identified AFTER hexa are created.
3475 : : // right now, I see only a cosmetic use of these lateral faces
3476 : : // the whole idea of volume maybe is overrated
3477 : : // volume should be just quads extruded from bottom ?
3478 : : //
3479 : 8 : return MB_SUCCESS;
3480 : : }
3481 : : // this will be used during creation of the volume, to set unique
3482 : : // tags on all surfaces
3483 : : // it is changing original tags, so do not use it if you want to preserve
3484 : : // initial neumann tags
3485 : 1 : ErrorCode FBEngine::set_default_neumann_tags()
3486 : : {
3487 : : // these are for debugging purposes only
3488 : : // check the initial tag, then
3489 : : Tag ntag;
3490 [ + - ]: 1 : ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag );
3491 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get tag handle" );
[ # # ]
3492 : : // get all surfaces in the model now
3493 [ + - ][ + + ]: 11 : Range sets[5];
[ # # ]
3494 [ + - ]: 1 : rval = _my_geomTopoTool->find_geomsets( sets );
3495 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't get geo sets" );
[ # # ]
3496 [ + - ]: 1 : int nfaces = (int)sets[2].size();
3497 [ + - ][ + - ]: 1 : int* vals = new int[nfaces];
3498 [ + + ]: 9 : for( int i = 0; i < nfaces; i++ )
3499 : : {
3500 : 8 : vals[i] = i + 1;
3501 : : }
3502 [ + - ]: 1 : rval = _mbImpl->tag_set_data( ntag, sets[2], (void*)vals );
3503 [ - + ][ # # ]: 1 : MBERRORR( rval, "can't set tag values for neumann sets" );
[ # # ]
3504 : :
3505 [ + - ]: 1 : delete[] vals;
3506 : :
3507 [ + + ]: 6 : return MB_SUCCESS;
[ # # # # ]
3508 : : }
3509 : : // a reverse operation for splitting an gedge at a mesh node
3510 : 4 : ErrorCode FBEngine::chain_edges( double min_dot )
3511 : : {
3512 [ + - ][ + + ]: 44 : Range sets[5];
[ # # ]
3513 : : ErrorCode rval;
3514 : : while( 1 ) // break out only if no edges are chained
3515 : : {
3516 [ + - ]: 6 : rval = _my_geomTopoTool->find_geomsets( sets );
3517 [ - + ][ # # ]: 6 : MBERRORR( rval, "can't get geo sets" );
[ # # ]
3518 : : // these gentities are "always" current, while the ones in this-> _my_gsets[0:4] are
3519 : : // the "originals" before FBEngine modifications
3520 [ + - ]: 6 : int nedges = (int)sets[1].size();
3521 : : // as long as we have chainable edges, continue;
3522 : 6 : bool chain = false;
3523 [ + + ]: 37 : for( int i = 0; i < nedges; i++ )
3524 : : {
3525 [ + - ]: 33 : EntityHandle edge = sets[1][i];
3526 : : EntityHandle next_edge;
3527 : 33 : bool chainable = false;
3528 [ + - ]: 33 : rval = chain_able_edge( edge, min_dot, next_edge, chainable );
3529 [ - + ][ # # ]: 33 : MBERRORR( rval, "can't determine chain-ability" );
[ # # ]
3530 [ + + ]: 33 : if( chainable )
3531 : : {
3532 [ + - ]: 2 : rval = chain_two_edges( edge, next_edge );
3533 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't chain 2 edges" );
[ # # ]
3534 : 2 : chain = true;
3535 : 2 : break; // interrupt for loop
3536 : : }
3537 : : }
3538 [ + + ]: 6 : if( !chain )
3539 : : {
3540 : 6 : break; // break out of while loop
3541 : : }
3542 : : }
3543 [ + + ]: 24 : return MB_SUCCESS;
[ # # # # ]
3544 : : }
3545 : :
3546 : : // determine if from the end of edge we can extend with another edge; return also the
3547 : : // extension edge (next_edge)
3548 : 33 : ErrorCode FBEngine::chain_able_edge( EntityHandle edge, double min_dot, EntityHandle& next_edge, bool& chainable )
3549 : : {
3550 : : // get the end, then get all parents of end
3551 : : // see if some are the starts of
3552 : 33 : chainable = false;
3553 : : EntityHandle v1, v2;
3554 [ + - ]: 33 : ErrorCode rval = get_vert_edges( edge, v1, v2 );
3555 [ - + ][ # # ]: 33 : MBERRORR( rval, "can't get vertices" );
[ # # ]
3556 [ + + ]: 33 : if( v1 == v2 ) return MB_SUCCESS; // it is periodic, can't chain it with another edge!
3557 : :
3558 : : // v2 is a set, get its parents, which should be edges
3559 [ + - ]: 27 : Range edges;
3560 [ + - ]: 27 : rval = _mbImpl->get_parent_meshsets( v2, edges );
3561 [ - + ][ # # ]: 27 : MBERRORR( rval, "can't get parents of vertex set" );
[ # # ]
3562 : : // get parents of current edge (faces)
3563 [ + - ]: 54 : Range faces;
3564 [ + - ]: 27 : rval = _mbImpl->get_parent_meshsets( edge, faces );
3565 [ - + ][ # # ]: 27 : MBERRORR( rval, "can't get parents of edge set" );
[ # # ]
3566 : : // get also the last edge "tangent" at the vertex
3567 [ + - ]: 54 : std::vector< EntityHandle > mesh_edges;
3568 [ + - ]: 27 : rval = _mbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges );
3569 [ - + ][ # # ]: 27 : MBERRORR( rval, "can't get mesh edges from edge set" );
[ # # ]
3570 [ + - ]: 27 : EntityHandle lastMeshEdge = mesh_edges[mesh_edges.size() - 1];
3571 : 27 : const EntityHandle* conn2 = NULL;
3572 : 27 : int len = 0;
3573 [ + - ]: 27 : rval = _mbImpl->get_connectivity( lastMeshEdge, conn2, len );
3574 [ - + ][ # # ]: 27 : MBERRORR( rval, "can't connectivity of last mesh edge" );
[ # # ]
3575 : : // get the coordinates of last edge
3576 [ - + ][ # # ]: 27 : if( len != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
[ # # ]
3577 [ + - ][ + + ]: 81 : CartVect P[2];
3578 [ + - ]: 27 : rval = _mbImpl->get_coords( conn2, len, (double*)&P[0] );
3579 [ - + ][ # # ]: 27 : MBERRORR( rval, "Failed to get coordinates" );
[ # # ]
3580 : :
3581 [ + - ]: 27 : CartVect vec1( P[1] - P[0] );
3582 [ + - ]: 27 : vec1.normalize();
3583 [ + - ][ + - ]: 85 : for( Range::iterator edgeIter = edges.begin(); edgeIter != edges.end(); ++edgeIter )
[ + - ][ + - ]
[ + + ]
3584 : : {
3585 [ + - ]: 58 : EntityHandle otherEdge = *edgeIter;
3586 [ + + ]: 58 : if( edge == otherEdge ) continue;
3587 : : // get faces parents of this edge
3588 [ + - ]: 31 : Range faces2;
3589 [ + - ]: 31 : rval = _mbImpl->get_parent_meshsets( otherEdge, faces2 );
3590 [ - + ][ # # ]: 31 : MBERRORR( rval, "can't get parents of other edge set" );
[ # # ]
3591 [ + - ][ + + ]: 31 : if( faces != faces2 ) continue;
3592 : : // now, if the first mesh edge is within given angle, we can go on
3593 [ + - ]: 23 : std::vector< EntityHandle > mesh_edges2;
3594 [ + - ]: 23 : rval = _mbImpl->get_entities_by_type( otherEdge, MBEDGE, mesh_edges2 );
3595 [ - + ][ # # ]: 23 : MBERRORR( rval, "can't get mesh edges from other edge set" );
[ # # ]
3596 [ + - ]: 23 : EntityHandle firstMeshEdge = mesh_edges2[0];
3597 : : const EntityHandle* conn22;
3598 : : int len2;
3599 [ + - ]: 23 : rval = _mbImpl->get_connectivity( firstMeshEdge, conn22, len2 );
3600 [ - + ][ # # ]: 23 : MBERRORR( rval, "can't connectivity of first mesh edge" );
[ # # ]
3601 [ - + ][ # # ]: 23 : if( len2 != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
[ # # ]
3602 [ - + ]: 23 : if( conn2[1] != conn22[0] ) continue; // the mesh edges are not one after the other
3603 : : // get the coordinates of first edge
3604 : :
3605 : : // CartVect P2[2];
3606 [ + - ]: 23 : rval = _mbImpl->get_coords( conn22, len, (double*)&P[0] );
3607 [ + - ]: 23 : CartVect vec2( P[1] - P[0] );
3608 [ + - ]: 23 : vec2.normalize();
3609 [ + - ][ + + ]: 23 : if( vec1 % vec2 < min_dot ) continue;
3610 : : // we found our edge, victory! we can get out
3611 : 2 : next_edge = otherEdge;
3612 : 2 : chainable = true;
3613 [ + + ][ + + ]: 31 : return MB_SUCCESS;
3614 : 0 : }
3615 : :
3616 : 58 : return MB_SUCCESS; // in general, hard to come by chain-able edges
3617 : : }
3618 : 2 : ErrorCode FBEngine::chain_two_edges( EntityHandle edge, EntityHandle next_edge )
3619 : : {
3620 : : // the biggest thing is to see the sense tags; or maybe not...
3621 : : // they should be correct :)
3622 : : // get the vertex sets
3623 : : EntityHandle v11, v12, v21, v22;
3624 [ + - ]: 2 : ErrorCode rval = get_vert_edges( edge, v11, v12 );
3625 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't get vert sets" );
[ # # ]
3626 [ + - ]: 2 : rval = get_vert_edges( next_edge, v21, v22 );
3627 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't get vert sets" );
[ # # ]
3628 [ - + ]: 2 : assert( v12 == v21 );
3629 [ + - ]: 2 : std::vector< EntityHandle > mesh_edges;
3630 [ + - ]: 2 : rval = MBI->get_entities_by_type( next_edge, MBEDGE, mesh_edges );
3631 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't get mesh edges" );
[ # # ]
3632 : :
3633 [ + - ][ + - ]: 2 : rval = _mbImpl->add_entities( edge, &mesh_edges[0], (int)mesh_edges.size() );
3634 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't add new mesh edges" );
[ # # ]
3635 : : // remove the child - parent relation for second vertex of first edge
3636 [ + - ]: 2 : rval = _mbImpl->remove_parent_child( edge, v12 );
3637 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove parent - child relation between first edge and middle vertex" );
[ # # ]
3638 : :
3639 [ + - ]: 2 : if( v22 != v11 ) // the edge would become periodic, do not add again the relationship
3640 : : {
3641 [ + - ]: 2 : rval = _mbImpl->add_parent_child( edge, v22 );
3642 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't add second vertex to edge " );
[ # # ]
3643 : : }
3644 : : // we can now safely eliminate next_edge
3645 [ + - ]: 2 : rval = _mbImpl->remove_parent_child( next_edge, v21 );
3646 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove child - parent relation " );
[ # # ]
3647 : :
3648 [ + - ]: 2 : rval = _mbImpl->remove_parent_child( next_edge, v22 );
3649 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove child - parent relation " );
[ # # ]
3650 : :
3651 : : // remove the next_edge relation to the faces
3652 [ + - ]: 4 : Range faces;
3653 [ + - ]: 2 : rval = _mbImpl->get_parent_meshsets( next_edge, faces );
3654 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't get parent faces " );
[ # # ]
3655 : :
3656 [ + - ][ + - ]: 6 : for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
[ + - ][ + - ]
[ + + ]
3657 : : {
3658 [ + - ]: 4 : EntityHandle ff = *it;
3659 [ + - ]: 4 : rval = _mbImpl->remove_parent_child( ff, next_edge );
3660 [ - + ][ # # ]: 4 : MBERRORR( rval, "can't remove parent-edge rel " );
[ # # ]
3661 : : }
3662 : :
3663 [ + - ]: 2 : rval = _mbImpl->delete_entities( &next_edge, 1 );
3664 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove edge set " );
[ # # ]
3665 : :
3666 : : // delete the vertex set that is idle now (v12 = v21)
3667 [ + - ]: 2 : rval = _mbImpl->delete_entities( &v12, 1 );
3668 [ - + ][ # # ]: 2 : MBERRORR( rval, "can't remove edge set " );
[ # # ]
3669 : 4 : return MB_SUCCESS;
3670 : : }
3671 : 37 : ErrorCode FBEngine::get_vert_edges( EntityHandle edge, EntityHandle& v1, EntityHandle& v2 )
3672 : : {
3673 : : // need to decide first or second vertex
3674 : : // important for moab
3675 : :
3676 [ + - ]: 37 : Range children;
3677 : : // EntityHandle v1, v2;
3678 [ + - ]: 37 : ErrorCode rval = _mbImpl->get_child_meshsets( edge, children );
3679 [ - + ][ # # ]: 37 : MBERRORR( rval, "can't get child meshsets" );
[ # # ]
3680 [ + - ][ + + ]: 37 : if( children.size() == 1 )
3681 : : {
3682 : : // this is periodic edge, get out early
3683 [ + - ]: 6 : v1 = children[0];
3684 : 6 : v2 = v1;
3685 : 6 : return MB_SUCCESS;
3686 : : }
3687 [ + - ][ - + ]: 31 : else if( children.size() > 2 )
3688 [ # # ][ # # ]: 0 : MBERRORR( MB_FAILURE, "too many vertices in one edge" );
3689 : : // edge: get one vertex as part of the vertex set
3690 [ + - ]: 62 : Range entities;
3691 [ + - ][ + - ]: 31 : rval = MBI->get_entities_by_type( children[0], MBVERTEX, entities );
3692 [ - + ][ # # ]: 31 : MBERRORR( rval, "can't get entities from vertex set" );
[ # # ]
3693 [ + - ][ - + ]: 31 : if( entities.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh nodes in vertex set" );
[ # # ][ # # ]
3694 [ + - ]: 31 : EntityHandle node0 = entities[0]; // the first vertex
3695 [ + - ]: 31 : entities.clear();
3696 : :
3697 : : // now get the edges, and get the first node and the last node in sequence of edges
3698 : : // the order is important...
3699 : : // these are ordered sets !!
3700 [ + - ]: 62 : std::vector< EntityHandle > ents;
3701 [ + - ]: 31 : rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
3702 [ - + ][ # # ]: 31 : MBERRORR( rval, "can't get mesh edges" );
[ # # ]
3703 [ - + ][ # # ]: 31 : if( ents.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh edges in edge set" );
[ # # ]
3704 : :
3705 : 31 : const EntityHandle* conn = NULL;
3706 : : int len;
3707 [ + - ][ + - ]: 31 : rval = MBI->get_connectivity( ents[0], conn, len );
3708 [ - + ][ # # ]: 31 : MBERRORR( rval, "can't connectivity of first mesh edge" );
[ # # ]
3709 : :
3710 [ + + ]: 31 : if( conn[0] == node0 )
3711 : : {
3712 [ + - ]: 21 : v1 = children[0];
3713 [ + - ]: 21 : v2 = children[1];
3714 : : }
3715 : : else // the other way around, although we should check (we are paranoid)
3716 : : {
3717 [ + - ]: 10 : v2 = children[0];
3718 [ + - ]: 10 : v1 = children[1];
3719 : : }
3720 : :
3721 : 68 : return MB_SUCCESS;
3722 : : }
3723 [ + - ][ + - ]: 44 : } // namespace moab
|
