Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Intx2MeshOnSphere.cpp
Go to the documentation of this file.
1 /*
2  * Intx2MeshOnSphere.cpp
3  *
4  * Created on: Oct 3, 2012
5  */
6 
7 #ifdef _MSC_VER /* windows */
8 #define _USE_MATH_DEFINES // For M_PI
9 #endif
10 
13 #include "moab/GeomUtil.hpp"
14 #include "moab/BoundBox.hpp"
15 #include "moab/MeshTopoUtil.hpp"
16 #ifdef MOAB_HAVE_MPI
17 #include "moab/ParallelComm.hpp"
18 #endif
19 #include "MBTagConventions.hpp"
20 
21 #include <cassert>
22 
23 // #define ENABLE_DEBUG
24 #define CHECK_CONVEXITY
25 namespace moab
26 {
27 
29  : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
30 {
31 }
32 
34 
35 /*
36  * return also the area for robustness verification
37  */
39 {
40 
41  // get coordinates of the target quad, to decide the gnomonic plane
42  double cellArea = 0;
43 
44  int num_nodes;
45  ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
46 
47  nsTgt = num_nodes;
48  // account for possible padded polygons
49  while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
50  nsTgt--;
51 
52  // CartVect coords[4];
53  rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
54 
55  CartVect middle = tgtCoords[0];
56  for( int i = 1; i < nsTgt; i++ )
57  middle += tgtCoords[i];
58  middle = 1. / nsTgt * middle;
59 
60  IntxUtils::decide_gnomonic_plane( middle, plane ); // output the plane
61  for( int j = 0; j < nsTgt; j++ )
62  {
63  // populate coords in the plane for intersection
64  // they should be oriented correctly, positively
65  rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
66  }
67 
68  for( int j = 1; j < nsTgt - 1; j++ )
69  cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
70 
71  // take target coords in order and compute area in plane
72  return cellArea;
73 }
74 
75 /* the elements are convex for sure, then do a gnomonic projection of both,
76  * compute intersection in the plane, then go back to the sphere for the points
77  * */
79  EntityHandle src,
80  double* P,
81  int& nP,
82  double& area,
83  int markb[MAXEDGES],
84  int markr[MAXEDGES],
85  int& nsBlue,
86  int& nsTgt,
87  bool check_boxes_first )
88 {
89  // the area will be used from now on, to see how well we fill the target cell with polygons
90  // the points will be at most 40; they will describe a convex patch, after the points will be
91  // ordered and collapsed (eliminate doubles)
92 
93  // CartVect srccoords[4];
94  int num_nodes = 0;
95  ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
96  nsBlue = num_nodes;
97  // account for possible padded polygons
98  while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
99  nsBlue--;
100  rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
101 
102  area = 0.;
103  nP = 0; // number of intersection points we are marking the boundary of src!
104  if( check_boxes_first )
105  {
106  // look at the boxes formed with vertices; if they are far away, return false early
107  // make sure the target is setup already
108  setup_tgt_cell( tgt, nsTgt ); // we do not need area here
109  // use here gnomonic plane (plane) to see where source is
110  bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
111  int planeb;
112  CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
113  IntxUtils::decide_gnomonic_plane( mid3, planeb );
114  if( !overlap3d && ( plane != planeb ) ) // plane was set at setup_tgt_cell
115  return MB_SUCCESS; // no error, but no intersection, decide early to get out
116  // if same plane, still check for gnomonic plane in 2d
117  // if no overlap in 2d, get out
118  if( !overlap3d && plane == planeb ) // CHECK 2D too
119  {
120  for( int j = 0; j < nsBlue; j++ )
121  {
123  srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
124  }
125  bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
126  if( !overlap2d ) return MB_SUCCESS; // we are sure they are not overlapping in 2d , either
127  }
128  }
129 #ifdef ENABLE_DEBUG
130  if( dbg_1 )
131  {
132  std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
133  for( int j = 0; j < nsTgt; j++ )
134  {
135  std::cout << tgtCoords[j] << "\n";
136  }
137  std::cout << "src " << mb->id_from_handle( src ) << "\n";
138  for( int j = 0; j < nsBlue; j++ )
139  {
140  std::cout << srcCoords[j] << "\n";
141  }
142  mb->list_entities( &tgt, 1 );
143  mb->list_entities( &src, 1 );
144  }
145 #endif
146 
147  for( int j = 0; j < nsBlue; j++ )
148  {
149  rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
150  }
151 
152 #ifdef ENABLE_DEBUG
153  if( dbg_1 )
154  {
155  std::cout << "gnomonic plane: " << plane << "\n";
156  std::cout << " target src\n";
157  for( int j = 0; j < nsTgt; j++ )
158  {
159  std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
160  }
161  for( int j = 0; j < nsBlue; j++ )
162  {
163  std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
164  }
165  }
166 #endif
167 
168  rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
169 
170  int side[MAXEDGES] = { 0 }; // this refers to what side? source or tgt?
171  int extraPoints =
172  IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
173  if( extraPoints >= 1 )
174  {
175  for( int k = 0; k < nsBlue; k++ )
176  {
177  if( side[k] )
178  {
179  // this means that vertex k of source is inside convex tgt; mark edges k-1 and k in
180  // src,
181  // as being "intersected" by tgt; (even though they might not be intersected by
182  // other edges, the fact that their apex is inside, is good enough)
183  markb[k] = 1;
184  markb[( k + nsBlue - 1 ) % nsBlue] =
185  1; // it is the previous edge, actually, but instead of doing -1, it is
186  // better to do modulo +3 (modulo 4)
187  // null side b for next call
188  side[k] = 0;
189  }
190  }
191  }
192  nP += extraPoints;
193 
194  extraPoints =
195  IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
196  if( extraPoints >= 1 )
197  {
198  for( int k = 0; k < nsTgt; k++ )
199  {
200  if( side[k] )
201  {
202  // this is to mark that target edges k-1 and k are intersecting src
203  markr[k] = 1;
204  markr[( k + nsTgt - 1 ) % nsTgt] =
205  1; // it is the previous edge, actually, but instead of doing -1, it is
206  // better to do modulo +3 (modulo 4)
207  // null side b for next call
208  }
209  }
210  }
211  nP += extraPoints;
212 
213  // now sort and orient the points in P, such that they are forming a convex polygon
214  // this will be the foundation of our new mesh
215  // this works if the polygons are convex
216  IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 ); // nP should be at most 8 in the end ?
217  // if there are more than 3 points, some area will be positive
218 
219  if( nP >= 3 )
220  {
221  for( int k = 1; k < nP - 1; k++ )
222  area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
223 #ifdef CHECK_CONVEXITY
224  // each edge should be large enough that we can compute angles between edges
225  for( int k = 0; k < nP; k++ )
226  {
227  int k1 = ( k + 1 ) % nP;
228  int k2 = ( k1 + 1 ) % nP;
229  double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
230  if( orientedArea < 0 )
231  {
232  std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
233  << " " << src << " \n";
234  }
235  }
236 #endif
237  }
238 
239  return MB_SUCCESS; // no error
240 }
241 
242 // this method will also construct the triangles/quads/polygons in the new mesh
243 // if we accept planar polygons, we just save them
244 // also, we could just create new vertices every time, and merge only in the end;
245 // could be too expensive, and the tolerance for merging could be an
246 // interesting topic
247 ErrorCode Intx2MeshOnSphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsBlue, double* iP, int nP )
248 {
249 #ifdef ENABLE_DEBUG
250  // first of all, check against target and source vertices
251  //
252  if( dbg_1 )
253  {
254  std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
255  << "\n";
256  for( int n = 0; n < nP; n++ )
257  std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
258  }
259 #endif
260 
261  // get the edges for the target triangle; the extra points will be on those edges, saved as
262  // lists (unordered)
263 
264  // first get the list of edges adjacent to the target cell
265  // use the neighTgtEdgeTag
266  EntityHandle adjTgtEdges[MAXEDGES];
267  ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge target tag" );
268  // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
269  // potatoes some of them will be handles to the initial vertices from source or target meshes
270 
271  std::vector< EntityHandle > foundIds;
272  foundIds.resize( nP );
273 #ifdef CHECK_CONVEXITY
274  int npBefore1 = nP;
275  int oldNodes = 0;
276  int otherIntx = 0;
277  moab::IntxAreaUtils areaAdaptor;
278 #endif
279  for( int i = 0; i < nP; i++ )
280  {
281  double* pp = &iP[2 * i]; // iP+2*i
282  // project the point back on the sphere
283  CartVect pos;
284  IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
285  int found = 0;
286  // first, are they on vertices from target or src?
287  // priority is the target mesh (mb2?)
288  int j = 0;
289  EntityHandle outNode = (EntityHandle)0;
290  for( j = 0; j < nsTgt && !found; j++ )
291  {
292  // int node = tgtTri.v[j];
293  double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
294  if( d2 < epsilon_1 / 1000 ) // two orders of magnitude smaller than it should, to avoid concave polygons
295  {
296 
297  foundIds[i] = tgtConn[j]; // no new node
298  found = 1;
299 #ifdef CHECK_CONVEXITY
300  oldNodes++;
301 #endif
302 #ifdef ENABLE_DEBUG
303  if( dbg_1 )
304  std::cout << " target node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
305  << " 2d coords:" << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << " d2: " << d2
306  << " \n";
307 #endif
308  }
309  }
310 
311  for( j = 0; j < nsBlue && !found; j++ )
312  {
313  // int node = srcTri.v[j];
314  double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
315  if( d2 < epsilon_1 / 1000 )
316  {
317  // suspect is srcConn[j] corresponding in mbOut
318 
319  foundIds[i] = srcConn[j]; // no new node
320  found = 1;
321 #ifdef CHECK_CONVEXITY
322  oldNodes++;
323 #endif
324 #ifdef ENABLE_DEBUG
325  if( dbg_1 )
326  std::cout << " source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
327  << " \n";
328 #endif
329  }
330  }
331 
332  if( !found )
333  {
334  // find the edge it belongs, first, on the red element
335  // look at the minimum area, not at the first below some tolerance
336  double minArea = 1.e+38;
337  int index_min = -1;
338  for( j = 0; j < nsTgt; j++ )
339  {
340  int j1 = ( j + 1 ) % nsTgt;
341  double area = fabs( IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp ) );
342  // how to check if pp is between redCoords2D[j] and redCoords2D[j1] ?
343  // they should form a straight line; the sign should be -1
344  double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
345  IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
346  IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
347  if( area < minArea && checkx < 2 * epsilon_1 ) // round off error or not?
348  {
349  index_min = j;
350  minArea = area;
351  }
352  }
353  // verify that index_min is valid
354  assert( index_min >= 0 );
355 
356  if( minArea < epsilon_1 / 2 ) // we found the smallest area, so we think we found the
357  // target edge it belongs
358  {
359  // found the edge; now find if there is a point in the list here
360  // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
361  int indx = TgtEdges.index( adjTgtEdges[index_min] );
362  if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
363  {
364  std::cerr << " error in adjacent target edge: " << mb->id_from_handle( adjTgtEdges[index_min] )
365  << "\n";
366  return MB_FAILURE;
367  }
368  std::vector< EntityHandle >* expts = extraNodesVec[indx];
369  // if the points pp is between extra points, then just give that id
370  // if not, create a new point, (check the id)
371  // get the coordinates of the extra points so far
372  int nbExtraNodesSoFar = expts->size();
373  if( nbExtraNodesSoFar > 0 )
374  {
375  std::vector< CartVect > coords1;
376  coords1.resize( nbExtraNodesSoFar );
377  mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
378  // std::list<int>::iterator it;
379  for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
380  {
381  // int pnt = *it;
382  double d2 = ( pos - coords1[k] ).length();
383  if( d2 < 2 * epsilon_1 ) // is this below machine precision?
384  {
385  found = 1;
386  foundIds[i] = ( *expts )[k];
387 #ifdef CHECK_CONVEXITY
388  otherIntx++;
389 #endif
390  }
391  }
392  }
393  if( !found )
394  {
395  // create a new point in 2d (at the intersection)
396  // foundIds[i] = m_num2dPoints;
397  // expts.push_back(m_num2dPoints);
398  // need to create a new node in mbOut
399  // this will be on the edge, and it will be added to the local list
400  rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
401  ( *expts ).push_back( outNode );
402  // CID 181168; avoid leak storage error
403  rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
404  foundIds[i] = outNode;
405  found = 1;
406  }
407  }
408  }
409  if( !found )
410  {
411  std::cout << " target quad: ";
412  for( int j1 = 0; j1 < nsTgt; j1++ )
413  {
414  std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
415  }
416  std::cout << " a point pp is not on a target quad " << *pp << " " << pp[1] << " target quad "
417  << mb->id_from_handle( tgt ) << " \n";
418  return MB_FAILURE;
419  }
420  }
421 #ifdef ENABLE_DEBUG
422  if( dbg_1 )
423  {
424  std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
425  for( int i1 = 0; i1 < nP; i1++ )
426  std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
427  }
428 #endif
429  // first, find out if we have nodes collapsed; shrink them
430  // we may have to reduce nP
431  // it is possible that some nodes are collapsed after intersection only
432  // nodes will always be in order (convex intersection)
433 #ifdef CHECK_CONVEXITY
434  int npBefore2 = nP;
435 #endif
436  correct_polygon( &foundIds[0], nP );
437  // now we can build the triangles, from P array, with foundIds
438  // we will put them in the out set
439  if( nP >= 3 )
440  {
441  EntityHandle polyNew;
442  rval = mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew );MB_CHK_ERR( rval );
443  rval = mb->add_entities( outSet, &polyNew, 1 );MB_CHK_ERR( rval );
444 
445  // tag it with the global ids from target and source elements
446  int globalID;
447  rval = mb->tag_get_data( gid, &src, 1, &globalID );MB_CHK_ERR( rval );
448  rval = mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
449  // if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
450  // " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";
451  rval = mb->tag_get_data( gid, &tgt, 1, &globalID );MB_CHK_ERR( rval );
452  rval = mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
453  // if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
454  // " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";
455 
456  counting++;
457  rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
458  if( orgSendProcTag )
459  {
460  int org_proc = -1;
461  rval = mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc );MB_CHK_ERR( rval );
462  rval = mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc );MB_CHK_ERR( rval ); // yet another tag
463  }
464 #ifdef CHECK_CONVEXITY
465  // each edge should be large enough that we can compute angles between edges
466  std::vector< double > coords;
467  coords.resize( 3 * nP );
468  rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
469  std::vector< CartVect > posi( nP );
470  rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
471 
472  for( int k = 0; k < nP; k++ )
473  {
474  int k1 = ( k + 1 ) % nP;
475  int k2 = ( k1 + 1 ) % nP;
476  double orientedArea = areaAdaptor. area_spherical_triangle( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
477  if( orientedArea < 0 )
478  {
479  std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
480  for( int i = 0; i < nP; i++ )
481  {
482  int nexti = ( i + 1 ) % nP;
483  double lengthEdge = ( posi[i] - posi[nexti] ).length();
484  std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n";
485  }
486  std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n";
487 
488  std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
489  << " target, src:" << tgt << " " << src << " \n";
490  }
491  }
492 #endif
493 
494 #ifdef ENABLE_DEBUG
495  if( dbg_1 )
496  {
497  std::cout << "Counting: " << counting << "\n";
498  std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
499  for( int i1 = 0; i1 < nP; i1++ )
500  std::cout << " " << mb->id_from_handle( foundIds[i1] );
501  std::cout << " plane: " << plane << "\n";
502  std::vector< CartVect > posi( nP );
503  mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
504  for( int i1 = 0; i1 < nP; i1++ )
505  std::cout << foundIds[i1] << " " << posi[i1] << "\n";
506 
507  std::stringstream fff;
508  fff << "file0" << counting << ".vtk";
509  rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
510  }
511 #endif
512  }
513  // else {
514  // std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";
515  // }
516  // disable_debug();
517  return MB_SUCCESS;
518 }
519 
521 {
522  EntityHandle dum = 0;
523 
524  Tag corrTag;
526  &dum ); // it should have been created
527  MB_CHK_SET_ERR( rval, "can't get correlation tag" );
528 
529  // get all polygons out of out_set; then see where are they coming from
530  Range polys;
531  rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
532 
533  // rs2 is the target range, arrival; rs1 is src, departure;
534  // there is a connection between rs1 and rs2, through the corrTag
535  // corrTag is __correlation
536  // basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly);
537  // also, mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly);
538  // we start from rs2 existing, then we have to update something
539 
540  // tagElem will have multiple tracers
541  int numTracers = 0;
542  rval = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
543  if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
544 
545  std::vector< double > currentVals( rs2.size() * numTracers );
546  rval = mb->tag_get_data( tagElem, rs2, &currentVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" );
547 
548  // create new tuple list for tracers to other processors, from remote_cells
549 #ifdef MOAB_HAVE_MPI
550  if( remote_cells )
551  {
552  int n = remote_cells->get_n();
553  if( n > 0 )
554  {
555  remote_cells_with_tracers = new TupleList();
556  remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
557  n ); // tracers are in these tuples
558  remote_cells_with_tracers->enableWriteAccess();
559  for( int i = 0; i < n; i++ )
560  {
561  remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i];
562  remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
563  // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
564  remote_cells_with_tracers->vul_wr[i] =
565  remote_cells->vul_wr[i]; // this is the corresponding target cell (arrival)
566  for( int k = 0; k < numTracers; k++ )
567  remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0; // initialize tracers to be transported
568  remote_cells_with_tracers->inc_n();
569  }
570  }
571  delete remote_cells;
572  remote_cells = NULL;
573  }
574 #endif
575  // for each polygon, we have 2 indices: target and source parents
576  // we need index source to update index tgt?
577  std::vector< double > newValues( rs2.size() * numTracers,
578  0. ); // initialize with 0 all of them
579  // area of the polygon * conc on target (old) current quantity
580  // finally, divide by the area of the tgt
581  double check_intx_area = 0.;
582  moab::IntxAreaUtils intxAreas( this->areaMethod ); // use_lHuiller = true
583  for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
584  {
585  EntityHandle poly = *it;
586  int srcIndex, tgtIndex;
587  rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
588 
589  EntityHandle src = rs1[srcIndex - 1]; // big assumption, it should work for meshes where global id is the same
590  // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes
591  rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" );
592  // EntityHandle target = rs2[tgtIndex];
593  // big assumption here, target and source are "parallel" ;we should have an index from
594  // source to target (so a deformed source corresponds to an arrival tgt)
595  /// TODO: VSM: Its unclear whether we need the source or destination radius here.
596  double radius = Rsrc;
597  double areap = intxAreas.area_spherical_element( mb, poly, radius );
598  check_intx_area += areap;
599  // so the departure cell at time t (srcIndex) covers a portion of a tgtCell
600  // that quantity will be transported to the tgtCell at time t+dt
601  // the source corresponds to a target arrival
602  EntityHandle tgtArr;
603  rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
604  if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
605  {
606 #ifdef MOAB_HAVE_MPI
607  if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
608  // maybe the element is remote, from another processor
609  int global_id_src;
610  rval = mb->tag_get_data( gid, &src, 1, &global_id_src );MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding source gid" );
611  // find the
612  int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
613  if( index_in_remote == -1 )
614  MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
615  for( int k = 0; k < numTracers; k++ )
616  remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
617  currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
618 #endif
619  }
620  else if( MB_SUCCESS == rval )
621  {
622  int arrTgtIndex = rs2.index( tgtArr );
623  if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
624  for( int k = 0; k < numTracers; k++ )
625  newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
626  }
627 
628  else
629  MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " );
630  }
631  // now, send back the remote_cells_with_tracers to the processors they came from, with the
632  // updated values for the tracer mass in a cell
633 #ifdef MOAB_HAVE_MPI
634  if( remote_cells_with_tracers )
635  {
636  // so this means that some cells will be sent back with tracer info to the procs they were
637  // sent from
638  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
639  // now, look at the global id, find the proper "tgt" cell with that index and update its
640  // mass
641  // remote_cells->print("remote cells after routing");
642  int n = remote_cells_with_tracers->get_n();
643  for( int j = 0; j < n; j++ )
644  {
645  EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j]; // entity handle sent back
646  int arrTgtIndex = rs2.index( tgtCell );
647  if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
648  for( int k = 0; k < numTracers; k++ )
649  newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
650  }
651  }
652 #endif /* MOAB_HAVE_MPI */
653  // now divide by target area (current)
654  int j = 0;
655  Range::iterator iter = rs2.begin();
656  void* data = NULL; // used for stored area
657  int count = 0;
658  std::vector< double > total_mass_local( numTracers, 0. );
659  while( iter != rs2.end() )
660  {
661  rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
662  double* ptrArea = (double*)data;
663  for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
664  {
665  for( int k = 0; k < numTracers; k++ )
666  {
667  total_mass_local[k] += newValues[j * numTracers + k];
668  newValues[j * numTracers + k] /= ( *ptrArea );
669  }
670  }
671  }
672  rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" );
673 
674 #ifdef MOAB_HAVE_MPI
675  std::vector< double > total_mass( numTracers, 0. );
676  double total_intx_area = 0;
677  int mpi_err =
678  MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
679  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
680  // now reduce total area
681  mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
682  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
683  if( my_rank == 0 )
684  {
685  for( int k = 0; k < numTracers; k++ )
686  std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n";
687  std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
688  << total_intx_area << "\n";
689  }
690 
691  if( remote_cells_with_tracers )
692  {
693  delete remote_cells_with_tracers;
694  remote_cells_with_tracers = NULL;
695  }
696 #else
697  for( int k = 0; k < numTracers; k++ )
698  std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n";
699  std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
700  << check_intx_area << "\n";
701 #endif
702  return MB_SUCCESS;
703 }
704 
705 #ifdef MOAB_HAVE_MPI
706 ErrorCode Intx2MeshOnSphere::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic )
707 {
708  if( !gnomonic )
709  {
710  return Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts, gnomonic );
711  }
712  localEnts.clear();
713  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );MB_CHK_SET_ERR( rval, "can't get local ents" );
714 
715  rval = mb->get_connectivity( localEnts, local_verts );MB_CHK_SET_ERR( rval, "can't get connectivity" );
716  int num_local_verts = (int)local_verts.size();
717 
718  assert( parcomm != NULL );
719 
720  if( num_local_verts == 0 )
721  {
722  // it is probably point cloud, get the local vertices from set
723  rval = mb->get_entities_by_dimension( euler_set, 0, local_verts );MB_CHK_SET_ERR( rval, "can't get local vertices from set" );
724  num_local_verts = (int)local_verts.size();
725  localEnts = local_verts;
726  }
727  // will use 6 gnomonic planes to decide boxes for each gnomonic plane
728  // each gnomonic box will be 2d, min, max
729  double gnom_box[24];
730  for( int i = 0; i < 6; i++ )
731  {
732  gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
733  gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
734  }
735 
736  // there are 6 gnomonic planes; some elements could be on the corners, and affect multiple
737  // planes decide what gnomonic planes will be affected by each cell some elements could appear
738  // in multiple gnomonic planes !
739  std::vector< double > coords( 3 * num_local_verts );
740  rval = mb->get_coords( local_verts, &coords[0] );MB_CHK_SET_ERR( rval, "can't get vertex coords" );ERRORR( rval, "can't get coords of vertices " );
741  // decide each local vertex to what gnomonic plane it belongs
742 
743  std::vector< int > gnplane;
744  gnplane.resize( num_local_verts );
745  for( int i = 0; i < num_local_verts; i++ )
746  {
747  CartVect pos( &coords[3 * i] );
748  int pl;
750  gnplane[i] = pl;
751  }
752 
753  for( Range::iterator it = localEnts.begin(); it != localEnts.end(); it++ )
754  {
755  EntityHandle cell = *it;
756  EntityType typeCell = mb->type_from_handle( cell ); // could be vertex, for point cloud
757  // get coordinates, and decide gnomonic planes for it
758  int nnodes;
759  const EntityHandle* conn = NULL;
760  EntityHandle c[1];
761  if( typeCell != MBVERTEX )
762  {
763  rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
764  }
765  else
766  {
767  nnodes = 1;
768  c[0] = cell; // actual node
769  conn = &c[0];
770  }
771  // get coordinates of vertices involved with this
772  std::vector< double > elco( 3 * nnodes );
773  std::set< int > planes;
774  for( int i = 0; i < nnodes; i++ )
775  {
776  int ix = local_verts.index( conn[i] );
777  planes.insert( gnplane[ix] );
778  for( int j = 0; j < 3; j++ )
779  {
780  elco[3 * i + j] = coords[3 * ix + j];
781  }
782  }
783  // now, augment the boxes for all planes involved
784  for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
785  {
786  int pl = *st;
787  for( int i = 0; i < nnodes; i++ )
788  {
789  CartVect pos( &elco[3 * i] );
790  double c2[2];
791  IntxUtils::gnomonic_projection( pos, Rdest, pl, c2[0],
792  c2[1] ); // 2 coordinates
793  //
794  for( int k = 0; k < 2; k++ )
795  {
796  double val = c2[k];
797  if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val; // min in k direction
798  if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
799  gnom_box[4 * ( pl - 1 ) + 2 + k] = val; // max in k direction
800  }
801  }
802  }
803  }
804 
805  int numprocs = parcomm->proc_config().proc_size();
806  allBoxes.resize( 24 * numprocs ); // 6 gnomonic planes , 4 doubles for each for 2d box
807 
808  my_rank = parcomm->proc_config().proc_rank();
809  for( int k = 0; k < 24; k++ )
810  allBoxes[24 * my_rank + k] = gnom_box[k];
811 
812  // now communicate to get all boxes
813  int mpi_err;
814 #if( MPI_VERSION >= 2 )
815  // use "in place" option
816  mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 24, MPI_DOUBLE,
817  parcomm->proc_config().proc_comm() );
818 #else
819  {
820  std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
821  mpi_err = MPI_Allgather( &allBoxes[24 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
822  parcomm->proc_config().proc_comm() );
823  allBoxes = allBoxes_tmp;
824  }
825 #endif
826  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
827 
828 #ifdef VERBOSE
829  if( my_rank == 0 )
830  {
831  std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
832  << " on second mesh \n";
833  for( int i = 0; i < numprocs; i++ )
834  {
835  std::cout << "task: " << i << " \n";
836  for( int pl = 1; pl <= 6; pl++ )
837  {
838  std::cout << " plane " << pl << " min: \t" << allBoxes[24 * i + 4 * ( pl - 1 )] << " \t"
839  << allBoxes[24 * i + 4 * ( pl - 1 ) + 1] << "\n";
840  std::cout << " \t max: \t" << allBoxes[24 * i + 4 * ( pl - 1 ) + 2] << " \t"
841  << allBoxes[24 * i + 4 * ( pl - 1 ) + 3] << "\n";
842  }
843  }
844  }
845 #endif
846 
847  return MB_SUCCESS;
848 }
849 //#define VERBOSE
850 // this will use the bounding boxes for the (euler)/ fix mesh that are already established
851 // will distribute the mesh to other procs, so that on each task, the covering set covers the local
852 // bounding box this means it will cover the second (local) mesh set; So the covering set will cover
853 // completely the second local mesh set (in intersection)
854 ErrorCode Intx2MeshOnSphere::construct_covering_set( EntityHandle& initial_distributed_set,
855  EntityHandle& covering_set,
856  bool gnomonic,
857  int order )
858 {
859  // primary element came from, in the joint communicator ; this will be forwarded by coverage
860  // mesh needed for tag migrate later on
861  int defaultInt = -1; // no processor, so it was not migrated from somewhere else
862  ErrorCode rval = mb->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag,
863  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" );
864 
865  assert( parcomm != NULL );
866  Range meshCells;
867  rval = mb->get_entities_by_dimension( initial_distributed_set, 2, meshCells );MB_CHK_SET_ERR( rval, "can't get cells by dimension from mesh set" );
868 
869  bool extraWork = (order >= 2);
870  if( 1 == parcomm->proc_config().proc_size() )
871  {
872  // move all initial cells to coverage set
873  rval = mb->add_entities( covering_set, meshCells );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
874  // if point cloud source, add vertices
875  if( 0 == meshCells.size() || max_edges_1 == 0 )
876  {
877  // add vertices from the source set
878  Range verts;
879  rval = mb->get_entities_by_dimension( initial_distributed_set, 0, verts );MB_CHK_SET_ERR( rval, "can't get vertices from mesh set" );
880  rval = mb->add_entities( covering_set, verts );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
881  }
882  return MB_SUCCESS;
883  }
884 
885  // mark on the coverage mesh where this element came from
886  Tag sendProcTag; /// for coverage mesh, will store the sender
887  rval = mb->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag, MB_TAG_DENSE | MB_TAG_CREAT,
888  &defaultInt );MB_CHK_SET_ERR( rval, "can't create sending processor tag" );
889 
890  // this information needs to be forwarded to coverage mesh, if this mesh was already migrated
891  // from somewhere else
892  // look at the value of orgSendProcTag for one mesh cell; if -1, no need to forward that; if
893  // !=-1, we know that this mesh was migrated, we need to find out more about origin of cell
894  int orig_sender = -1;
895  EntityHandle oneCell = 0;
896  // decide if we need to transfer global DOFs info attached to each HOMME coarse cell; first we
897  // need to decide if the mesh has that tag; will affect the size of the tuple list involved in
898  // the crystal routing
899  int size_gdofs_tag = 0;
900  std::vector< int > valsDOFs;
901  Tag gdsTag;
902  rval = mb->tag_get_handle( "GLOBAL_DOFS", gdsTag );
903 
904  if( meshCells.size() > 0 )
905  {
906  oneCell = meshCells[0]; // it is possible we do not have any cells, even after migration
907  rval = mb->tag_get_data( orgSendProcTag, &oneCell, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sending processor value" );
908  if( gdsTag )
909  {
910  DataType dtype;
911  rval = mb->tag_get_data_type( gdsTag, dtype );
912  if( MB_SUCCESS == rval && MB_TYPE_INTEGER == dtype )
913  {
914  // find the values on first cell
915  int lenTag = 0;
916  rval = mb->tag_get_length( gdsTag, lenTag );
917  if( MB_SUCCESS == rval && lenTag > 0 )
918  {
919  valsDOFs.resize( lenTag );
920  rval = mb->tag_get_data( gdsTag, &oneCell, 1, &valsDOFs[0] );
921  if( MB_SUCCESS == rval && valsDOFs[0] > 0 )
922  {
923  // first value positive means we really need to transport this data during
924  // coverage
925  size_gdofs_tag = lenTag;
926  }
927  }
928  }
929  }
930  }
931 
932  // another collective call, to see if the mesh is migrated and if the GLOBAL_DOFS tag need to be
933  // transferred over to the coverage mesh. It is possible that there is no initial mesh source
934  // mesh on the task, so we do not know that info from the tag but TupleList needs to be sized
935  // uniformly for all tasks. Do a collective MPI_MAX to see if it is migrated and if we have
936  // (collectively) a GLOBAL_DOFS task
937 
938  int local_int_array[2], global_int_array[2];
939  local_int_array[0] = orig_sender;
940  local_int_array[1] = size_gdofs_tag;
941  // now reduce over all processors
942  int mpi_err =
943  MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
944  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
945  orig_sender = global_int_array[0];
946  size_gdofs_tag = global_int_array[1];
947 #ifdef VERBOSE
948  std::cout << "proc: " << my_rank << " size_gdofs_tag:" << size_gdofs_tag << "\n";
949 #endif
950  valsDOFs.resize( size_gdofs_tag );
951 
952  // finally, we have correct global info, we can decide if mesh was migrated and if we have
953  // global dofs tag that need to be sent with coverage info
954  int migrated_mesh = 0;
955  if( orig_sender != -1 ) migrated_mesh = 1; //
956 
957  int ghost_info = 0;
958  if (extraWork)
959  ghost_info = 1; // an extra field for ghost cell owner; we need it
960  // if size_gdofs_tag>0, we are sure valsDOFs got resized to what we need
961 
962  // get all mesh verts1
963  Range mesh_verts;
964  rval = mb->get_connectivity( meshCells, mesh_verts );MB_CHK_SET_ERR( rval, "can't get mesh vertices" );
965  size_t num_mesh_verts = mesh_verts.size();
966 
967  // now see the mesh points positions; to what boxes should we send them?
968  std::vector< double > coords_mesh( 3 * num_mesh_verts );
969  rval = mb->get_coords( mesh_verts, &coords_mesh[0] );MB_CHK_SET_ERR( rval, "can't get mesh points position" );
970 
971  // decide gnomonic plane for each vertex, as in the compute boxes
972  std::vector< int > gnplane;
973  if( gnomonic )
974  {
975  gnplane.resize( num_mesh_verts );
976  for( size_t i = 0; i < num_mesh_verts; i++ )
977  {
978  CartVect pos( &coords_mesh[3 * i] );
979  int pl;
981  gnplane[i] = pl;
982  }
983  }
984 
985  std::vector< int > gids( num_mesh_verts );
986  rval = mb->tag_get_data( gid, mesh_verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get vertices gids" );
987 
988  // ranges to send to each processor; will hold vertices and elements (quads/ polygons)
989  // will look if the box of the mesh cell covers bounding box(es) (within tolerances)
990  std::map< int, Range > Rto;
991  int numprocs = parcomm->proc_config().proc_size();
992 
993  for( Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
994  {
995  EntityHandle q = *eit;
996  const EntityHandle* conn;
997  int num_nodes;
998  rval = mb->get_connectivity( q, conn, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity on cell" );
999 
1000  // first decide what planes need to consider
1001  std::set< int > planes; // if this list contains more than 3 planes, we have a very bad mesh!!!
1002  std::vector< double > elco( 3 * num_nodes );
1003  for( int i = 0; i < num_nodes; i++ )
1004  {
1005  EntityHandle v = conn[i];
1006  int index = mesh_verts.index( v );
1007  if( gnomonic ) planes.insert( gnplane[index] );
1008  for( int j = 0; j < 3; j++ )
1009  {
1010  elco[3 * i + j] = coords_mesh[3 * index + j]; // extract from coords
1011  }
1012  }
1013  if( gnomonic )
1014  {
1015  // now loop over all planes that need to be considered for this element
1016  for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
1017  {
1018  int pl = *st; // gnomonic plane considered
1019  double qmin[2] = { DBL_MAX, DBL_MAX };
1020  double qmax[2] = { -DBL_MAX, -DBL_MAX };
1021  for( int i = 0; i < num_nodes; i++ )
1022  {
1023  CartVect dp( &elco[3 * i] ); // uses constructor for CartVect that takes a
1024  // pointer to double
1025  // gnomonic projection
1026  double c2[2];
1027  IntxUtils::gnomonic_projection( dp, Rsrc, pl, c2[0], c2[1] ); // 2 coordinates
1028  for( int j = 0; j < 2; j++ )
1029  {
1030  if( qmin[j] > c2[j] ) qmin[j] = c2[j];
1031  if( qmax[j] < c2[j] ) qmax[j] = c2[j];
1032  }
1033  }
1034  // now decide if processor p should be interested in this cell, by looking at plane pl
1035  // 2d box this is one of the few size n loops;
1036  for( int p = 0; p < numprocs; p++ ) // each cell q can be sent to more than one processor
1037  {
1038  double procMin1 = allBoxes[24 * p + 4 * ( pl - 1 )]; // these were determined before
1039  //
1040  if( procMin1 >= DBL_MAX ) // the processor has no targets on this plane
1041  continue;
1042  double procMin2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
1043  double procMax1 = allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
1044  double procMax2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
1045  // test overlap of 2d boxes
1046  if( procMin1 > qmax[0] + box_error || procMin2 > qmax[1] + box_error ) continue; //
1047  if( qmin[0] > procMax1 + box_error || qmin[1] > procMax2 + box_error ) continue;
1048  // good to be inserted
1049  Rto[p].insert( q );
1050  }
1051  }
1052  }
1053  else // regular 3d box; one box per processor
1054  {
1055  for( int p = 0; p < numprocs; p++ )
1056  {
1057  BoundBox box( &allBoxes[6 * p] );
1058  bool insert = false;
1059  for( int i = 0; i < num_nodes; i++ )
1060  {
1061  if( box.contains_point( &elco[3 * i], box_error ) )
1062  {
1063  insert = true;
1064  break;
1065  }
1066  }
1067  if( insert ) Rto[p].insert( q );
1068  }
1069  }
1070  }
1071 
1072  // now, for higher order, we might need to send to processor p, not only the Range Rto[p] cells;
1073  // we need to augment those ranges with adjacent cells, according to the order passed
1074  // if order is 1, not do anything
1075  // if order is 2, for example, we need to add all adj cells level 1, by edge
1076  if (extraWork)
1077  {
1078  for (int p = 0; p < numprocs; p++ )
1079  {
1080  Range originalSend = Rto[p];
1081  // determine all adjacent cells for order 2 and higher
1082  // Need to get layers of bridge-adj entities
1083  if( originalSend.empty() ) continue;
1084  Range extraCells;
1085  rval = MeshTopoUtil( mb ).get_bridge_adjacencies( originalSend, 0, 2, extraCells, order - 1 ); MB_CHK_SET_ERR( rval, "Failed to get bridge adjacencies" );
1086  // big miss : need to merge only cells from initial source (ghost) set;
1087  // get_bridge adj will get all cells adjacent to a vertex
1088  extraCells = intersect(extraCells, meshCells);
1089  Rto[p].merge(extraCells);
1090  }
1091  }
1092 
1093  // here, we will use crystal router to send each cell to designated tasks (mesh migration)
1094 
1095  // a better implementation would be to use pcomm send / recv entities; a good test case
1096  // pcomm send / receives uses point to point communication, not global gather / scatter
1097 
1098  // now, build TLv and TLq (tuple list for vertices and cells, separately sent)
1099  size_t numq = 0;
1100  size_t numv = 0;
1101 
1102  // merge the list of vertices to be sent
1103  for( int p = 0; p < numprocs; p++ )
1104  {
1105  if( p == (int)my_rank ) continue; // do not "send" it to current task, because it is already here
1106  Range& range_to_P = Rto[p];
1107  // add the vertices to it
1108  if( range_to_P.empty() ) continue; // nothing to send to proc p
1109 #ifdef VERBOSE
1110  std::cout << " proc : " << my_rank << " to proc " << p << " send " << range_to_P.size() << " cells "
1111  << " psize: " << range_to_P.psize() << "\n";
1112 #endif
1113  Range vertsToP;
1114  rval = mb->get_connectivity( range_to_P, vertsToP );MB_CHK_SET_ERR( rval, "can't get connectivity" );
1115  numq = numq + range_to_P.size();
1116  numv = numv + vertsToP.size();
1117  range_to_P.merge( vertsToP );
1118  }
1119 
1120  TupleList TLv; // send vertices with a different tuple list
1121  TupleList TLq;
1122  TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, 3 real coordinates
1123  TLv.enableWriteAccess();
1124 
1125  // add also GLOBAL_DOFS info, if found on the mesh cell; it should be found only on HOMME cells!
1126  int sizeTuple =
1127  2 + max_edges_1 + migrated_mesh + size_gdofs_tag + ghost_info; // max edges could be up to MAXEDGES :) for polygons
1128  TLq.initialize( sizeTuple, 0, 0, 0,
1129  numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v), plus
1130  // original sender if set (migrated mesh case)
1131  // we will not send the entity handle, global ID should be more than enough
1132  // we will not need more than 2B vertices TODO 2B vertices or cells
1133  // if we need more than 2B, we will need to use a different marker anyway (GLOBAL ID is not
1134  // enough then)
1135 
1136  TLq.enableWriteAccess();
1137 #ifdef VERBOSE
1138  std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1139 #endif
1140 
1141  for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1142  {
1143  if( to_proc == (int)my_rank ) continue;
1144  Range& range_to_P = Rto[to_proc];
1145  Range V = range_to_P.subset_by_type( MBVERTEX );
1146 
1147  for( Range::iterator it = V.begin(); it != V.end(); ++it )
1148  {
1149  EntityHandle v = *it;
1150  int index = mesh_verts.index( v ); //
1151  assert( -1 != index );
1152  int n = TLv.get_n(); // current size of tuple list
1153  TLv.vi_wr[2 * n] = to_proc; // send to processor
1154  TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the second_mesh_verts range
1155  TLv.vr_wr[3 * n] = coords_mesh[3 * index]; // departure position, of the node local_verts[i]
1156  TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
1157  TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
1158  TLv.inc_n(); // increment tuple list size
1159  }
1160  // also, prep the 2d cells for sending ...
1161  Range Q = range_to_P.subset_by_dimension( 2 );
1162  for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
1163  {
1164  EntityHandle q = *it; // this is a second mesh cell (or src, lagrange set)
1165  int global_id;
1166  rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_SET_ERR( rval, "can't get gid for polygon" );
1167  int n = TLq.get_n(); // current size
1168  TLq.vi_wr[sizeTuple * n] = to_proc; //
1169  TLq.vi_wr[sizeTuple * n + 1] =
1170  global_id; // global id of element, used to identify it for debug purposes only
1171  const EntityHandle* conn4;
1172  int num_nodes; // could be up to MAXEDGES; max_edges?;
1173  rval = mb->get_connectivity( q, conn4, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity for cell" );
1174  if( num_nodes > max_edges_1 )
1175  {
1176  mb->list_entities( &q, 1 );MB_CHK_SET_ERR( MB_FAILURE, "too many nodes in a cell (" << num_nodes << "," << max_edges_1 << ")" );
1177  }
1178  for( int i = 0; i < num_nodes; i++ )
1179  {
1180  EntityHandle v = conn4[i];
1181  int index = mesh_verts.index( v );
1182  assert( -1 != index );
1183  TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1184  }
1185  for( int k = num_nodes; k < max_edges_1; k++ )
1186  {
1187  TLq.vi_wr[sizeTuple * n + 2 + k] =
1188  0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1189  }
1190  int currentIndexIntTuple = 2 + max_edges_1;
1191  // is the mesh migrated before or not?
1192  if( migrated_mesh )
1193  {
1194  // case of extra work, maybe need to check if it is ghost ? yes, next loop !
1195  rval = mb->tag_get_data( orgSendProcTag, &q, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sender for polygon, in migrate scenario" );
1196  TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender; // should be different than -1
1197  currentIndexIntTuple++;
1198  }
1199  if (ghost_info > 0) // extraWork
1200  {
1201  // case of ghost
1202  int owner = my_rank;
1203  // this could happen if extra work, real owner is different ?
1204  rval = parcomm->get_owner(q, owner); MB_CHK_SET_ERR( rval, "can't get owner for cell" );
1205  TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = owner; // should be different than -1
1206  currentIndexIntTuple++;
1207  }
1208  // GLOBAL_DOFS info, if available
1209  if( size_gdofs_tag )
1210  {
1211  rval = mb->tag_get_data( gdsTag, &q, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't get gdofs data in HOMME" );
1212  for( int i = 0; i < size_gdofs_tag; i++ )
1213  {
1214  TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
1215  valsDOFs[i]; // should be different than 0 or -1
1216  }
1217  }
1218 
1219  TLq.inc_n(); // increment tuple list size
1220  }
1221  } // end for loop over total number of processors
1222 
1223  // now we are done populating the tuples; route them to the appropriate processors
1224  // this does the communication magic
1225  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1226  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1227 
1228  // the first mesh elements are in localEnts; we do not need them at all
1229 
1230  // maps from global ids to new vertex and cell handles, that are added
1231 
1232  std::map< int, EntityHandle > globalID_to_vertex_handle;
1233  // we already have some vertices from second mesh set; they are already in the processor, even
1234  // before receiving other verts from neighbors this is an inverse map from gid to vertex handle,
1235  // which is local here, we do not want to duplicate vertices their identifier is the global ID!!
1236  // it must be unique per mesh ! (I mean, first mesh, source); gid for second mesh is not needed here
1237  int k = 0;
1238  for( Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
1239  {
1240  globalID_to_vertex_handle[gids[k]] = *vit;
1241  }
1242  /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one?
1243  globalID_to_eh.clear(); // we need it now in case of extra work, to not duplicate cells
1244 
1245  // now, look at every TLv, and see if we have to create a vertex there or not
1246  int n = TLv.get_n(); // the size of the points received
1247  for( int i = 0; i < n; i++ )
1248  {
1249  int globalId = TLv.vi_rd[2 * i + 1];
1250  if( globalID_to_vertex_handle.find( globalId ) ==
1251  globalID_to_vertex_handle.end() ) // we do not have locally this vertex (yet)
1252  // so we have to create it, and add to the inverse map
1253  {
1254  EntityHandle new_vert;
1255  double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1256  rval = mb->create_vertex( dp_pos, new_vert );MB_CHK_SET_ERR( rval, "can't create new vertex " );
1257  globalID_to_vertex_handle[globalId] = new_vert; // now add it to the map
1258  // set the GLOBAL ID tag on the new vertex
1259  rval = mb->tag_set_data( gid, &new_vert, 1, &globalId );MB_CHK_SET_ERR( rval, "can't set global ID tag on new vertex " );
1260  }
1261  }
1262 
1263  // now, all necessary vertices should be created
1264  // look in the local list of 2d cells for this proc, and add all those cells to covering set
1265  // also
1266 
1267  Range& local = Rto[my_rank];
1268  Range local_q = local.subset_by_dimension( 2 );
1269 
1270  for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1271  {
1272  EntityHandle q = *it; // these are from source cells, local
1273  int gid_el;
1274  rval = mb->tag_get_data( gid, &q, 1, &gid_el );MB_CHK_SET_ERR( rval, "can't get global id of cell " );
1275  assert( gid_el >= 0 );
1276  globalID_to_eh[gid_el] = q; // do we need this? yes, now we do; parent tags are now using it heavily
1277  if (extraWork)
1278  {
1279  int owner = my_rank;
1280  // this could happen if extra work, real owner is different ?
1281  rval = parcomm->get_owner(q, owner); MB_CHK_SET_ERR( rval, "can't get owner for cell" );
1282  rval = mb->tag_set_data( sendProcTag, &q, 1, &owner );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
1283  }
1284  else
1285  {
1286  rval = mb->tag_set_data( sendProcTag, &q, 1, &my_rank );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
1287  }
1288  }
1289 
1290  // now look at all elements received through; we do not want to duplicate them
1291  n = TLq.get_n(); // number of elements received by this processor
1292  // a cell should be received from one proc only; so why are we so worried about duplicated
1293  // elements? a vertex can be received from multiple sources, that is fine
1294 
1295  for( int i = 0; i < n; i++ )
1296  {
1297  int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1298  // int from_proc=TLq.vi_rd[sizeTuple * i ]; // we do not need from_proc anymore
1299 
1300  // do we already have a cell with this global ID, represented?
1301  // yes, it could happen for extraWork !
1302  if (globalID_to_eh.find(globalIdEl) != globalID_to_eh.end())
1303  continue;
1304  // construct the conn triangle , quad or polygon
1305  EntityHandle new_conn[MAXEDGES]; // we should use std::vector with max_edges_1
1306  int nnodes = -1;
1307  for( int j = 0; j < max_edges_1; j++ )
1308  {
1309  int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1310  if( vgid == 0 )
1311  new_conn[j] = 0; // this can actually happen for polygon mesh (when we have less
1312  // number of vertices than max_edges)
1313  else
1314  {
1315  assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
1316  new_conn[j] = globalID_to_vertex_handle[vgid];
1317  nnodes = j + 1; // nodes are at the beginning, and are variable number
1318  }
1319  }
1320  EntityHandle new_element;
1321  //
1322  EntityType entType = MBQUAD;
1323  if( nnodes > 4 ) entType = MBPOLYGON;
1324  if( nnodes < 4 ) entType = MBTRI;
1325  rval = mb->create_element( entType, new_conn, nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element for second mesh " );
1326 
1327  globalID_to_eh[globalIdEl] = new_element;
1328  local_q.insert( new_element );
1329  rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set gid for cell " );
1330  int currentIndexIntTuple = 2 + max_edges_1;
1331  if( migrated_mesh )
1332  {
1333  orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1334  rval = mb->tag_set_data( orgSendProcTag, &new_element, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't set original sender for cell, in migrate scenario" );
1335  currentIndexIntTuple++; // add one more
1336  }
1337  // store also the processor this coverage element came from
1338  int from_proc = TLq.vi_rd[sizeTuple * i];
1339  if( ghost_info ) // ghost info will have the original owner of the coverage cell
1340  {
1341  // case of ghost
1342  from_proc = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1343  currentIndexIntTuple++; // add one more
1344  }
1345  rval = mb->tag_set_data( sendProcTag, &new_element, 1, &from_proc );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
1346 
1347  // check if we need to retrieve and set GLOBAL_DOFS data
1348  if( size_gdofs_tag )
1349  {
1350  for( int j = 0; j < size_gdofs_tag; j++ )
1351  {
1352  valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
1353  }
1354  rval = mb->tag_set_data( gdsTag, &new_element, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't set GLOBAL_DOFS data on coverage mesh" );
1355  }
1356 
1357  }
1358 
1359  // now, add to the covering_set the elements created in the local_q range
1360  rval = mb->add_entities( covering_set, local_q );MB_CHK_SET_ERR( rval, "can't add entities to new mesh set " );
1361 #ifdef VERBOSE
1362  std::cout << " proc " << my_rank << " add " << local_q.size() << " cells to covering set \n";
1363 #endif
1364  return MB_SUCCESS;
1365 }
1366 
1367 #endif // MOAB_HAVE_MPI
1368 //#undef VERBOSE
1369 #undef CHECK_CONVEXITY
1370 
1371 } /* namespace moab */