Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  11 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp" 12 #include "moab/IntxMesh/IntxUtils.hpp" 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  28 Intx2MeshOnSphere::Intx2MeshOnSphere( Interface* mbimpl, IntxAreaUtils::AreaMethod amethod ) 29  : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 ) 30 { 31 } 32  33 Intx2MeshOnSphere::~Intx2MeshOnSphere() {} 34  35 /* 36  * return also the area for robustness verification 37  */ 38 double Intx2MeshOnSphere::setup_tgt_cell( EntityHandle tgt, int& nsTgt ) 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  * */ 78 ErrorCode Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt, 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  { 122  rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], 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 = 477  areaAdaptor.area_spherical_triangle( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest ); 478  if( orientedArea < 0 ) 479  { 480  std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n"; 481  for( int i = 0; i < nP; i++ ) 482  { 483  int nexti = ( i + 1 ) % nP; 484  double lengthEdge = ( posi[i] - posi[nexti] ).length(); 485  std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n"; 486  } 487  std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n"; 488  489  std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k 490  << " target, src:" << tgt << " " << src << " \n"; 491  } 492  } 493 #endif 494  495 #ifdef ENABLE_DEBUG 496  if( dbg_1 ) 497  { 498  std::cout << "Counting: " << counting << "\n"; 499  std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :"; 500  for( int i1 = 0; i1 < nP; i1++ ) 501  std::cout << " " << mb->id_from_handle( foundIds[i1] ); 502  std::cout << " plane: " << plane << "\n"; 503  std::vector< CartVect > posi( nP ); 504  mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) ); 505  for( int i1 = 0; i1 < nP; i1++ ) 506  std::cout << foundIds[i1] << " " << posi[i1] << "\n"; 507  508  std::stringstream fff; 509  fff << "file0" << counting << ".vtk"; 510  rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval ); 511  } 512 #endif 513  } 514  // else { 515  // std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n"; 516  // } 517  // disable_debug(); 518  return MB_SUCCESS; 519 } 520  521 ErrorCode Intx2MeshOnSphere::update_tracer_data( EntityHandle out_set, Tag& tagElem, Tag& tagArea ) 522 { 523  EntityHandle dum = 0; 524  525  Tag corrTag; 526  ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, 527  &dum ); // it should have been created 528  MB_CHK_SET_ERR( rval, "can't get correlation tag" ); 529  530  // get all polygons out of out_set; then see where are they coming from 531  Range polys; 532  rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" ); 533  534  // rs2 is the target range, arrival; rs1 is src, departure; 535  // there is a connection between rs1 and rs2, through the corrTag 536  // corrTag is __correlation 537  // basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly); 538  // also, mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly); 539  // we start from rs2 existing, then we have to update something 540  541  // tagElem will have multiple tracers 542  int numTracers = 0; 543  rval = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" ); 544  if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" ); 545  546  std::vector< double > currentVals( rs2.size() * numTracers ); 547  rval = mb->tag_get_data( tagElem, rs2, &currentVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" ); 548  549  // create new tuple list for tracers to other processors, from remote_cells 550 #ifdef MOAB_HAVE_MPI 551  if( remote_cells ) 552  { 553  int n = remote_cells->get_n(); 554  if( n > 0 ) 555  { 556  remote_cells_with_tracers = new TupleList(); 557  remote_cells_with_tracers->initialize( 2, 0, 1, numTracers, 558  n ); // tracers are in these tuples 559  remote_cells_with_tracers->enableWriteAccess(); 560  for( int i = 0; i < n; i++ ) 561  { 562  remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i]; 563  remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1]; 564  // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication 565  remote_cells_with_tracers->vul_wr[i] = 566  remote_cells->vul_wr[i]; // this is the corresponding target cell (arrival) 567  for( int k = 0; k < numTracers; k++ ) 568  remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0; // initialize tracers to be transported 569  remote_cells_with_tracers->inc_n(); 570  } 571  } 572  delete remote_cells; 573  remote_cells = NULL; 574  } 575 #endif 576  // for each polygon, we have 2 indices: target and source parents 577  // we need index source to update index tgt? 578  std::vector< double > newValues( rs2.size() * numTracers, 579  0. ); // initialize with 0 all of them 580  // area of the polygon * conc on target (old) current quantity 581  // finally, divide by the area of the tgt 582  double check_intx_area = 0.; 583  moab::IntxAreaUtils intxAreas( this->areaMethod ); // use_lHuiller = true 584  for( Range::iterator it = polys.begin(); it != polys.end(); ++it ) 585  { 586  EntityHandle poly = *it; 587  int srcIndex, tgtIndex; 588  rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" ); 589  590  EntityHandle src = rs1[srcIndex - 1]; // big assumption, it should work for meshes where global id is the same 591  // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes 592  rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" ); 593  // EntityHandle target = rs2[tgtIndex]; 594  // big assumption here, target and source are "parallel" ;we should have an index from 595  // source to target (so a deformed source corresponds to an arrival tgt) 596  /// TODO: VSM: Its unclear whether we need the source or destination radius here. 597  double radius = Rsrc; 598  double areap = intxAreas.area_spherical_element( mb, poly, radius ); 599  check_intx_area += areap; 600  // so the departure cell at time t (srcIndex) covers a portion of a tgtCell 601  // that quantity will be transported to the tgtCell at time t+dt 602  // the source corresponds to a target arrival 603  EntityHandle tgtArr; 604  rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr ); 605  if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval ) 606  { 607 #ifdef MOAB_HAVE_MPI 608  if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" ); 609  // maybe the element is remote, from another processor 610  int global_id_src; 611  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" ); 612  // find the 613  int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src ); 614  if( index_in_remote == -1 ) 615  MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" ); 616  for( int k = 0; k < numTracers; k++ ) 617  remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] += 618  currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap; 619 #endif 620  } 621  else if( MB_SUCCESS == rval ) 622  { 623  int arrTgtIndex = rs2.index( tgtArr ); 624  if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" ); 625  for( int k = 0; k < numTracers; k++ ) 626  newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap; 627  } 628  629  else 630  MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " ); 631  } 632  // now, send back the remote_cells_with_tracers to the processors they came from, with the 633  // updated values for the tracer mass in a cell 634 #ifdef MOAB_HAVE_MPI 635  if( remote_cells_with_tracers ) 636  { 637  // so this means that some cells will be sent back with tracer info to the procs they were 638  // sent from 639  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 ); 640  // now, look at the global id, find the proper "tgt" cell with that index and update its 641  // mass 642  // remote_cells->print("remote cells after routing"); 643  int n = remote_cells_with_tracers->get_n(); 644  for( int j = 0; j < n; j++ ) 645  { 646  EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j]; // entity handle sent back 647  int arrTgtIndex = rs2.index( tgtCell ); 648  if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" ); 649  for( int k = 0; k < numTracers; k++ ) 650  newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k]; 651  } 652  } 653 #endif /* MOAB_HAVE_MPI */ 654  // now divide by target area (current) 655  int j = 0; 656  Range::iterator iter = rs2.begin(); 657  void* data = NULL; // used for stored area 658  int count = 0; 659  std::vector< double > total_mass_local( numTracers, 0. ); 660  while( iter != rs2.end() ) 661  { 662  rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" ); 663  double* ptrArea = (double*)data; 664  for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ ) 665  { 666  for( int k = 0; k < numTracers; k++ ) 667  { 668  total_mass_local[k] += newValues[j * numTracers + k]; 669  newValues[j * numTracers + k] /= ( *ptrArea ); 670  } 671  } 672  } 673  rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" ); 674  675 #ifdef MOAB_HAVE_MPI 676  std::vector< double > total_mass( numTracers, 0. ); 677  double total_intx_area = 0; 678  int mpi_err = 679  MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); 680  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 681  // now reduce total area 682  mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); 683  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 684  if( my_rank == 0 ) 685  { 686  for( int k = 0; k < numTracers; k++ ) 687  std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n"; 688  std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " " 689  << total_intx_area << "\n"; 690  } 691  692  if( remote_cells_with_tracers ) 693  { 694  delete remote_cells_with_tracers; 695  remote_cells_with_tracers = NULL; 696  } 697 #else 698  for( int k = 0; k < numTracers; k++ ) 699  std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n"; 700  std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " " 701  << check_intx_area << "\n"; 702 #endif 703  return MB_SUCCESS; 704 } 705  706 #ifdef MOAB_HAVE_MPI 707 ErrorCode Intx2MeshOnSphere::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic ) 708 { 709  if( !gnomonic ) 710  { 711  return Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts, gnomonic ); 712  } 713  // so here, we know that the logic is for gnomonic == true 714  localEnts.clear(); 715  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );MB_CHK_SET_ERR( rval, "can't get local ents" ); 716  717  rval = mb->get_connectivity( localEnts, local_verts );MB_CHK_SET_ERR( rval, "can't get connectivity" ); 718  int num_local_verts = (int)local_verts.size(); 719  720  assert( parcomm != NULL ); 721  722  if( num_local_verts == 0 ) 723  { 724  // it is probably point cloud, get the local vertices from set 725  rval = mb->get_entities_by_dimension( euler_set, 0, local_verts );MB_CHK_SET_ERR( rval, "can't get local vertices from set" ); 726  num_local_verts = (int)local_verts.size(); 727  localEnts = local_verts; 728  } 729  // will use 6 gnomonic planes to decide boxes for each gnomonic plane 730  // each gnomonic box will be 2d, min, max 731  double gnom_box[24]; 732  for( int i = 0; i < 6; i++ ) 733  { 734  gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX; 735  gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX; 736  } 737  738  // there are 6 gnomonic planes; some elements could be on the corners, and affect multiple 739  // planes decide what gnomonic planes will be affected by each cell some elements could appear 740  // in multiple gnomonic planes ! 741  std::vector< double > coords( 3 * num_local_verts ); 742  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 " ); 743  // decide each local vertex to what gnomonic plane it belongs 744  745  std::vector< int > gnplane; 746  gnplane.resize( num_local_verts ); 747  for( int i = 0; i < num_local_verts; i++ ) 748  { 749  CartVect pos( &coords[3 * i] ); 750  int pl; 751  IntxUtils::decide_gnomonic_plane( pos, pl ); 752  gnplane[i] = pl; 753  } 754  755  for( Range::iterator it = localEnts.begin(); it != localEnts.end(); it++ ) 756  { 757  EntityHandle cell = *it; 758  EntityType typeCell = mb->type_from_handle( cell ); // could be vertex, for point cloud 759  // get coordinates, and decide gnomonic planes for it 760  int nnodes; 761  const EntityHandle* conn = NULL; 762  EntityHandle c[1]; 763  if( typeCell != MBVERTEX ) 764  { 765  rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" ); 766  } 767  else 768  { 769  nnodes = 1; 770  c[0] = cell; // actual node 771  conn = &c[0]; 772  } 773  // get coordinates of vertices involved with this 774  std::vector< double > elco( 3 * nnodes ); 775  std::set< int > planes; 776  for( int i = 0; i < nnodes; i++ ) 777  { 778  int ix = local_verts.index( conn[i] ); 779  planes.insert( gnplane[ix] ); 780  for( int j = 0; j < 3; j++ ) 781  { 782  elco[3 * i + j] = coords[3 * ix + j]; 783  } 784  } 785  // now, augment the boxes for all planes involved 786  for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ ) 787  { 788  int pl = *st; 789  for( int i = 0; i < nnodes; i++ ) 790  { 791  CartVect pos( &elco[3 * i] ); 792  double c2[2]; 793  IntxUtils::gnomonic_projection( pos, Rdest, pl, c2[0], 794  c2[1] ); // 2 coordinates 795  // 796  for( int k = 0; k < 2; k++ ) 797  { 798  double val = c2[k]; 799  if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val; // min in k direction 800  if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] ) 801  gnom_box[4 * ( pl - 1 ) + 2 + k] = val; // max in k direction 802  } 803  } 804  } 805  } 806  807  int numprocs = parcomm->proc_config().proc_size(); 808  allBoxes.resize( 24 * numprocs ); // 6 gnomonic planes , 4 doubles for each for 2d box 809  810  my_rank = parcomm->proc_config().proc_rank(); 811  for( int k = 0; k < 24; k++ ) 812  allBoxes[24 * my_rank + k] = gnom_box[k]; 813  814  // now communicate to get all boxes 815  int mpi_err; 816 #if( MPI_VERSION >= 2 ) 817  // use "in place" option 818  mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 24, MPI_DOUBLE, 819  parcomm->proc_config().proc_comm() ); 820 #else 821  { 822  std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() ); 823  mpi_err = MPI_Allgather( &allBoxes[24 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE, 824  parcomm->proc_config().proc_comm() ); 825  allBoxes = allBoxes_tmp; 826  } 827 #endif 828  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 829  830 #ifdef VERBOSE 831  if( my_rank == 0 ) 832  { 833  std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2 834  << " on second mesh \n"; 835  for( int i = 0; i < numprocs; i++ ) 836  { 837  std::cout << "task: " << i << " \n"; 838  for( int pl = 1; pl <= 6; pl++ ) 839  { 840  std::cout << " plane " << pl << " min: \t" << allBoxes[24 * i + 4 * ( pl - 1 )] << " \t" 841  << allBoxes[24 * i + 4 * ( pl - 1 ) + 1] << "\n"; 842  std::cout << " \t max: \t" << allBoxes[24 * i + 4 * ( pl - 1 ) + 2] << " \t" 843  << allBoxes[24 * i + 4 * ( pl - 1 ) + 3] << "\n"; 844  } 845  } 846  } 847 #endif 848  849  return MB_SUCCESS; 850 } 851 //#define VERBOSE 852 // this will use the bounding boxes for the (euler)/ fix mesh that are already established 853 // will distribute the mesh to other procs, so that on each task, the covering set covers the local 854 // bounding box this means it will cover the second (local) mesh set; So the covering set will cover 855 // completely the second local mesh set (in intersection) 856 // now, when covering set needs to have extra layers, we will increase dramatically the box_eps, from something close to 0, 857 // to something larger than the "source mesh size" * sqrt(3) for each layer needed 858 // so the first step is finding the global diagonal mesh size in the source mesh 859 ErrorCode Intx2MeshOnSphere::construct_covering_set( EntityHandle& initial_distributed_set, 860  EntityHandle& covering_set, 861  bool gnomonic, 862  int nb_ghost_layers ) 863 { 864  // primary element came from, in the joint communicator ; this will be forwarded by coverage 865  // mesh needed for tag migrate later on 866  int defaultInt = -1; // no processor, so it was not migrated from somewhere else 867  ErrorCode rval = mb->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag, 868  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" ); 869  870  assert( parcomm != NULL ); 871  Range meshCells; 872  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" ); 873  874  if( 1 == parcomm->proc_config().proc_size() ) 875  { 876  // move all initial cells to coverage set 877  rval = mb->add_entities( covering_set, meshCells );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" ); 878  // if point cloud source, add vertices 879  if( 0 == meshCells.size() || max_edges_1 == 0 ) 880  { 881  // add vertices from the source set 882  Range verts; 883  rval = mb->get_entities_by_dimension( initial_distributed_set, 0, verts );MB_CHK_SET_ERR( rval, "can't get vertices from mesh set" ); 884  rval = mb->add_entities( covering_set, verts );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" ); 885  } 886  return MB_SUCCESS; 887  } 888  889  // mark on the coverage mesh where this element came from 890  Tag sendProcTag; /// for coverage mesh, will store the sender 891  rval = mb->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag, MB_TAG_DENSE | MB_TAG_CREAT, 892  &defaultInt );MB_CHK_SET_ERR( rval, "can't create sending processor tag" ); 893  894  // this information needs to be forwarded to coverage mesh, if this mesh was already migrated 895  // from somewhere else 896  // look at the value of orgSendProcTag for one mesh cell; if -1, no need to forward that; if 897  // !=-1, we know that this mesh was migrated, we need to find out more about origin of cell 898  int orig_sender = -1; 899  EntityHandle oneCell = 0; 900  // decide if we need to transfer global DOFs info attached to each HOMME coarse cell; first we 901  // need to decide if the mesh has that tag; will affect the size of the tuple list involved in 902  // the crystal routing 903  int size_gdofs_tag = 0; 904  std::vector< int > valsDOFs; 905  Tag gdsTag; 906  rval = mb->tag_get_handle( "GLOBAL_DOFS", gdsTag ); 907  908  if( meshCells.size() > 0 ) 909  { 910  oneCell = meshCells[0]; // it is possible we do not have any cells, even after migration 911  rval = mb->tag_get_data( orgSendProcTag, &oneCell, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sending processor value" ); 912  if( gdsTag ) 913  { 914  DataType dtype; 915  rval = mb->tag_get_data_type( gdsTag, dtype ); 916  if( MB_SUCCESS == rval && MB_TYPE_INTEGER == dtype ) 917  { 918  // find the values on first cell 919  int lenTag = 0; 920  rval = mb->tag_get_length( gdsTag, lenTag ); 921  if( MB_SUCCESS == rval && lenTag > 0 ) 922  { 923  valsDOFs.resize( lenTag ); 924  rval = mb->tag_get_data( gdsTag, &oneCell, 1, &valsDOFs[0] ); 925  if( MB_SUCCESS == rval && valsDOFs[0] > 0 ) 926  { 927  // first value positive means we really need to transport this data during 928  // coverage 929  size_gdofs_tag = lenTag; 930  } 931  } 932  } 933  } 934  } 935  936  // another collective call, to see if the mesh is migrated and if the GLOBAL_DOFS tag need to be 937  // transferred over to the coverage mesh. It is possible that there is no initial mesh source 938  // mesh on the task, so we do not know that info from the tag but TupleList needs to be sized 939  // uniformly for all tasks. Do a collective MPI_MAX to see if it is migrated and if we have 940  // (collectively) a GLOBAL_DOFS task 941  942  int local_int_array[2], global_int_array[2]; 943  local_int_array[0] = orig_sender; 944  local_int_array[1] = size_gdofs_tag; 945  // now reduce over all processors 946  int mpi_err = 947  MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() ); 948  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 949  orig_sender = global_int_array[0]; 950  size_gdofs_tag = global_int_array[1]; 951 #ifdef VERBOSE 952  std::cout << "proc: " << my_rank << " size_gdofs_tag:" << size_gdofs_tag << "\n"; 953 #endif 954  valsDOFs.resize( size_gdofs_tag ); 955  956  // finally, we have correct global info, we can decide if mesh was migrated and if we have 957  // global dofs tag that need to be sent with coverage info 958  int migrated_mesh = 0; 959  if( orig_sender != -1 ) migrated_mesh = 1; // 960  961  // if size_gdofs_tag>0, we are sure valsDOFs got resized to what we need 962  963  // get all mesh verts1 964  Range mesh_verts; 965  rval = mb->get_connectivity( meshCells, mesh_verts );MB_CHK_SET_ERR( rval, "can't get mesh vertices" ); 966  size_t num_mesh_verts = mesh_verts.size(); 967  968  // now see the mesh points positions; to what boxes should we send them? 969  std::vector< double > coords_mesh( 3 * num_mesh_verts ); 970  rval = mb->get_coords( mesh_verts, &coords_mesh[0] );MB_CHK_SET_ERR( rval, "can't get mesh points position" ); 971  972  // decide gnomonic plane for each vertex, as in the compute boxes 973  std::vector< int > gnplane; 974  if( gnomonic ) 975  { 976  gnplane.resize( num_mesh_verts ); 977  for( size_t i = 0; i < num_mesh_verts; i++ ) 978  { 979  CartVect pos( &coords_mesh[3 * i] ); 980  int pl; 981  IntxUtils::decide_gnomonic_plane( pos, pl ); 982  gnplane[i] = pl; 983  } 984  } 985  986  std::vector< int > gids( num_mesh_verts ); 987  rval = mb->tag_get_data( gid, mesh_verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get vertices gids" ); 988  989  // ranges to send to each processor; will hold vertices and elements (quads/ polygons) 990  // will look if the box of the mesh cell covers bounding box(es) (within tolerances) 991  std::map< int, Range > Rto; 992  int numprocs = parcomm->proc_config().proc_size(); 993  994  // now, box error is pretty small, in general 995  // for bilinear mesh, we need an extra layer, which we will get by increasing the epsilon to catch the extra layer 996  // it will depend on the size of the source mesh 997  // so we will compute the max diagonal length for each cell, on the sphere, so we will modify box_error 998  if( nb_ghost_layers > 0 ) 999  { 1000  double diagonal; 1001  rval = IntxUtils::max_diagonal( mb, meshCells, max_edges_1, diagonal );MB_CHK_SET_ERR( rval, "can't get max diagonal" ); 1002  // 1003  double global_diag = 0; 1004  mpi_err = MPI_Allreduce( &diagonal, &global_diag, 1, MPI_DOUBLE, MPI_MAX, parcomm->proc_config().proc_comm() ); 1005  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 1006  double extra_thickness = global_diag * nb_ghost_layers; 1007  if( gnomonic ) extra_thickness *= sqrt( 3. ); 1008  box_error += extra_thickness; // 1009  if( !my_rank ) 1010  std::cout << "ghost_layers:" << nb_ghost_layers << " max diagonal:" << global_diag 1011  << " extra thickness:" << extra_thickness << " box_error:" << box_error << "\n"; 1012  } 1013  for( Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit ) 1014  { 1015  EntityHandle q = *eit; 1016  const EntityHandle* conn; 1017  int num_nodes; 1018  rval = mb->get_connectivity( q, conn, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity on cell" ); 1019  1020  // first decide what planes need to consider 1021  std::set< int > planes; // if this list contains more than 3 planes, we have a very bad mesh!!! 1022  std::vector< double > elco( 3 * num_nodes ); 1023  for( int i = 0; i < num_nodes; i++ ) 1024  { 1025  EntityHandle v = conn[i]; 1026  int index = mesh_verts.index( v ); 1027  if( gnomonic ) planes.insert( gnplane[index] ); 1028  for( int j = 0; j < 3; j++ ) 1029  { 1030  elco[3 * i + j] = coords_mesh[3 * index + j]; // extract from coords 1031  } 1032  } 1033  if( gnomonic ) 1034  { 1035  // now loop over all planes that need to be considered for this element 1036  for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ ) 1037  { 1038  int pl = *st; // gnomonic plane considered 1039  double qmin[2] = { DBL_MAX, DBL_MAX }; 1040  double qmax[2] = { -DBL_MAX, -DBL_MAX }; 1041  for( int i = 0; i < num_nodes; i++ ) 1042  { 1043  CartVect dp( &elco[3 * i] ); // uses constructor for CartVect that takes a 1044  // pointer to double 1045  // gnomonic projection 1046  double c2[2]; 1047  IntxUtils::gnomonic_projection( dp, Rsrc, pl, c2[0], c2[1] ); // 2 coordinates 1048  for( int j = 0; j < 2; j++ ) 1049  { 1050  if( qmin[j] > c2[j] ) qmin[j] = c2[j]; 1051  if( qmax[j] < c2[j] ) qmax[j] = c2[j]; 1052  } 1053  } 1054  // now decide if processor p should be interested in this cell, by looking at plane pl 1055  // 2d box this is one of the few size n loops; 1056  for( int p = 0; p < numprocs; p++ ) // each cell q can be sent to more than one processor 1057  { 1058  double procMin1 = allBoxes[24 * p + 4 * ( pl - 1 )]; // these were determined before 1059  // 1060  if( procMin1 >= DBL_MAX ) // the processor has no targets on this plane 1061  continue; 1062  double procMin2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 1]; 1063  double procMax1 = allBoxes[24 * p + 4 * ( pl - 1 ) + 2]; 1064  double procMax2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 3]; 1065  // test overlap of 2d boxes 1066  if( procMin1 > qmax[0] + box_error || procMin2 > qmax[1] + box_error ) continue; // 1067  if( qmin[0] > procMax1 + box_error || qmin[1] > procMax2 + box_error ) continue; 1068  // good to be inserted 1069  Rto[p].insert( q ); 1070  } 1071  } 1072  } 1073  else // regular 3d box; one box per processor 1074  { 1075  for( int p = 0; p < numprocs; p++ ) 1076  { 1077  BoundBox box( &allBoxes[6 * p] ); 1078  bool insert = false; 1079  for( int i = 0; i < num_nodes; i++ ) 1080  { 1081  if( box.contains_point( &elco[3 * i], box_error ) ) 1082  { 1083  insert = true; 1084  break; 1085  } 1086  } 1087  if( insert ) Rto[p].insert( q ); 1088  } 1089  } 1090  } 1091  1092  // here, we will use crystal router to send each cell to designated tasks (mesh migration) 1093  1094  // a better implementation would be to use pcomm send / recv entities; a good test case 1095  // pcomm send / receives uses point to point communication, not global gather / scatter 1096  1097  // now, build TLv and TLq (tuple list for vertices and cells, separately sent) 1098  size_t numq = 0; 1099  size_t numv = 0; 1100  1101  // merge the list of vertices to be sent 1102  for( int p = 0; p < numprocs; p++ ) 1103  { 1104  if( p == (int)my_rank ) continue; // do not "send" it to current task, because it is already here 1105  Range& range_to_P = Rto[p]; 1106  // add the vertices to it 1107  if( range_to_P.empty() ) continue; // nothing to send to proc p 1108 #ifdef VERBOSE 1109  std::cout << " proc : " << my_rank << " to proc " << p << " send " << range_to_P.size() << " cells " 1110  << " psize: " << range_to_P.psize() << "\n"; 1111 #endif 1112  Range vertsToP; 1113  rval = mb->get_connectivity( range_to_P, vertsToP );MB_CHK_SET_ERR( rval, "can't get connectivity" ); 1114  numq = numq + range_to_P.size(); 1115  numv = numv + vertsToP.size(); 1116  range_to_P.merge( vertsToP ); 1117  } 1118  1119  TupleList TLv; // send vertices with a different tuple list 1120  TupleList TLq; 1121  TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, 3 real coordinates 1122  TLv.enableWriteAccess(); 1123  1124  // add also GLOBAL_DOFS info, if found on the mesh cell; it should be found only on HOMME cells! 1125  int sizeTuple = 1126  2 + max_edges_1 + migrated_mesh + size_gdofs_tag; // max edges could be up to MAXEDGES :) for polygons 1127  TLq.initialize( sizeTuple, 0, 0, 0, 1128  numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v), plus 1129  // original sender if set (migrated mesh case) 1130  // we will not send the entity handle, global ID should be more than enough 1131  // we will not need more than 2B vertices TODO 2B vertices or cells 1132  // if we need more than 2B, we will need to use a different marker anyway (GLOBAL ID is not 1133  // enough then) 1134  1135  TLq.enableWriteAccess(); 1136 #ifdef VERBOSE 1137  std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n"; 1138 #endif 1139  1140  for( int to_proc = 0; to_proc < numprocs; to_proc++ ) 1141  { 1142  if( to_proc == (int)my_rank ) continue; 1143  Range& range_to_P = Rto[to_proc]; 1144  Range V = range_to_P.subset_by_type( MBVERTEX ); 1145  1146  for( Range::iterator it = V.begin(); it != V.end(); ++it ) 1147  { 1148  EntityHandle v = *it; 1149  int index = mesh_verts.index( v ); // 1150  assert( -1 != index ); 1151  int n = TLv.get_n(); // current size of tuple list 1152  TLv.vi_wr[2 * n] = to_proc; // send to processor 1153  TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the second_mesh_verts range 1154  TLv.vr_wr[3 * n] = coords_mesh[3 * index]; // departure position, of the node local_verts[i] 1155  TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1]; 1156  TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2]; 1157  TLv.inc_n(); // increment tuple list size 1158  } 1159  // also, prep the 2d cells for sending ... 1160  Range Q = range_to_P.subset_by_dimension( 2 ); 1161  for( Range::iterator it = Q.begin(); it != Q.end(); ++it ) 1162  { 1163  EntityHandle q = *it; // this is a second mesh cell (or src, lagrange set) 1164  int global_id; 1165  rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_SET_ERR( rval, "can't get gid for polygon" ); 1166  int n = TLq.get_n(); // current size 1167  TLq.vi_wr[sizeTuple * n] = to_proc; // 1168  TLq.vi_wr[sizeTuple * n + 1] = 1169  global_id; // global id of element, used to identify it for debug purposes only 1170  const EntityHandle* conn4; 1171  int num_nodes; // could be up to MAXEDGES; max_edges?; 1172  rval = mb->get_connectivity( q, conn4, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity for cell" ); 1173  if( num_nodes > max_edges_1 ) 1174  { 1175  mb->list_entities( &q, 1 ); 1176  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  // GLOBAL_DOFS info, if available 1200  if( size_gdofs_tag ) 1201  { 1202  rval = mb->tag_get_data( gdsTag, &q, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't get gdofs data in HOMME" ); 1203  for( int i = 0; i < size_gdofs_tag; i++ ) 1204  { 1205  TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] = 1206  valsDOFs[i]; // should be different than 0 or -1 1207  } 1208  } 1209  1210  TLq.inc_n(); // increment tuple list size 1211  } 1212  } // end for loop over total number of processors 1213  1214  // now we are done populating the tuples; route them to the appropriate processors 1215  // this does the communication magic 1216  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 ); 1217  ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 ); 1218  1219  // the first mesh elements are in localEnts; we do not need them at all 1220  1221  // maps from global ids to new vertex and cell handles, that are added 1222  1223  std::map< int, EntityHandle > globalID_to_vertex_handle; 1224  // we already have some vertices from second mesh set; they are already in the processor, even 1225  // before receiving other verts from neighbors this is an inverse map from gid to vertex handle, 1226  // which is local here, we do not want to duplicate vertices their identifier is the global ID!! 1227  // it must be unique per mesh ! (I mean, first mesh, source); gid for second mesh is not needed here 1228  int k = 0; 1229  for( Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ ) 1230  { 1231  globalID_to_vertex_handle[gids[k]] = *vit; 1232  } 1233  /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one? 1234  globalID_to_eh.clear(); // we need it now in case of extra work, to not duplicate cells 1235  1236  // now, look at every TLv, and see if we have to create a vertex there or not 1237  int n = TLv.get_n(); // the size of the points received 1238  for( int i = 0; i < n; i++ ) 1239  { 1240  int globalId = TLv.vi_rd[2 * i + 1]; 1241  if( globalID_to_vertex_handle.find( globalId ) == 1242  globalID_to_vertex_handle.end() ) // we do not have locally this vertex (yet) 1243  // so we have to create it, and add to the inverse map 1244  { 1245  EntityHandle new_vert; 1246  double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] }; 1247  rval = mb->create_vertex( dp_pos, new_vert );MB_CHK_SET_ERR( rval, "can't create new vertex " ); 1248  globalID_to_vertex_handle[globalId] = new_vert; // now add it to the map 1249  // set the GLOBAL ID tag on the new vertex 1250  rval = mb->tag_set_data( gid, &new_vert, 1, &globalId );MB_CHK_SET_ERR( rval, "can't set global ID tag on new vertex " ); 1251  } 1252  } 1253  1254  // now, all necessary vertices should be created 1255  // look in the local list of 2d cells for this proc, and add all those cells to covering set 1256  // also 1257  1258  Range& local = Rto[my_rank]; 1259  Range local_q = local.subset_by_dimension( 2 ); 1260  1261  for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it ) 1262  { 1263  EntityHandle q = *it; // these are from source cells, local 1264  int gid_el; 1265  rval = mb->tag_get_data( gid, &q, 1, &gid_el );MB_CHK_SET_ERR( rval, "can't get global id of cell " ); 1266  assert( gid_el >= 0 ); 1267  globalID_to_eh[gid_el] = q; // do we need this? yes, now we do; parent tags are now using it heavily 1268  rval = mb->tag_set_data( sendProcTag, &q, 1, &my_rank );MB_CHK_SET_ERR( rval, "can't set sender for cell" ); 1269  } 1270  1271  // now look at all elements received through; we do not want to duplicate them 1272  n = TLq.get_n(); // number of elements received by this processor 1273  // a cell should be received from one proc only; so why are we so worried about duplicated 1274  // elements? a vertex can be received from multiple sources, that is fine 1275  1276  for( int i = 0; i < n; i++ ) 1277  { 1278  int globalIdEl = TLq.vi_rd[sizeTuple * i + 1]; 1279  // int from_proc=TLq.vi_rd[sizeTuple * i ]; // we do not need from_proc anymore 1280  1281  // do we already have a cell with this global ID, represented? 1282  // yes, it could happen for extraWork ! 1283  if( globalID_to_eh.find( globalIdEl ) != globalID_to_eh.end() ) continue; 1284  // construct the conn triangle , quad or polygon 1285  EntityHandle new_conn[MAXEDGES]; // we should use std::vector with max_edges_1 1286  int nnodes = -1; 1287  for( int j = 0; j < max_edges_1; j++ ) 1288  { 1289  int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID 1290  if( vgid == 0 ) 1291  new_conn[j] = 0; // this can actually happen for polygon mesh (when we have less 1292  // number of vertices than max_edges) 1293  else 1294  { 1295  assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() ); 1296  new_conn[j] = globalID_to_vertex_handle[vgid]; 1297  nnodes = j + 1; // nodes are at the beginning, and are variable number 1298  } 1299  } 1300  EntityHandle new_element; 1301  // 1302  EntityType entType = MBQUAD; 1303  if( nnodes > 4 ) entType = MBPOLYGON; 1304  if( nnodes < 4 ) entType = MBTRI; 1305  rval = mb->create_element( entType, new_conn, nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element for second mesh " ); 1306  1307  globalID_to_eh[globalIdEl] = new_element; 1308  local_q.insert( new_element ); 1309  rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set gid for cell " ); 1310  int currentIndexIntTuple = 2 + max_edges_1; 1311  if( migrated_mesh ) 1312  { 1313  orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple]; 1314  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" ); 1315  currentIndexIntTuple++; // add one more 1316  } 1317  // store also the processor this coverage element came from 1318  int from_proc = TLq.vi_rd[sizeTuple * i]; 1319  rval = mb->tag_set_data( sendProcTag, &new_element, 1, &from_proc );MB_CHK_SET_ERR( rval, "can't set sender for cell" ); 1320  1321  // check if we need to retrieve and set GLOBAL_DOFS data 1322  if( size_gdofs_tag ) 1323  { 1324  for( int j = 0; j < size_gdofs_tag; j++ ) 1325  { 1326  valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j]; 1327  } 1328  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" ); 1329  } 1330  } 1331  1332  // now, add to the covering_set the elements created in the local_q range 1333  rval = mb->add_entities( covering_set, local_q );MB_CHK_SET_ERR( rval, "can't add entities to new mesh set " ); 1334 #ifdef VERBOSE 1335  std::cout << " proc " << my_rank << " add " << local_q.size() << " cells to covering set \n"; 1336 #endif 1337  return MB_SUCCESS; 1338 } 1339  1340 #endif // MOAB_HAVE_MPI 1341 //#undef VERBOSE 1342 #undef CHECK_CONVEXITY 1343  1344 } /* namespace moab */