Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Intx2Mesh.cpp
Go to the documentation of this file.
1 /*
2  * Intx2Mesh.cpp
3  *
4  * Created on: Oct 2, 2012
5  */
6 
7 #include <limits>
8 #include <queue>
9 #include <sstream>
10 //
12 #ifdef MOAB_HAVE_MPI
13 #include "moab/ParallelComm.hpp"
14 #include "MBParallelConventions.h"
16 #endif /* MOAB_HAVE_MPI */
17 #include "MBTagConventions.hpp"
18 #include "moab/GeomUtil.hpp"
19 #include "moab/AdaptiveKDTree.hpp"
20 
21 namespace moab
22 {
23 
24 #ifdef ENABLE_DEBUG
25 int Intx2Mesh::dbg_1 = 0;
26 #endif
27 
29  : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
30  countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), imaskTag( 0 ),
31  tgtConn( NULL ), srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ),
32  my_rank( 0 )
33 #ifdef MOAB_HAVE_MPI
34  ,
35  parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
36 #endif
37  ,
38  max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
39 {
40  gid = mbimpl->globalId_tag();
41 }
42 
44 {
45  // TODO Auto-generated destructor stub
46 #ifdef MOAB_HAVE_MPI
47  if( remote_cells )
48  {
49  delete remote_cells;
50  remote_cells = NULL;
51  }
52 #endif
53 }
54 
56 {
57  Range cells;
58  ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );
59 
60  max_edges = 0; // can be 0 for point clouds
61  for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
62  {
63  EntityHandle cell = *cit;
64  const EntityHandle* conn4;
65  int nnodes = 3;
66  rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
67  if( nnodes > max_edges ) max_edges = nnodes;
68  }
69  // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
70  // global) but it is better to be consistent
71 #ifdef MOAB_HAVE_MPI
72  if( parcomm )
73  {
74  int local_max_edges = max_edges;
75  // now reduce max_edges over all processors
76  int mpi_err =
77  MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
78  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
79  }
80 #endif
81 
82  return MB_SUCCESS;
83 }
84 
86 {
87  ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" );
88  rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" );
89 
90  return MB_SUCCESS;
91 }
92 
94 {
97  if( countTag ) mb->tag_delete( countTag );
98 
99  unsigned char def_data_bit = 0; // unused by default
100  // maybe the tgt tag is better to be deleted every time, and recreated;
101  // or is it easy to set all values to something again? like 0?
102  ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" );
103  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
104  // this is the only call that is potentially NlogN, in the whole method
105  rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
106 
107  // now, create a map from each edge to a list of potential new nodes on a tgt edge
108  // this memory has to be cleaned up
109  // change it to a vector, and use the index in range of tgt edges
110  int indx = 0;
111  extraNodesVec.resize( TgtEdges.size() );
112  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
113  {
114  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
115  extraNodesVec[indx] = nv;
116  }
117 
118  int defaultInt = -1;
119 
120  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
121  &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
122 
123  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
124  &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
125 
126  rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
127 
128  // for each cell in set 1, determine its neigh in set 1 (could be null too)
129  // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
130  rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
131  rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" );
132 
133  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
134  // them each time edges were for sure created before (tgtEdges)
135  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
136  // if we have a tag with this name, it could be of a different size, so delete it
137  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
139  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
140  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
141  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
142  {
143  EntityHandle tgtCell = *rit;
144  int num_nodes = 0;
145  rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
146  // account for padded polygons
147  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
148  num_nodes--;
149 
150  int i = 0;
151  for( i = 0; i < num_nodes; i++ )
152  {
153  EntityHandle v[2] = { tgtConn[i],
154  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
155  std::vector< EntityHandle > adj_entities;
156  rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
157  if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
158  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
159  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
160  // the tag is dense
161  }
162  // now set the value of the tag
163  rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
164  }
165  return MB_SUCCESS;
166 }
167 
168 ErrorCode Intx2Mesh::DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag )
169 {
170  Range cells;
171  ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" );
172 
173  std::vector< EntityHandle > neighbors( max_edges );
174  std::vector< EntityHandle > zeroh( max_edges, 0 );
175  // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
176  // one on src mesh
177  rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" );
178 
179  for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
180  {
181  EntityHandle cell = *cit;
182  int nnodes = 3;
183  // will get the nnodes ordered neighbors;
184  // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
185  const EntityHandle* conn4;
186  rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
187  int nsides = nnodes;
188  // account for possible padded polygons
189  while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
190  nsides--;
191 
192  for( int i = 0; i < nsides; i++ )
193  {
194  EntityHandle v[2];
195  v[0] = conn4[i];
196  v[1] = conn4[( i + 1 ) % nsides];
197  // get all cells adjacent to these 2 vertices on the edge
198  std::vector< EntityHandle > adjcells;
199  std::vector< EntityHandle > cellsInSet;
200  rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" );
201  // now look for the cells contained in the input set;
202  // the input set should be a correct mesh, not overlapping cells, and manifold
203  size_t siz = adjcells.size();
204  for( size_t j = 0; j < siz; j++ )
205  if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
206  siz = cellsInSet.size();
207 
208  if( siz > 2 )
209  {
210  std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
211  << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" ); // non-manifold
212  }
213  if( siz == 1 )
214  {
215  // it must be the border of the input mesh;
216  neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
217  // borders do not appear for a sphere in serial, but they do appear for
218  // parallel processing anyway
219  continue;
220  }
221  // here siz ==2, it is either the first or second
222  if( cell == cellsInSet[0] )
223  neighbors[i] = cellsInSet[1];
224  else
225  neighbors[i] = cellsInSet[0];
226  }
227  // fill the rest with 0
228  for( int i = nsides; i < max_edges; i++ )
229  neighbors[i] = 0;
230  // now simply set the neighbors tag; the last few positions will not be used, but for
231  // simplicity will keep them all (MAXEDGES)
232  rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
233  }
234  return MB_SUCCESS;
235 }
236 
237 // slow interface; this will not do the advancing front trick
238 // some are triangles, some are quads, some are polygons ...
240 {
241  ErrorCode rval;
242  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
243  mbs2 = mbset2;
244  outSet = outputSet;
245  rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
246  rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
247  // from create tags, copy relevant ones
250  if( countTag ) mb->tag_delete( countTag );
251 
252  // filter rs1 and rs2 by mask; remove everything with 0 mask
253  // get the mask tag if it exists; if not, leave it uninitialized (NULL)
254  rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
255  if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
256  rval = filterByMask( rs1 );MB_CHK_ERR( rval );
257  rval = filterByMask( rs2 );MB_CHK_ERR( rval );
258  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
259  // this is the only call that is potentially NlogN, in the whole method
260  rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
261 
262  int index = 0;
263  extraNodesVec.resize( TgtEdges.size() );
264  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, index++ )
265  {
266  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
267  extraNodesVec[index] = nv;
268  }
269 
270  int defaultInt = -1;
271  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
272  &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
273 
274  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
275  &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
276 
277  rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
278 
279  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
280  // them each time edges were for sure created before (tgtEdges)
281  // if we have a tag with this name, it could be of a different size, so delete it
282  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
284  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
285  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
286  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
287 
288  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
289  {
290  EntityHandle tgtCell = *rit;
291  int num_nodes = 0;
292  rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
293  // account for padded polygons
294  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
295  num_nodes--;
296 
297  for( int i = 0; i < num_nodes; i++ )
298  {
299  EntityHandle v[2] = { tgtConn[i],
300  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
301  std::vector< EntityHandle > adj_entities;
302  rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
303  if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
304  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
305  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
306  // the tag is dense
307  }
308  // now set the value of the tag
309  rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
310  }
311 
312  // find out max edge on source mesh;
313  double max_length = 0;
314  {
315  std::vector< double > coords( 3 * max_edges_1, 0.0 );
316  for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
317  {
318  const EntityHandle* conn = NULL;
319  int nnodes;
320  rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
321  while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
322  nnodes--;
323  rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" );
324  for( int j = 0; j < nnodes; j++ )
325  {
326  int next = ( j + 1 ) % nnodes;
327  double edge_length =
328  ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
329  ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
330  ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
331  if( edge_length > max_length ) max_length = edge_length;
332  }
333  }
334  max_length = std::sqrt( max_length );
335  }
336 
337  // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
338  double tolerance = 1.e-15;
339  if( max_length < 1. )
340  {
341  // basically, the sag for an arc of length max_length on a circle of radius 1
342  tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
344  tolerance = 3 * tolerance; // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
345  // be more generous, use 1.5 ~= sqrt(2.)
346 
347  if( !my_rank )
348  {
349  std::cout << " max edge length: " << max_length << " tolerance for kd tree: " << tolerance << "\n";
350  std::cout << " box overlap tolerance: " << box_error << "\n";
351  }
352  }
353 #ifdef MOAB_HAVE_MPI
354  // reduce box tolerance on every task, if needed
355  double min_box_eps;
356  MPI_Allreduce( &box_error, &min_box_eps, 1, MPI_DOUBLE, MPI_MIN, parcomm->comm() );
357  box_error = min_box_eps;
358 #endif
359 
360  // create the kd tree on source cells, and intersect all targets in an expensive loop
361  // build a kd tree with the rs1 (source) cells
362  FileOptions kdOpts("PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;");
363  AdaptiveKDTree kd( mb );
364  kd.parse_options( kdOpts );
365  EntityHandle tree_root = 0;
366  rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
367 
368  for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
369  {
370  EntityHandle tcell = *it;
371  // find vertex positions
372  const EntityHandle* conn = NULL;
373  int nnodes = 0;
374  rval = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
375  // find leaves close to those positions
376  double areaTgtCell = setup_tgt_cell( tcell, nnodes ); // this is the area in the gnomonic plane
377  double recoveredArea = 0;
378  std::vector< double > positions;
379  positions.resize( nnodes * 3 );
380  rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
381 
382  // distance to search will be based on average edge length
383  double av_len = 0;
384  for( int k = 0; k < nnodes; k++ )
385  {
386  int ik = ( k + 1 ) % nnodes;
387  double len1 = 0;
388  for( int j = 0; j < 3; j++ )
389  {
390  double len2 = positions[3 * k + j] - positions[3 * ik + j];
391  len1 += len2 * len2;
392  }
393  av_len += sqrt( len1 );
394  }
395  if( nnodes > 0 ) av_len /= nnodes;
396  // find leaves within a distance from each vertex of target
397  // in those leaves, collect all cells; we will try for an intx in there
398  Range close_source_cells;
399  std::vector< EntityHandle > leaves;
400  for( int i = 0; i < nnodes; i++ )
401  {
402  leaves.clear();
403  rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
404 
405  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
406  {
407  Range tmp;
408  rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
409 
410  close_source_cells.merge( tmp.begin(), tmp.end() );
411  }
412  }
413 #ifdef VERBOSE
414  if( close_source_cells.empty() )
415  {
416  std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
417  << mb->id_from_handle( tcell ) << "\n";
418  }
419 #endif
420  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
421  {
422  EntityHandle startSrc = *it2;
423  double area = 0;
424  // if area is > 0 , we have intersections
425  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
426  //
427  int nP = 0;
428  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
429  int nsTgt, nsSrc;
430  rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
431  if( area > 0 )
432  {
433  if( nP > 1 )
434  { // this will also construct triangles/polygons in the new mesh, if needed
435  rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
436  }
437  recoveredArea += area;
438  }
439  }
440  recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fract
441  }
442  // before cleaning up , we need to settle the position of the intersection points
443  // on the boundary edges
444  // this needs to be collective, so we should maybe wait something
445 #ifdef MOAB_HAVE_MPI
446  rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
447 #endif
448 
449  this->clean();
450  return MB_SUCCESS;
451 }
452 
453 // main interface; this will do the advancing front trick
454 // some are triangles, some are quads, some are polygons ...
456 {
457  ErrorCode rval;
458  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
459  mbs2 = mbset2;
460  outSet = outputSet;
461 #ifdef VERBOSE
462  std::stringstream ffs, fft;
463  ffs << "source_rank0" << my_rank << ".vtk";
464  rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
465  fft << "target_rank0" << my_rank << ".vtk";
466  rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
467 
468 #endif
469  // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
470 
471  EntityHandle startSrc = 0, startTgt = 0;
472 
473  rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
474  rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
475 
476  // filter rs1 and rs2 by mask; remove everything with 0 mask
477  // get the mask tag if it exists; if not, leave it uninitialized (NULL)
478  mb->tag_get_handle( "GRID_IMASK", imaskTag );
479  if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
480 
481  rval = filterByMask( rs1 );MB_CHK_ERR( rval );
482  rval = filterByMask( rs2 );MB_CHK_ERR( rval );
483 
484  createTags(); // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
485 
486  Range rs22 = rs2; // a copy of the initial range; we will remove from it elements as we
487  // advance ; rs2 is needed for marking the polygon to the tgt parent
488 
489  // create the local kdd tree with source elements; will use it to search
490  // more efficiently for the seeds in advancing front;
491  // some of the target cells will not be covered by source cells, and they need to be eliminated
492  // early from contention
493 
494  // build a kd tree with the rs1 (source) cells
495  FileOptions kdOpts("PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;");
496  AdaptiveKDTree kd( mb );
497  kd.parse_options( kdOpts );
498  EntityHandle tree_root = 0;
499  rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
500 
501  while( !rs22.empty() )
502  {
503 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
504  if( rs22.size() < rs2.size() )
505  {
506  std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
507  << "\n";
508  std::stringstream ffo;
509  ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
510  rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
511  }
512 #endif
513  bool seedFound = false;
514  for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
515  {
516  startTgt = *it;
517  int found = 0;
518  // find vertex positions
519  const EntityHandle* conn = NULL;
520  int nnodes = 0;
521  rval = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
522  // find leaves close to those positions
523  std::vector< double > positions;
524  positions.resize( nnodes * 3 );
525  rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
526  // find leaves within a distance from each vertex of target
527  // in those leaves, collect all cells; we will try for an intx in there, instead of
528  // looping over all rs1 cells, as before
529  Range close_source_cells;
530  std::vector< EntityHandle > leaves;
531  for( int i = 0; i < nnodes; i++ )
532  {
533  leaves.clear();
534  rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
535 
536  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
537  {
538  Range tmp;
539  rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
540 
541  close_source_cells.merge( tmp.begin(), tmp.end() );
542  }
543  }
544 
545  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
546  {
547  startSrc = *it2;
548  double area = 0;
549  // if area is > 0 , we have intersections
550  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
551  //
552  int nP = 0;
553  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
554  int nsTgt, nsSrc;
555  rval =
556  computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
557  if( area > 0 )
558  {
559  found = 1;
560  seedFound = true;
561  break; // found 2 elements that intersect; these will be the seeds
562  }
563  }
564  if( found )
565  break;
566  else
567  {
568 #if defined( VERBOSE )
569  std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
570  << " not intx with any source\n";
571 #endif
572  rs22.erase( startTgt );
573  }
574  }
575  if( !seedFound ) continue; // continue while(!rs22.empty())
576 
577  std::queue< EntityHandle > srcQueue; // these are corresponding to Ta,
578  srcQueue.push( startSrc );
579  std::queue< EntityHandle > tgtQueue;
580  tgtQueue.push( startTgt );
581 
582  Range toResetSrcs; // will be used to reset src flags for every tgt element processed
583 
584  /*if (my_rank==0)
585  dbg_1 = 1;*/
586  unsigned char used = 1;
587  // mark the start tgt quad as used, so it will not come back again
588  rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
589  while( !tgtQueue.empty() )
590  {
591  // flags for the side : 0 means a src cell not found on side
592  // a paired src not found yet for the neighbors of tgt
593  Range nextSrc[MAXEDGES]; // there are new ranges of possible next src cells for
594  // seeding the side j of tgt cell
595 
596  EntityHandle currentTgt = tgtQueue.front();
597  tgtQueue.pop();
598  int nsidesTgt; // will be initialized now
599  double areaTgtCell = setup_tgt_cell( currentTgt, nsidesTgt ); // this is the area in the gnomonic plane
600  double recoveredArea = 0;
601  // get the neighbors of tgt, and if they are solved already, do not bother with that
602  // side of tgt
603  EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
604  rval = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
605 #ifdef ENABLE_DEBUG
606  if( dbg_1 )
607  {
608  std::cout << "Next: neighbors for current tgt ";
609  for( int kk = 0; kk < nsidesTgt; kk++ )
610  {
611  if( tgtNeighbors[kk] > 0 )
612  std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
613  else
614  std::cout << 0 << " ";
615  }
616  std::cout << std::endl;
617  }
618 #endif
619  // now get the status of neighbors; if already solved, make them 0, so not to bother
620  // anymore on that side of tgt
621  for( int j = 0; j < nsidesTgt; j++ )
622  {
623  EntityHandle tgtNeigh = tgtNeighbors[j];
624  unsigned char status = 1;
625  if( tgtNeigh == 0 ) continue;
626  rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval ); // status 0 is unused
627  if( 1 == status ) tgtNeighbors[j] = 0; // so will not look anymore on this side of tgt
628  }
629 
630 #ifdef ENABLE_DEBUG
631  if( dbg_1 )
632  {
633  std::cout << "reset sources: ";
634  for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
635  std::cout << mb->id_from_handle( *itr ) << " ";
636  std::cout << std::endl;
637  }
638 #endif
639  EntityHandle currentSrc = srcQueue.front();
640  // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
641  // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
642  // pairs;
643  // at every moment, the queue contains pairs of cells that intersect, and they form the
644  // "advancing front"
645  srcQueue.pop();
646  toResetSrcs.clear(); // empty the range of used srcs, will have to be set unused again,
647  // at the end of tgt element processing
648  toResetSrcs.insert( currentSrc );
649  // mb2->set_tag_data
650  std::queue< EntityHandle > localSrc;
651  localSrc.push( currentSrc );
652 #ifdef VERBOSE
653  int countingStart = counting;
654 #endif
655  // will advance-front search in the neighborhood of tgt cell, until we finish processing
656  // all
657  // possible src cells; localSrc queue will contain all possible src cells that cover
658  // the current tgt cell
659  while( !localSrc.empty() )
660  {
661  //
662  EntityHandle srcT = localSrc.front();
663  localSrc.pop();
664  double P[10 * MAXEDGES], area; //
665  int nP = 0;
666  int nb[MAXEDGES] = { 0 };
667  int nr[MAXEDGES] = { 0 };
668 
669  int nsidesSrc; ///
670  // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
671  // intersection points could include the vertices of initial elements
672  // nb [j] = 0 means no intersection on the side j for element src (markers)
673  // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
674  // tgt poly. A potential next poly in the tgt queue is the tgt poly that is
675  // adjacent to this side
676  rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
677  nsidesTgt );MB_CHK_ERR( rval );
678  if( nP > 0 )
679  {
680 #ifdef ENABLE_DEBUG
681  if( dbg_1 )
682  {
683  for( int k = 0; k < 3; k++ )
684  {
685  std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
686  }
687  }
688 #endif
689 
690  // intersection found: output P and original triangles if nP > 2
691  EntityHandle neighbors[MAXEDGES] = { 0 };
692  rval = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
693  if( rval != MB_SUCCESS )
694  {
695  std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
696  return MB_FAILURE;
697  }
698 
699  // add neighbors to the localSrc queue, if they are not marked
700  for( int nn = 0; nn < nsidesSrc; nn++ )
701  {
702  EntityHandle neighbor = neighbors[nn];
703  if( neighbor > 0 && nb[nn] > 0 ) // advance across src boundary nn
704  {
705  if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
706  {
707  localSrc.push( neighbor );
708 #ifdef ENABLE_DEBUG
709  if( dbg_1 )
710  {
711  std::cout << " local src elem " << mb->id_from_handle( neighbor )
712  << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
713  mb->list_entities( &neighbor, 1 );
714  }
715 #endif
716  toResetSrcs.insert( neighbor );
717  }
718  }
719  }
720  // n(find(nc>0))=ac; % ac is starting candidate for neighbor
721  for( int nn = 0; nn < nsidesTgt; nn++ )
722  {
723  if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
724  nextSrc[nn].insert( srcT ); // potential src cell that can intersect
725  // the tgt neighbor nn
726  }
727  if( nP > 1 )
728  { // this will also construct triangles/polygons in the new mesh, if needed
729  rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
730  }
731 
732  recoveredArea += area;
733  }
734 #ifdef ENABLE_DEBUG
735  else if( dbg_1 )
736  {
737  std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
738  << mb->id_from_handle( srcT ) << "\n";
739  }
740 #endif
741  } // end while (!localSrc.empty())
742  recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fraction
743 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
744  if( fabs( recoveredArea ) > epsilon_1 )
745  {
746 #ifdef VERBOSE
747  std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
748  << " fraction error recovery:" << recoveredArea
749  << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
750  << "\n";
751 #endif
752  }
753 #endif
754  // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
755  // mesh)
756  rs22.erase( currentTgt );
757  // also, look at its neighbors, and add to the seeds a next one
758 
759  for( int j = 0; j < nsidesTgt; j++ )
760  {
761  EntityHandle tgtNeigh = tgtNeighbors[j];
762  if( tgtNeigh == 0 || nextSrc[j].size() == 0 ) // if tgt is bigger than src, there could be no src
763  // to advance on that side
764  continue;
765  int nsidesTgt2 = 0;
766  setup_tgt_cell( tgtNeigh,
767  nsidesTgt2 ); // find possible intersection with src cell from nextSrc
768  for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
769  {
770  EntityHandle nextB = *nit;
771  // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
772  // src quad
773  double P[10 * MAXEDGES], area; //
774  int nP = 0;
775  int nb[MAXEDGES] = { 0 };
776  int nr[MAXEDGES] = { 0 };
777 
778  int nsidesSrc; ///
780  /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
781  if( area > 0 )
782  {
783  tgtQueue.push( tgtNeigh );
784  srcQueue.push( nextB );
785 #ifdef ENABLE_DEBUG
786  if( dbg_1 )
787  std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
788  << mb->id_from_handle( nextB ) << std::endl;
789 #endif
790  rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
791  break; // so we are done with this side of tgt, we have found a proper next
792  // seed
793  }
794  }
795  }
796 
797  } // end while (!tgtQueue.empty())
798  }
799 #ifdef ENABLE_DEBUG
800  if( dbg_1 )
801  {
802  for( int k = 0; k < 6; k++ )
803  mout_1[k].close();
804  }
805 #endif
806  // before cleaning up , we need to settle the position of the intersection points
807  // on the boundary edges
808  // this needs to be collective, so we should maybe wait something
809 #ifdef MOAB_HAVE_MPI
810  rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
811 #endif
812 
813  this->clean();
814  return MB_SUCCESS;
815 }
816 
818 {
819  if( !imaskTag ) return MB_SUCCESS; // nothing to do
820  size_t sz = cells.size();
821  std::vector< int > masks( sz );
822 
823  ErrorCode rval = mb->tag_get_data( imaskTag, cells, &masks[0] );MB_CHK_ERR( rval );
824  Range cellsToRemove;
825  size_t indx = 0;
826  for( Range::iterator eit = cells.begin(); eit != cells.end(); ++eit, ++indx )
827  {
828  if( masks[indx] ) continue;
829  cellsToRemove.insert( *eit );
830  }
831  cells = subtract( cells, cellsToRemove );
832  return MB_SUCCESS;
833 }
834 
835 // clean some memory allocated
837 {
838  //
839  int indx = 0;
840  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
841  {
842  delete extraNodesVec[indx];
843  }
844  // extraNodesMap.clear();
845  extraNodesVec.clear();
846  // also, delete some bit tags, used to mark processed tgts and srcs
848  counting = 0; // reset counting to original value
849 }
850 
851 // this method will reduce number of nodes, collapse edges that are of length 0
852 // so a polygon like 428 431 431 will become a line 428 431
853 // or something like 428 431 431 531 -> 428 431 531
855 {
856  int i = 0;
857  while( i < nP )
858  {
859  int nextIndex = ( i + 1 ) % nP;
860  if( nodes[i] == nodes[nextIndex] )
861  {
862 #ifdef ENABLE_DEBUG
863  // we need to reduce nP, and collapse nodes
864  if( dbg_1 )
865  {
866  std::cout << " nodes duplicated in list: ";
867  for( int j = 0; j < nP; j++ )
868  std::cout << nodes[j] << " ";
869  std::cout << "\n";
870  std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
871  << "\n";
872  }
873 #endif
874  // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
875  // then next thing does nothing
876  // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
877  for( int k = i; k < nP - 1; k++ )
878  nodes[k] = nodes[k + 1];
879  nP--; // decrease the number of nodes; also, decrease i, just if we may need to check
880  // again
881  i--;
882  }
883  i++;
884  }
885  return;
886 }
887 
888 #ifdef MOAB_HAVE_MPI
889 
890 ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic )
891 {
892  // if it comes here, we want regular 3d boxes
893  // need to refactor this code
894  if( gnomonic ) gnomonic = false;
895  localEnts.clear();
896  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
897 
898  rval = mb->get_connectivity( localEnts, local_verts );
899  int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
900 
901  assert( parcomm != NULL );
902 
903  // get the position of local vertices, and decide local boxes (allBoxes...)
904  double bmin[3] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
905  double bmax[3] = { -std::numeric_limits<double>::max(), -std::numeric_limits<double>::max(), -std::numeric_limits<double>::max() };
906 
907  std::vector< double > coords( 3 * num_local_verts );
908  rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " );
909 
910  for( int i = 0; i < num_local_verts; i++ )
911  {
912  for( int k = 0; k < 3; k++ )
913  {
914  double val = coords[3 * i + k];
915  if( val < bmin[k] ) bmin[k] = val;
916  if( val > bmax[k] ) bmax[k] = val;
917  }
918  }
919  int numprocs = parcomm->proc_config().proc_size();
920  allBoxes.resize( 6 * numprocs );
921 
922  my_rank = parcomm->proc_config().proc_rank();
923  for( int k = 0; k < 3; k++ )
924  {
925  allBoxes[6 * my_rank + k] = bmin[k];
926  allBoxes[6 * my_rank + 3 + k] = bmax[k];
927  }
928 
929  // now communicate to get all boxes
930  int mpi_err;
931 #if( MPI_VERSION >= 2 )
932  // use "in place" option
933  mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
934  parcomm->proc_config().proc_comm() );
935 #else
936  {
937  std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
938  mpi_err = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
939  parcomm->proc_config().proc_comm() );
940  allBoxes = allBoxes_tmp;
941  }
942 #endif
943  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
944 
945 #ifdef VERBOSE
946  if( my_rank == 0 )
947  {
948  std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
949  << " on second mesh \n";
950  for( int i = 0; i < numprocs; i++ )
951  {
952  std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " "
953  << allBoxes[6 * i + 2] << " \n";
954  std::cout << " box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " "
955  << allBoxes[6 * i + 5] << " \n";
956  }
957  }
958 #endif
959 
960  return MB_SUCCESS;
961 }
962 
964 {
965  // compute the bounding box on each proc
966  assert( parcomm != NULL );
967 
968  localEnts.clear();
969  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
970 
971  Tag dpTag = 0;
972  std::string tag_name( "DP" );
973  rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" );
974 
975  EntityHandle dum = 0;
976  Tag corrTag;
977  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );ERRORR( rval, "can't get CORR tag" );
978  // get all local verts
979  Range local_verts;
980  rval = mb->get_connectivity( localEnts, local_verts );
981  int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
982 
983  rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" );
984 
985  std::vector< int > gids( num_local_verts );
986  rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
987 
988  // now see the departure points; to what boxes should we send them?
989  std::vector< double > dep_points( 3 * num_local_verts );
990  rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" );
991  // ranges to send to each processor; will hold vertices and elements (quads?)
992  // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
993  std::map< int, Range > Rto;
994  int numprocs = parcomm->proc_config().proc_size();
995 
996  for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit )
997  {
998  EntityHandle q = *eit;
999  const EntityHandle* conn4;
1000  int num_nodes;
1001  rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
1002  CartVect qbmin( std::numeric_limits<double>::max() );
1003  CartVect qbmax( -std::numeric_limits<double>::max() );
1004  for( int i = 0; i < num_nodes; i++ )
1005  {
1006  EntityHandle v = conn4[i];
1007  size_t index = local_verts.find( v ) - local_verts.begin();
1008  CartVect dp( &dep_points[3 * index] ); // will use constructor
1009  for( int j = 0; j < 3; j++ )
1010  {
1011  if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
1012  if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
1013  }
1014  }
1015  for( int p = 0; p < numprocs; p++ )
1016  {
1017  CartVect bbmin( &allBoxes[6 * p] );
1018  CartVect bbmax( &allBoxes[6 * p + 3] );
1019  if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
1020  {
1021  Rto[p].insert( q );
1022  }
1023  }
1024  }
1025 
1026  // now, build TLv and TLq, for each p
1027  size_t numq = 0;
1028  size_t numv = 0;
1029  for( int p = 0; p < numprocs; p++ )
1030  {
1031  if( p == (int)my_rank ) continue; // do not "send" it, because it is already here
1032  Range& range_to_P = Rto[p];
1033  // add the vertices to it
1034  if( range_to_P.empty() ) continue; // nothing to send to proc p
1035  Range vertsToP;
1036  rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
1037  numq = numq + range_to_P.size();
1038  numv = numv + vertsToP.size();
1039  range_to_P.merge( vertsToP );
1040  }
1041  TupleList TLv;
1042  TupleList TLq;
1043  TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points
1044  TLv.enableWriteAccess();
1045 
1046  int sizeTuple = 2 + max_edges_1; // determined earlier, for src, first mesh
1047  TLq.initialize( 2 + max_edges_1, 0, 1, 0,
1048  numq ); // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh
1049  TLq.enableWriteAccess();
1050 #ifdef VERBOSE
1051  std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1052 #endif
1053  for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1054  {
1055  if( to_proc == (int)my_rank ) continue;
1056  Range& range_to_P = Rto[to_proc];
1057  Range V = range_to_P.subset_by_type( MBVERTEX );
1058 
1059  for( Range::iterator it = V.begin(); it != V.end(); ++it )
1060  {
1061  EntityHandle v = *it;
1062  unsigned int index = local_verts.find( v ) - local_verts.begin();
1063  int n = TLv.get_n();
1064  TLv.vi_wr[2 * n] = to_proc; // send to processor
1065  TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range
1066  TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i]
1067  TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1068  TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1069  TLv.inc_n();
1070  }
1071  // also, prep the quad for sending ...
1072  Range Q = range_to_P.subset_by_dimension( 2 );
1073  for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
1074  {
1075  EntityHandle q = *it;
1076  int global_id;
1077  rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
1078  int n = TLq.get_n();
1079  TLq.vi_wr[sizeTuple * n] = to_proc; //
1080  TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ...
1081  const EntityHandle* conn4;
1082  int num_nodes;
1083  rval = mb->get_connectivity( q, conn4,
1084  num_nodes ); // could be up to MAXEDGES, but it is limited by max_edges_1
1085  ERRORR( rval, "can't get connectivity for cell" );
1086  if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
1087  for( int i = 0; i < num_nodes; i++ )
1088  {
1089  EntityHandle v = conn4[i];
1090  unsigned int index = local_verts.find( v ) - local_verts.begin();
1091  TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1092  }
1093  for( int k = num_nodes; k < max_edges_1; k++ )
1094  {
1095  TLq.vi_wr[sizeTuple * n + 2 + k] =
1096  0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1097  }
1098  TLq.vul_wr[n] = q; // save here the entity handle, it will be communicated back
1099  // maybe we should forget about global ID
1100  TLq.inc_n();
1101  }
1102  }
1103 
1104  // now we are done populating the tuples; route them to the appropriate processors
1105  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1106  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1107  // the elements are already in localEnts;
1108 
1109  // maps from global ids to new vertex and quad handles, that are added
1110  std::map< int, EntityHandle > globalID_to_handle;
1111  /*std::map<int, EntityHandle> globalID_to_eh;*/
1112  globalID_to_eh.clear(); // need for next iteration
1113  // now, look at every TLv, and see if we have to create a vertex there or not
1114  int n = TLv.get_n(); // the size of the points received
1115  for( int i = 0; i < n; i++ )
1116  {
1117  int globalId = TLv.vi_rd[2 * i + 1];
1118  if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1119  {
1120  EntityHandle new_vert;
1121  double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1122  rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
1123  globalID_to_handle[globalId] = new_vert;
1124  }
1125  }
1126 
1127  // now, all dep points should be at their place
1128  // look in the local list of q for this proc, and create all those quads and vertices if needed
1129  // it may be an overkill, but because it does not involve communication, we do it anyway
1130  Range& local = Rto[my_rank];
1131  Range local_q = local.subset_by_dimension( 2 );
1132  // the local should have all the vertices in local_verts
1133  for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1134  {
1135  EntityHandle q = *it;
1136  int nnodes;
1137  const EntityHandle* conn4;
1138  rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " );
1139  EntityHandle new_conn[MAXEDGES];
1140  for( int i = 0; i < nnodes; i++ )
1141  {
1142  EntityHandle v1 = conn4[i];
1143  unsigned int index = local_verts.find( v1 ) - local_verts.begin();
1144  int globalId = gids[index];
1145  if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1146  {
1147  // we need to create that vertex, at this position dep_points
1148  double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
1149  EntityHandle new_vert;
1150  rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
1151  globalID_to_handle[globalId] = new_vert;
1152  }
1153  new_conn[i] = globalID_to_handle[gids[index]];
1154  }
1155  EntityHandle new_element;
1156  //
1157  EntityType entType = MBQUAD;
1158  if( nnodes > 4 ) entType = MBPOLYGON;
1159  if( nnodes < 4 ) entType = MBTRI;
1160 
1161  rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " );
1162  rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
1163  int gid_el;
1164  // get the global ID of the initial quad
1165  rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
1166  globalID_to_eh[gid_el] = new_element;
1167  // is this redundant or not?
1168  rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" );
1169  // set the global id on new elem
1170  rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" );
1171  }
1172  // now look at all elements received through; we do not want to duplicate them
1173  n = TLq.get_n(); // number of elements received by this processor
1174  // form the remote cells, that will be used to send the tracer info back to the originating proc
1175  remote_cells = new TupleList();
1176  remote_cells->initialize( 2, 0, 1, 0, n ); // will not have tracer data anymore
1177  remote_cells->enableWriteAccess();
1178  for( int i = 0; i < n; i++ )
1179  {
1180  int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1181  int from_proc = TLq.vi_wr[sizeTuple * i];
1182  // do we already have a quad with this global ID, represented?
1183  if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
1184  {
1185  // construct the conn quad
1186  EntityHandle new_conn[MAXEDGES];
1187  int nnodes = -1;
1188  for( int j = 0; j < max_edges_1; j++ )
1189  {
1190  int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1191  if( vgid == 0 )
1192  new_conn[j] = 0;
1193  else
1194  {
1195  assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
1196  new_conn[j] = globalID_to_handle[vgid];
1197  nnodes = j + 1; // nodes are at the beginning, and are variable number
1198  }
1199  }
1200  EntityHandle new_element;
1201  //
1202  EntityType entType = MBQUAD;
1203  if( nnodes > 4 ) entType = MBPOLYGON;
1204  if( nnodes < 4 ) entType = MBTRI;
1205  rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
1206  globalID_to_eh[globalIdEl] = new_element;
1207  rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
1208  /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/
1209  remote_cells->vi_wr[2 * i] = from_proc;
1210  remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1211  // remote_cells->vr_wr[i] = 0.; // no contribution yet sent back
1212  remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival)
1213  remote_cells->inc_n();
1214  // set the global id on new elem
1215  rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" );
1216  }
1217  }
1218  // order the remote cells tuple list, with the global id, because we will search in it
1219  // remote_cells->print("remote_cells before sorting");
1220  moab::TupleList::buffer sort_buffer;
1221  sort_buffer.buffer_init( n );
1222  remote_cells->sort( 1, &sort_buffer );
1223  sort_buffer.reset();
1224  return MB_SUCCESS;
1225 }
1226 
1227 // this algorithm assumes lagr set is already created, and some elements will be coming from
1228 // other procs, and populate the covering_set
1229 // we need to keep in a tuple list the remote cells from other procs, because we need to send back
1230 // the intersection info (like area of the intx polygon, and the current concentration) maybe total
1231 // mass in that intx
1233 {
1234  EntityHandle dum = 0;
1235 
1236  Tag corrTag;
1238  // start copy from 2nd alg
1239  // compute the bounding box on each proc
1240  assert( parcomm != NULL );
1241  if( 1 == parcomm->proc_config().proc_size() )
1242  {
1243  covering_set = lagr_set; // nothing to communicate, it must be serial
1244  return MB_SUCCESS;
1245  }
1246 
1247  // get all local verts
1248  Range local_verts;
1249  rval = mb->get_connectivity( localEnts, local_verts );
1250  int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
1251 
1252  std::vector< int > gids( num_local_verts );
1253  rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
1254 
1255  Range localDepCells;
1256  rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" );
1257 
1258  // get all lagr verts (departure vertices)
1259  Range lagr_verts;
1260  rval = mb->get_connectivity( localDepCells, lagr_verts ); // they should be created in
1261  // the same order as the euler vertices
1262  int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" );
1263 
1264  // now see the departure points position; to what boxes should we send them?
1265  std::vector< double > dep_points( 3 * num_lagr_verts );
1266  rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" );
1267  // ranges to send to each processor; will hold vertices and elements (quads?)
1268  // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
1269  std::map< int, Range > Rto;
1270  int numprocs = parcomm->proc_config().proc_size();
1271 
1272  for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
1273  {
1274  EntityHandle q = *eit;
1275  const EntityHandle* conn4;
1276  int num_nodes;
1277  rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
1278  CartVect qbmin( std::numeric_limits<double>::max() );
1279  CartVect qbmax( -std::numeric_limits<double>::max() );
1280  for( int i = 0; i < num_nodes; i++ )
1281  {
1282  EntityHandle v = conn4[i];
1283  int index = lagr_verts.index( v );
1284  assert( -1 != index );
1285  CartVect dp( &dep_points[3 * index] ); // will use constructor
1286  for( int j = 0; j < 3; j++ )
1287  {
1288  if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
1289  if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
1290  }
1291  }
1292  for( int p = 0; p < numprocs; p++ )
1293  {
1294  CartVect bbmin( &allBoxes[6 * p] );
1295  CartVect bbmax( &allBoxes[6 * p + 3] );
1296  if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
1297  {
1298  Rto[p].insert( q );
1299  }
1300  }
1301  }
1302 
1303  // now, build TLv and TLq, for each p
1304  size_t numq = 0;
1305  size_t numv = 0;
1306  for( int p = 0; p < numprocs; p++ )
1307  {
1308  if( p == (int)my_rank ) continue; // do not "send" it, because it is already here
1309  Range& range_to_P = Rto[p];
1310  // add the vertices to it
1311  if( range_to_P.empty() ) continue; // nothing to send to proc p
1312  Range vertsToP;
1313  rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
1314  numq = numq + range_to_P.size();
1315  numv = numv + vertsToP.size();
1316  range_to_P.merge( vertsToP );
1317  }
1318  TupleList TLv;
1319  TupleList TLq;
1320  TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points
1321  TLv.enableWriteAccess();
1322 
1323  int sizeTuple = 2 + max_edges_1; // max edges could be up to MAXEDGES :) for polygons
1324  TLq.initialize( 2 + max_edges_1, 0, 1, 0,
1325  numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v)
1326  // send also the corresponding tgt cell it will come to
1327  TLq.enableWriteAccess();
1328 #ifdef VERBOSE
1329  std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1330 #endif
1331 
1332  for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1333  {
1334  if( to_proc == (int)my_rank ) continue;
1335  Range& range_to_P = Rto[to_proc];
1336  Range V = range_to_P.subset_by_type( MBVERTEX );
1337 
1338  for( Range::iterator it = V.begin(); it != V.end(); ++it )
1339  {
1340  EntityHandle v = *it;
1341  int index = lagr_verts.index( v ); // will be the same index as the corresponding vertex in euler verts
1342  assert( -1 != index );
1343  int n = TLv.get_n();
1344  TLv.vi_wr[2 * n] = to_proc; // send to processor
1345  TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range
1346  TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i]
1347  TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1348  TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1349  TLv.inc_n();
1350  }
1351  // also, prep the 2d cells for sending ...
1352  Range Q = range_to_P.subset_by_dimension( 2 );
1353  for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
1354  {
1355  EntityHandle q = *it; // this is a src cell
1356  int global_id;
1357  rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
1358  int n = TLq.get_n();
1359  TLq.vi_wr[sizeTuple * n] = to_proc; //
1360  TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ...
1361  const EntityHandle* conn4;
1362  int num_nodes;
1363  rval = mb->get_connectivity(
1364  q, conn4, num_nodes ); // could be up to 10;ERRORR( rval, "can't get connectivity for quad" );
1365  if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
1366  for( int i = 0; i < num_nodes; i++ )
1367  {
1368  EntityHandle v = conn4[i];
1369  int index = lagr_verts.index( v );
1370  assert( -1 != index );
1371  TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1372  }
1373  for( int k = num_nodes; k < max_edges_1; k++ )
1374  {
1375  TLq.vi_wr[sizeTuple * n + 2 + k] =
1376  0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1377  }
1378  EntityHandle tgtCell;
1379  rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" );
1380  TLq.vul_wr[n] = tgtCell; // this will be sent to remote_cells, to be able to come back
1381  TLq.inc_n();
1382  }
1383  }
1384  // now we can route them to each processor
1385  // now we are done populating the tuples; route them to the appropriate processors
1386  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1387  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1388  // the elements are already in localEnts;
1389 
1390  // maps from global ids to new vertex and quad handles, that are added
1391  std::map< int, EntityHandle > globalID_to_handle;
1392  // we already have vertices from lagr set; they are already in the processor, even before
1393  // receiving other verts from neighbors
1394  int k = 0;
1395  for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
1396  {
1397  globalID_to_handle[gids[k]] = *vit; // a little bit of overkill
1398  // we do know that the global ids between euler and lagr verts are parallel
1399  }
1400  /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one?
1401  globalID_to_eh.clear();
1402  // now, look at every TLv, and see if we have to create a vertex there or not
1403  int n = TLv.get_n(); // the size of the points received
1404  for( int i = 0; i < n; i++ )
1405  {
1406  int globalId = TLv.vi_rd[2 * i + 1];
1407  if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1408  {
1409  EntityHandle new_vert;
1410  double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1411  rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
1412  globalID_to_handle[globalId] = new_vert;
1413  }
1414  }
1415 
1416  // now, all dep points should be at their place
1417  // look in the local list of 2d cells for this proc, and create all those cells if needed
1418  // it may be an overkill, but because it does not involve communication, we do it anyway
1419  Range& local = Rto[my_rank];
1420  Range local_q = local.subset_by_dimension( 2 );
1421  // the local should have all the vertices in lagr_verts
1422  for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1423  {
1424  EntityHandle q = *it; // these are from lagr cells, local
1425  int gid_el;
1426  rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
1427  globalID_to_eh[gid_el] = q; // do we need this? maybe to just mark the ones on this processor
1428  // maybe a range of global cell ids is fine?
1429  }
1430  // now look at all elements received through; we do not want to duplicate them
1431  n = TLq.get_n(); // number of elements received by this processor
1432  // a cell should be received from one proc only; so why are we so worried about duplicated
1433  // elements? a vertex can be received from multiple sources, that is fine
1434 
1435  remote_cells = new TupleList();
1436  remote_cells->initialize( 2, 0, 1, 0, n ); // no tracers anymore in these tuples
1437  remote_cells->enableWriteAccess();
1438  for( int i = 0; i < n; i++ )
1439  {
1440  int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1441  int from_proc = TLq.vi_rd[sizeTuple * i];
1442  // do we already have a quad with this global ID, represented?
1443  if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
1444  {
1445  // construct the conn quad
1446  EntityHandle new_conn[MAXEDGES];
1447  int nnodes = -1;
1448  for( int j = 0; j < max_edges_1; j++ )
1449  {
1450  int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1451  if( vgid == 0 )
1452  new_conn[j] = 0;
1453  else
1454  {
1455  assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
1456  new_conn[j] = globalID_to_handle[vgid];
1457  nnodes = j + 1; // nodes are at the beginning, and are variable number
1458  }
1459  }
1460  EntityHandle new_element;
1461  //
1462  EntityType entType = MBQUAD;
1463  if( nnodes > 4 ) entType = MBPOLYGON;
1464  if( nnodes < 4 ) entType = MBTRI;
1465  rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
1466  globalID_to_eh[globalIdEl] = new_element;
1467  local_q.insert( new_element );
1468  rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " );
1469  }
1470  remote_cells->vi_wr[2 * i] = from_proc;
1471  remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1472  // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
1473  remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival)
1474  remote_cells->inc_n();
1475  }
1476  // now, create a new set, covering_set
1477  rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " );
1478  rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " );
1479  // order the remote cells tuple list, with the global id, because we will search in it
1480  // remote_cells->print("remote_cells before sorting");
1481  moab::TupleList::buffer sort_buffer;
1482  sort_buffer.buffer_init( n );
1483  remote_cells->sort( 1, &sort_buffer );
1484  sort_buffer.reset();
1485  return MB_SUCCESS;
1486  // end copy
1487 }
1488 
1489 ErrorCode Intx2Mesh::resolve_intersection_sharing()
1490 {
1491  if( parcomm && parcomm->size() > 1 )
1492  {
1493  /*
1494  moab::ParallelMergeMesh pm(parcomm, epsilon_1);
1495  ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip
1496  local merge, use dim 2 ERRORR(rval, "can't merge intersection ");
1497  */
1498  // look at non-owned shared vertices, that could be part of original source set
1499  // they should be removed from intx set reference, because they might not have a
1500  // correspondent on the other task
1501  Range nonOwnedVerts;
1502  Range vertsInIntx;
1503  Range intxCells;
1504  ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval );
1505  rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval );
1506 
1507  rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval );
1508 
1509  // some of these vertices can be in original set 1, which was covered, transported;
1510  // but they should not be "shared" from the intx point of view, because they are not shared
1511  // with another task they might have come from coverage as a plain vertex, so losing the
1512  // sharing property ?
1513 
1514  Range coverVerts;
1515  rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval );
1516  // find out those that are on the interface
1517  Range vertsCovInterface;
1518  rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval );
1519  // how many of these are in
1520  Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts );
1521  // first, get all cells connected to these vertices, from intxCells
1522 
1523  Range connectedCells;
1524  rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval );
1525  // only those in intx set:
1526  connectedCells = intersect( connectedCells, intxCells );
1527  // first duplicate vertices in question:
1528  std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
1529  for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
1530  {
1531  EntityHandle vertex = *vit;
1532  double coords[3];
1533  rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval );
1534  EntityHandle newVertex;
1535  rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval );
1536  duplicatedVerticesMap[vertex] = newVertex;
1537  }
1538 
1539  // look now at connectedCells, and change their connectivities:
1540  for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
1541  {
1542  EntityHandle intxCell = *eit;
1543  // replace connectivity
1544  std::vector< EntityHandle > connectivity;
1545  rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval );
1546  for( size_t i = 0; i < connectivity.size(); i++ )
1547  {
1548  EntityHandle currentVertex = connectivity[i];
1549  std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
1550  if( mit != duplicatedVerticesMap.end() )
1551  {
1552  connectivity[i] = mit->second; // replace connectivity directly
1553  }
1554  }
1555  int nnodes = (int)connectivity.size();
1556  rval = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval );
1557  }
1558  }
1559  return MB_SUCCESS;
1560 }
1561 #endif /* MOAB_HAVE_MPI */
1562 } /* namespace moab */