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