Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
FBEngine.cpp
Go to the documentation of this file.
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"
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
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 
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  {
199  _t_created = true;
200  }
201  // should this be part of the constructor or not?
202  // Init();
203 }
205 {
206  clean();
207  _smooth = false;
208 }
209 
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 
238 {
239  if( !_initialized )
240  {
241  if( !_my_geomTopoTool ) return MB_FAILURE;
242 
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 
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 }
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
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 
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 
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 
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 =
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 
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 }
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 }
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 
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 
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 
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 
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 
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 
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 
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 }
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 
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 
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  */
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!)
1610 
1611  MBERRORR( rval, "can't delete triangles" );
1613  // delete edges that are broke up in 2
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)
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( &currentTriangle, 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 
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
1990 {
1991  std::cout << "FBEngine::BreakTriangle not implemented yet\n";
1992  return MB_FAILURE;
1993 }
1994 // case 2, n1 and n2 on boundary
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 
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( &currentBoundary, 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 
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 
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 
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 
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 
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 }
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 
2945 {
2946  // these are for debugging purposes only
2947  // check the initial tag, then
2948  Tag 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
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
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 }
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
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
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 
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;
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 }
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
3603 {
3604  // these are for debugging purposes only
3605  // check the initial tag, then
3606  Tag 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
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 }
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 }
3789 {
3790  // need to decide first or second vertex
3791  // important for moab
3792 
3793  Range children;
3794  // EntityHandle v1, v2;
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