MOAB: Mesh Oriented datABase  (version 5.4.1)
Intx2MeshOnSphere.cpp
Go to the documentation of this file.
00001 /*
00002  * Intx2MeshOnSphere.cpp
00003  *
00004  *  Created on: Oct 3, 2012
00005  */
00006 
00007 #ifdef _MSC_VER            /* windows */
00008 #define _USE_MATH_DEFINES  // For M_PI
00009 #endif
00010 
00011 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00012 #include "moab/IntxMesh/IntxUtils.hpp"
00013 #include "moab/GeomUtil.hpp"
00014 #ifdef MOAB_HAVE_MPI
00015 #include "moab/ParallelComm.hpp"
00016 #endif
00017 #include "MBTagConventions.hpp"
00018 
00019 #include <cassert>
00020 
00021 // #define ENABLE_DEBUG
00022 //#define CHECK_CONVEXITY
00023 namespace moab
00024 {
00025 
00026 Intx2MeshOnSphere::Intx2MeshOnSphere( Interface* mbimpl, IntxAreaUtils::AreaMethod amethod )
00027     : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
00028 {
00029 }
00030 
00031 Intx2MeshOnSphere::~Intx2MeshOnSphere() {}
00032 
00033 /*
00034  * return also the area for robustness verification
00035  */
00036 double Intx2MeshOnSphere::setup_tgt_cell( EntityHandle tgt, int& nsTgt )
00037 {
00038 
00039     // get coordinates of the target quad, to decide the gnomonic plane
00040     double cellArea = 0;
00041 
00042     int num_nodes;
00043     ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
00044 
00045     nsTgt = num_nodes;
00046     // account for possible padded polygons
00047     while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
00048         nsTgt--;
00049 
00050     // CartVect coords[4];
00051     rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
00052 
00053     CartVect middle = tgtCoords[0];
00054     for( int i = 1; i < nsTgt; i++ )
00055         middle += tgtCoords[i];
00056     middle = 1. / nsTgt * middle;
00057 
00058     IntxUtils::decide_gnomonic_plane( middle, plane );  // output the plane
00059     for( int j = 0; j < nsTgt; j++ )
00060     {
00061         // populate coords in the plane for intersection
00062         // they should be oriented correctly, positively
00063         rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
00064     }
00065 
00066     for( int j = 1; j < nsTgt - 1; j++ )
00067         cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
00068 
00069     // take target coords in order and compute area in plane
00070     return cellArea;
00071 }
00072 
00073 /* the elements are convex for sure, then do a gnomonic projection of both,
00074  *  compute intersection in the plane, then go back to the sphere for the points
00075  *  */
00076 ErrorCode Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
00077                                                                   EntityHandle src,
00078                                                                   double* P,
00079                                                                   int& nP,
00080                                                                   double& area,
00081                                                                   int markb[MAXEDGES],
00082                                                                   int markr[MAXEDGES],
00083                                                                   int& nsBlue,
00084                                                                   int& nsTgt,
00085                                                                   bool check_boxes_first )
00086 {
00087     // the area will be used from now on, to see how well we fill the target cell with polygons
00088     // the points will be at most 40; they will describe a convex patch, after the points will be
00089     // ordered and collapsed (eliminate doubles)
00090 
00091     // CartVect srccoords[4];
00092     int num_nodes  = 0;
00093     ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
00094     nsBlue = num_nodes;
00095     // account for possible padded polygons
00096     while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
00097         nsBlue--;
00098     rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
00099 
00100     area = 0.;
00101     nP   = 0;  // number of intersection points we are marking the boundary of src!
00102     if( check_boxes_first )
00103     {
00104         // look at the boxes formed with vertices; if they are far away, return false early
00105         // make sure the target is setup already
00106         setup_tgt_cell( tgt, nsTgt );  // we do not need area here
00107         // use here gnomonic plane (plane) to see where source is
00108         bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
00109         int planeb;
00110         CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
00111         IntxUtils::decide_gnomonic_plane( mid3, planeb );
00112         if( !overlap3d && ( plane != planeb ) )  // plane was set at setup_tgt_cell
00113             return MB_SUCCESS;                   // no error, but no intersection, decide early to get out
00114         // if same plane, still check for gnomonic plane in 2d
00115         // if no overlap in 2d, get out
00116         if( !overlap3d && plane == planeb )  // CHECK 2D too
00117         {
00118             for( int j = 0; j < nsBlue; j++ )
00119             {
00120                 rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j],
00121                                                        srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
00122             }
00123             bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
00124             if( !overlap2d ) return MB_SUCCESS;  // we are sure they are not overlapping in 2d , either
00125         }
00126     }
00127 #ifdef ENABLE_DEBUG
00128     if( dbg_1 )
00129     {
00130         std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
00131         for( int j = 0; j < nsTgt; j++ )
00132         {
00133             std::cout << tgtCoords[j] << "\n";
00134         }
00135         std::cout << "src " << mb->id_from_handle( src ) << "\n";
00136         for( int j = 0; j < nsBlue; j++ )
00137         {
00138             std::cout << srcCoords[j] << "\n";
00139         }
00140         mb->list_entities( &tgt, 1 );
00141         mb->list_entities( &src, 1 );
00142     }
00143 #endif
00144 
00145     for( int j = 0; j < nsBlue; j++ )
00146     {
00147         rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
00148     }
00149 
00150 #ifdef ENABLE_DEBUG
00151     if( dbg_1 )
00152     {
00153         std::cout << "gnomonic plane: " << plane << "\n";
00154         std::cout << " target                                src\n";
00155         for( int j = 0; j < nsTgt; j++ )
00156         {
00157             std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
00158         }
00159         for( int j = 0; j < nsBlue; j++ )
00160         {
00161             std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
00162         }
00163     }
00164 #endif
00165 
00166     rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
00167 
00168     int side[MAXEDGES] = { 0 };  // this refers to what side? source or tgt?
00169     int extraPoints =
00170         IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
00171     if( extraPoints >= 1 )
00172     {
00173         for( int k = 0; k < nsBlue; k++ )
00174         {
00175             if( side[k] )
00176             {
00177                 // this means that vertex k of source is inside convex tgt; mark edges k-1 and k in
00178                 // src,
00179                 //   as being "intersected" by tgt; (even though they might not be intersected by
00180                 //   other edges, the fact that their apex is inside, is good enough)
00181                 markb[k] = 1;
00182                 markb[( k + nsBlue - 1 ) % nsBlue] =
00183                     1;  // it is the previous edge, actually, but instead of doing -1, it is
00184                 // better to do modulo +3 (modulo 4)
00185                 // null side b for next call
00186                 side[k] = 0;
00187             }
00188         }
00189     }
00190     nP += extraPoints;
00191 
00192     extraPoints =
00193         IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
00194     if( extraPoints >= 1 )
00195     {
00196         for( int k = 0; k < nsTgt; k++ )
00197         {
00198             if( side[k] )
00199             {
00200                 // this is to mark that target edges k-1 and k are intersecting src
00201                 markr[k] = 1;
00202                 markr[( k + nsTgt - 1 ) % nsTgt] =
00203                     1;  // it is the previous edge, actually, but instead of doing -1, it is
00204                 // better to do modulo +3 (modulo 4)
00205                 // null side b for next call
00206             }
00207         }
00208     }
00209     nP += extraPoints;
00210 
00211     // now sort and orient the points in P, such that they are forming a convex polygon
00212     // this will be the foundation of our new mesh
00213     // this works if the polygons are convex
00214     IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 );  // nP should be at most 8 in the end ?
00215     // if there are more than 3 points, some area will be positive
00216 
00217     if( nP >= 3 )
00218     {
00219         for( int k = 1; k < nP - 1; k++ )
00220             area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
00221 #ifdef CHECK_CONVEXITY
00222         // each edge should be large enough that we can compute angles between edges
00223         for( int k = 0; k < nP; k++ )
00224         {
00225             int k1              = ( k + 1 ) % nP;
00226             int k2              = ( k1 + 1 ) % nP;
00227             double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
00228             if( orientedArea < 0 )
00229             {
00230                 std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
00231                           << " " << src << " \n";
00232             }
00233         }
00234 #endif
00235     }
00236 
00237     return MB_SUCCESS;  // no error
00238 }
00239 
00240 // this method will also construct the triangles/quads/polygons in the new mesh
00241 // if we accept planar polygons, we just save them
00242 // also, we could just create new vertices every time, and merge only in the end;
00243 // could be too expensive, and the tolerance for merging could be an
00244 // interesting topic
00245 ErrorCode Intx2MeshOnSphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsBlue, double* iP, int nP )
00246 {
00247 #ifdef ENABLE_DEBUG
00248     // first of all, check against target and source vertices
00249     //
00250     if( dbg_1 )
00251     {
00252         std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
00253                   << "\n";
00254         for( int n = 0; n < nP; n++ )
00255             std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
00256     }
00257 #endif
00258 
00259     // get the edges for the target triangle; the extra points will be on those edges, saved as
00260     // lists (unordered)
00261 
00262     // first get the list of edges adjacent to the target cell
00263     // use the neighTgtEdgeTag
00264     EntityHandle adjTgtEdges[MAXEDGES];
00265     ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge target tag" );
00266     // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
00267     // potatoes some of them will be handles to the initial vertices from source or target meshes
00268 
00269     std::vector< EntityHandle > foundIds;
00270     foundIds.resize( nP );
00271 #ifdef CHECK_CONVEXITY
00272     int npBefore1 = nP;
00273     int oldNodes  = 0;
00274     int otherIntx = 0;
00275 #endif
00276     for( int i = 0; i < nP; i++ )
00277     {
00278         double* pp = &iP[2 * i];  // iP+2*i
00279         // project the point back on the sphere
00280         CartVect pos;
00281         IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
00282         int found = 0;
00283         // first, are they on vertices from target or src?
00284         // priority is the target mesh (mb2?)
00285         int j                = 0;
00286         EntityHandle outNode = (EntityHandle)0;
00287         for( j = 0; j < nsTgt && !found; j++ )
00288         {
00289             // int node = tgtTri.v[j];
00290             double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
00291             if( d2 < epsilon_1 )
00292             {
00293 
00294                 foundIds[i] = tgtConn[j];  // no new node
00295                 found       = 1;
00296 #ifdef CHECK_CONVEXITY
00297                 oldNodes++;
00298 #endif
00299 #ifdef ENABLE_DEBUG
00300                 if( dbg_1 )
00301                     std::cout << "  target node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
00302                               << " 2d coords:" << tgtCoords2D[2 * j] << "  " << tgtCoords2D[2 * j + 1] << " d2: " << d2
00303                               << " \n";
00304 #endif
00305             }
00306         }
00307 
00308         for( j = 0; j < nsBlue && !found; j++ )
00309         {
00310             // int node = srcTri.v[j];
00311             double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
00312             if( d2 < epsilon_1 )
00313             {
00314                 // suspect is srcConn[j] corresponding in mbOut
00315 
00316                 foundIds[i] = srcConn[j];  // no new node
00317                 found       = 1;
00318 #ifdef CHECK_CONVEXITY
00319                 oldNodes++;
00320 #endif
00321 #ifdef ENABLE_DEBUG
00322                 if( dbg_1 )
00323                     std::cout << "  source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
00324                               << " \n";
00325 #endif
00326             }
00327         }
00328 
00329         if( !found )
00330         {
00331             // find the edge it belongs, first, on the red element
00332             // look at the minimum area, not at the first below some tolerance
00333             double minArea = 1.e+38;
00334             int index_min  = -1;
00335             for( j = 0; j < nsTgt; j++ )
00336             {
00337                 int j1      = ( j + 1 ) % nsTgt;
00338                 double area = fabs( IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp ) );
00339                 // how to check if pp is between redCoords2D[j] and redCoords2D[j1] ?
00340                 // they should form a straight line; the sign should be -1
00341                 double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
00342                                 IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
00343                                 IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
00344                 if( area < minArea && checkx < 2 * epsilon_1 )  // round off error or not?
00345                 {
00346                     index_min = j;
00347                     minArea   = area;
00348                 }
00349             }
00350             // verify that index_min is valid
00351             assert( index_min >= 0 );
00352 
00353             if( minArea < epsilon_1 / 2 )  // we found the smallest area, so we think we found the
00354                                            // target edge it belongs
00355             {
00356                 // found the edge; now find if there is a point in the list here
00357                 // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
00358                 int indx = TgtEdges.index( adjTgtEdges[index_min] );
00359                 if( indx < 0 )  // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
00360                 {
00361                     std::cerr << " error in adjacent target edge: " << mb->id_from_handle( adjTgtEdges[index_min] )
00362                               << "\n";
00363                     return MB_FAILURE;
00364                 }
00365                 std::vector< EntityHandle >* expts = extraNodesVec[indx];
00366                 // if the points pp is between extra points, then just give that id
00367                 // if not, create a new point, (check the id)
00368                 // get the coordinates of the extra points so far
00369                 int nbExtraNodesSoFar = expts->size();
00370                 if( nbExtraNodesSoFar > 0 )
00371                 {
00372                     std::vector< CartVect > coords1;
00373                     coords1.resize( nbExtraNodesSoFar );
00374                     mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
00375                     // std::list<int>::iterator it;
00376                     for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
00377                     {
00378                         // int pnt = *it;
00379                         double d2 = ( pos - coords1[k] ).length();
00380                         if( d2 < 2 * epsilon_1 )  // is this below machine precision?
00381                         {
00382                             found       = 1;
00383                             foundIds[i] = ( *expts )[k];
00384 #ifdef CHECK_CONVEXITY
00385                             otherIntx++;
00386 #endif
00387                         }
00388                     }
00389                 }
00390                 if( !found )
00391                 {
00392                     // create a new point in 2d (at the intersection)
00393                     // foundIds[i] = m_num2dPoints;
00394                     // expts.push_back(m_num2dPoints);
00395                     // need to create a new node in mbOut
00396                     // this will be on the edge, and it will be added to the local list
00397                     rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
00398                     ( *expts ).push_back( outNode );
00399                     // CID 181168; avoid leak storage error
00400                     rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
00401                     foundIds[i] = outNode;
00402                     found       = 1;
00403                 }
00404             }
00405         }
00406         if( !found )
00407         {
00408             std::cout << " target quad: ";
00409             for( int j1 = 0; j1 < nsTgt; j1++ )
00410             {
00411                 std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
00412             }
00413             std::cout << " a point pp is not on a target quad " << *pp << " " << pp[1] << " target quad "
00414                       << mb->id_from_handle( tgt ) << " \n";
00415             return MB_FAILURE;
00416         }
00417     }
00418 #ifdef ENABLE_DEBUG
00419     if( dbg_1 )
00420     {
00421         std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
00422         for( int i1 = 0; i1 < nP; i1++ )
00423             std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
00424     }
00425 #endif
00426     // first, find out if we have nodes collapsed; shrink them
00427     // we may have to reduce nP
00428     // it is possible that some nodes are collapsed after intersection only
00429     // nodes will always be in order (convex intersection)
00430 #ifdef CHECK_CONVEXITY
00431     int npBefore2 = nP;
00432 #endif
00433     correct_polygon( &foundIds[0], nP );
00434     // now we can build the triangles, from P array, with foundIds
00435     // we will put them in the out set
00436     if( nP >= 3 )
00437     {
00438         EntityHandle polyNew;
00439         rval = mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew );MB_CHK_ERR( rval );
00440         rval = mb->add_entities( outSet, &polyNew, 1 );MB_CHK_ERR( rval );
00441 
00442         // tag it with the global ids from target and source elements
00443         int globalID;
00444         rval = mb->tag_get_data( gid, &src, 1, &globalID );MB_CHK_ERR( rval );
00445         rval = mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
00446         // if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
00447         // " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";
00448         rval = mb->tag_get_data( gid, &tgt, 1, &globalID );MB_CHK_ERR( rval );
00449         rval = mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
00450         // if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
00451         // " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";
00452 
00453         counting++;
00454         rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
00455         if( orgSendProcTag )
00456         {
00457             int org_proc = -1;
00458             rval         = mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc );MB_CHK_ERR( rval );
00459             rval = mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc );MB_CHK_ERR( rval );  // yet another tag
00460         }
00461 #ifdef CHECK_CONVEXITY
00462         // each edge should be large enough that we can compute angles between edges
00463         std::vector< double > coords;
00464         coords.resize( 3 * nP );
00465         rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
00466         std::vector< CartVect > posi( nP );
00467         rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
00468 
00469         for( int k = 0; k < nP; k++ )
00470         {
00471             int k1 = ( k + 1 ) % nP;
00472             int k2 = ( k1 + 1 ) % nP;
00473             double orientedArea =
00474                 area_spherical_triangle_lHuiller( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
00475             if( orientedArea < 0 )
00476             {
00477                 std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
00478                 for( int i = 0; i < nP; i++ )
00479                 {
00480                     int nexti         = ( i + 1 ) % nP;
00481                     double lengthEdge = ( posi[i] - posi[nexti] ).length();
00482                     std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n";
00483                 }
00484                 std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n";
00485 
00486                 std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
00487                           << " target, src:" << tgt << " " << src << " \n";
00488             }
00489         }
00490 #endif
00491 
00492 #ifdef ENABLE_DEBUG
00493         if( dbg_1 )
00494         {
00495             std::cout << "Counting: " << counting << "\n";
00496             std::cout << " polygon " << mb->id_from_handle( polyNew ) << "  nodes: " << nP << " :";
00497             for( int i1 = 0; i1 < nP; i1++ )
00498                 std::cout << " " << mb->id_from_handle( foundIds[i1] );
00499             std::cout << " plane: " << plane << "\n";
00500             std::vector< CartVect > posi( nP );
00501             mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
00502             for( int i1 = 0; i1 < nP; i1++ )
00503                 std::cout << foundIds[i1] << " " << posi[i1] << "\n";
00504 
00505             std::stringstream fff;
00506             fff << "file0" << counting << ".vtk";
00507             rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
00508         }
00509 #endif
00510     }
00511     // else {
00512     //   std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";
00513     // }
00514     // disable_debug();
00515     return MB_SUCCESS;
00516 }
00517 
00518 ErrorCode Intx2MeshOnSphere::update_tracer_data( EntityHandle out_set, Tag& tagElem, Tag& tagArea )
00519 {
00520     EntityHandle dum = 0;
00521 
00522     Tag corrTag;
00523     ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE,
00524                                          &dum );  // it should have been created
00525     MB_CHK_SET_ERR( rval, "can't get correlation tag" );
00526 
00527     // get all polygons out of out_set; then see where are they coming from
00528     Range polys;
00529     rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
00530 
00531     // rs2 is the target range, arrival; rs1 is src, departure;
00532     // there is a connection between rs1 and rs2, through the corrTag
00533     // corrTag is __correlation
00534     // basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly);
00535     // also,  mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly);
00536     // we start from rs2 existing, then we have to update something
00537 
00538     // tagElem will have multiple tracers
00539     int numTracers = 0;
00540     rval           = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
00541     if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
00542 
00543     std::vector< double > currentVals( rs2.size() * numTracers );
00544     rval = mb->tag_get_data( tagElem, rs2, &currentVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" );
00545 
00546     // create new tuple list for tracers to other processors, from remote_cells
00547 #ifdef MOAB_HAVE_MPI
00548     if( remote_cells )
00549     {
00550         int n = remote_cells->get_n();
00551         if( n > 0 )
00552         {
00553             remote_cells_with_tracers = new TupleList();
00554             remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
00555                                                    n );  // tracers are in these tuples
00556             remote_cells_with_tracers->enableWriteAccess();
00557             for( int i = 0; i < n; i++ )
00558             {
00559                 remote_cells_with_tracers->vi_wr[2 * i]     = remote_cells->vi_wr[2 * i];
00560                 remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
00561                 //    remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
00562                 remote_cells_with_tracers->vul_wr[i] =
00563                     remote_cells->vul_wr[i];  // this is the corresponding target cell (arrival)
00564                 for( int k = 0; k < numTracers; k++ )
00565                     remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0;  // initialize tracers to be transported
00566                 remote_cells_with_tracers->inc_n();
00567             }
00568         }
00569         delete remote_cells;
00570         remote_cells = NULL;
00571     }
00572 #endif
00573     // for each polygon, we have 2 indices: target and source parents
00574     // we need index source to update index tgt?
00575     std::vector< double > newValues( rs2.size() * numTracers,
00576                                      0. );  // initialize with 0 all of them
00577     // area of the polygon * conc on target (old) current quantity
00578     // finally, divide by the area of the tgt
00579     double check_intx_area = 0.;
00580     moab::IntxAreaUtils intxAreas( this->areaMethod );  // use_lHuiller = true
00581     for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
00582     {
00583         EntityHandle poly = *it;
00584         int srcIndex, tgtIndex;
00585         rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
00586 
00587         EntityHandle src = rs1[srcIndex - 1];  // big assumption, it should work for meshes where global id is the same
00588         // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes
00589         rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" );
00590         // EntityHandle target = rs2[tgtIndex];
00591         // big assumption here, target and source are "parallel" ;we should have an index from
00592         // source to target (so a deformed source corresponds to an arrival tgt)
00593         /// TODO: VSM: Its unclear whether we need the source or destination radius here.
00594         double radius = Rsrc;
00595         double areap  = intxAreas.area_spherical_element( mb, poly, radius );
00596         check_intx_area += areap;
00597         // so the departure cell at time t (srcIndex) covers a portion of a tgtCell
00598         // that quantity will be transported to the tgtCell at time t+dt
00599         // the source corresponds to a target arrival
00600         EntityHandle tgtArr;
00601         rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
00602         if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
00603         {
00604 #ifdef MOAB_HAVE_MPI
00605             if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
00606             // maybe the element is remote, from another processor
00607             int global_id_src;
00608             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" );
00609             // find the
00610             int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
00611             if( index_in_remote == -1 )
00612                 MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
00613             for( int k = 0; k < numTracers; k++ )
00614                 remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
00615                     currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
00616 #endif
00617         }
00618         else if( MB_SUCCESS == rval )
00619         {
00620             int arrTgtIndex = rs2.index( tgtArr );
00621             if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
00622             for( int k = 0; k < numTracers; k++ )
00623                 newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
00624         }
00625 
00626         else
00627             MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " );
00628     }
00629     // now, send back the remote_cells_with_tracers to the processors they came from, with the
00630     // updated values for the tracer mass in a cell
00631 #ifdef MOAB_HAVE_MPI
00632     if( remote_cells_with_tracers )
00633     {
00634         // so this means that some cells will be sent back with tracer info to the procs they were
00635         // sent from
00636         ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
00637         // now, look at the global id, find the proper "tgt" cell with that index and update its
00638         // mass
00639         // remote_cells->print("remote cells after routing");
00640         int n = remote_cells_with_tracers->get_n();
00641         for( int j = 0; j < n; j++ )
00642         {
00643             EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j];  // entity handle sent back
00644             int arrTgtIndex      = rs2.index( tgtCell );
00645             if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
00646             for( int k = 0; k < numTracers; k++ )
00647                 newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
00648         }
00649     }
00650 #endif /* MOAB_HAVE_MPI */
00651     // now divide by target area (current)
00652     int j                = 0;
00653     Range::iterator iter = rs2.begin();
00654     void* data           = NULL;  // used for stored area
00655     int count            = 0;
00656     std::vector< double > total_mass_local( numTracers, 0. );
00657     while( iter != rs2.end() )
00658     {
00659         rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
00660         double* ptrArea = (double*)data;
00661         for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
00662         {
00663             for( int k = 0; k < numTracers; k++ )
00664             {
00665                 total_mass_local[k] += newValues[j * numTracers + k];
00666                 newValues[j * numTracers + k] /= ( *ptrArea );
00667             }
00668         }
00669     }
00670     rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" );
00671 
00672 #ifdef MOAB_HAVE_MPI
00673     std::vector< double > total_mass( numTracers, 0. );
00674     double total_intx_area = 0;
00675     int mpi_err =
00676         MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00677     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00678     // now reduce total area
00679     mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00680     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00681     if( my_rank == 0 )
00682     {
00683         for( int k = 0; k < numTracers; k++ )
00684             std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n";
00685         std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
00686                   << total_intx_area << "\n";
00687     }
00688 
00689     if( remote_cells_with_tracers )
00690     {
00691         delete remote_cells_with_tracers;
00692         remote_cells_with_tracers = NULL;
00693     }
00694 #else
00695     for( int k = 0; k < numTracers; k++ )
00696         std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n";
00697     std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
00698               << check_intx_area << "\n";
00699 #endif
00700     return MB_SUCCESS;
00701 }
00702 
00703 #ifdef MOAB_HAVE_MPI
00704 ErrorCode Intx2MeshOnSphere::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
00705 {
00706     localEnts.clear();
00707     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );MB_CHK_SET_ERR( rval, "can't get local ents" );
00708 
00709     rval = mb->get_connectivity( localEnts, local_verts );MB_CHK_SET_ERR( rval, "can't get connectivity" );
00710     int num_local_verts = (int)local_verts.size();
00711 
00712     assert( parcomm != NULL );
00713 
00714     if( num_local_verts == 0 )
00715     {
00716         // it is probably point cloud, get the local vertices from set
00717         rval = mb->get_entities_by_dimension( euler_set, 0, local_verts );MB_CHK_SET_ERR( rval, "can't get local vertices from set" );
00718         num_local_verts = (int)local_verts.size();
00719         localEnts       = local_verts;
00720     }
00721     // will use 6 gnomonic planes to decide boxes for each gnomonic plane
00722     // each gnomonic box will be 2d, min, max
00723     double gnom_box[24];
00724     for( int i = 0; i < 6; i++ )
00725     {
00726         gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
00727         gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
00728     }
00729 
00730     // there are 6 gnomonic planes; some elements could be on the corners, and affect multiple
00731     // planes decide what gnomonic planes will be affected by each cell some elements could appear
00732     // in multiple gnomonic planes !
00733     std::vector< double > coords( 3 * num_local_verts );
00734     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 " );
00735     // decide each local vertex to what gnomonic plane it belongs
00736 
00737     std::vector< int > gnplane;
00738     gnplane.resize( num_local_verts );
00739     for( int i = 0; i < num_local_verts; i++ )
00740     {
00741         CartVect pos( &coords[3 * i] );
00742         int pl;
00743         IntxUtils::decide_gnomonic_plane( pos, pl );
00744         gnplane[i] = pl;
00745     }
00746 
00747     for( Range::iterator it = localEnts.begin(); it != localEnts.end(); it++ )
00748     {
00749         EntityHandle cell   = *it;
00750         EntityType typeCell = mb->type_from_handle( cell );  // could be vertex, for point cloud
00751         // get coordinates, and decide gnomonic planes for it
00752         int nnodes;
00753         const EntityHandle* conn = NULL;
00754         EntityHandle c[1];
00755         if( typeCell != MBVERTEX )
00756         {
00757             rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
00758         }
00759         else
00760         {
00761             nnodes = 1;
00762             c[0]   = cell;  // actual node
00763             conn   = &c[0];
00764         }
00765         // get coordinates of vertices involved with this
00766         std::vector< double > elco( 3 * nnodes );
00767         std::set< int > planes;
00768         for( int i = 0; i < nnodes; i++ )
00769         {
00770             int ix = local_verts.index( conn[i] );
00771             planes.insert( gnplane[ix] );
00772             for( int j = 0; j < 3; j++ )
00773             {
00774                 elco[3 * i + j] = coords[3 * ix + j];
00775             }
00776         }
00777         // now, augment the boxes for all planes involved
00778         for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
00779         {
00780             int pl = *st;
00781             for( int i = 0; i < nnodes; i++ )
00782             {
00783                 CartVect pos( &elco[3 * i] );
00784                 double c2[2];
00785                 IntxUtils::gnomonic_projection( pos, Rdest, pl, c2[0],
00786                                                 c2[1] );  // 2 coordinates
00787                 //
00788                 for( int k = 0; k < 2; k++ )
00789                 {
00790                     double val = c2[k];
00791                     if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val;  // min in k direction
00792                     if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
00793                         gnom_box[4 * ( pl - 1 ) + 2 + k] = val;  // max in k direction
00794                 }
00795             }
00796         }
00797     }
00798 
00799     int numprocs = parcomm->proc_config().proc_size();
00800     allBoxes.resize( 24 * numprocs );  // 6 gnomonic planes , 4 doubles for each for 2d box
00801 
00802     my_rank = parcomm->proc_config().proc_rank();
00803     for( int k = 0; k < 24; k++ )
00804         allBoxes[24 * my_rank + k] = gnom_box[k];
00805 
00806     // now communicate to get all boxes
00807     int mpi_err;
00808 #if( MPI_VERSION >= 2 )
00809     // use "in place" option
00810     mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 24, MPI_DOUBLE,
00811                              parcomm->proc_config().proc_comm() );
00812 #else
00813     {
00814         std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
00815         mpi_err  = MPI_Allgather( &allBoxes[24 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
00816                                   parcomm->proc_config().proc_comm() );
00817         allBoxes = allBoxes_tmp;
00818     }
00819 #endif
00820     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00821 
00822 #ifdef VERBOSE
00823     if( my_rank == 0 )
00824     {
00825         std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
00826                   << " on second mesh \n";
00827         for( int i = 0; i < numprocs; i++ )
00828         {
00829             std::cout << "task: " << i << " \n";
00830             for( int pl = 1; pl <= 6; pl++ )
00831             {
00832                 std::cout << "  plane " << pl << " min: \t" << allBoxes[24 * i + 4 * ( pl - 1 )] << " \t"
00833                           << allBoxes[24 * i + 4 * ( pl - 1 ) + 1] << "\n";
00834                 std::cout << " \t  max: \t" << allBoxes[24 * i + 4 * ( pl - 1 ) + 2] << " \t"
00835                           << allBoxes[24 * i + 4 * ( pl - 1 ) + 3] << "\n";
00836             }
00837         }
00838     }
00839 #endif
00840 
00841     return MB_SUCCESS;
00842 }
00843 //#define VERBOSE
00844 // this will use the bounding boxes for the (euler)/ fix  mesh that are already established
00845 // will distribute the mesh to other procs, so that on each task, the covering set covers the local
00846 // bounding box this means it will cover the second (local) mesh set; So the covering set will cover
00847 // completely the second local mesh set (in intersection)
00848 ErrorCode Intx2MeshOnSphere::construct_covering_set( EntityHandle& initial_distributed_set, EntityHandle& covering_set )
00849 {
00850     // primary element came from, in the joint communicator ; this will be forwarded by coverage
00851     // mesh needed for tag migrate later on
00852     int defaultInt = -1;  // no processor, so it was not migrated from somewhere else
00853     ErrorCode rval = mb->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag,
00854                                          MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" );
00855 
00856     assert( parcomm != NULL );
00857     Range meshCells;
00858     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" );
00859 
00860     if( 1 == parcomm->proc_config().proc_size() )
00861     {
00862         // move all initial cells to coverage set
00863         rval = mb->add_entities( covering_set, meshCells );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
00864         // if point cloud source, add vertices
00865         if( 0 == meshCells.size() || max_edges_1 == 0 )
00866         {
00867             // add vertices from the source set
00868             Range verts;
00869             rval = mb->get_entities_by_dimension( initial_distributed_set, 0, verts );MB_CHK_SET_ERR( rval, "can't get vertices from mesh set" );
00870             rval = mb->add_entities( covering_set, verts );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
00871         }
00872         return MB_SUCCESS;
00873     }
00874 
00875     // mark on the coverage mesh where this element came from
00876     Tag sendProcTag;  /// for coverage mesh, will store the sender
00877     rval = mb->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag, MB_TAG_DENSE | MB_TAG_CREAT,
00878                                &defaultInt );MB_CHK_SET_ERR( rval, "can't create sending processor tag" );
00879 
00880     // this information needs to be forwarded to coverage mesh, if this mesh was already migrated
00881     // from somewhere else
00882     // look at the value of orgSendProcTag for one mesh cell; if -1, no need to forward that; if
00883     // !=-1, we know that this mesh was migrated, we need to find out more about origin of cell
00884     int orig_sender      = -1;
00885     EntityHandle oneCell = 0;
00886     // decide if we need to transfer global DOFs info attached to each HOMME coarse cell; first we
00887     // need to decide if the mesh has that tag; will affect the size of the tuple list involved in
00888     // the crystal routing
00889     int size_gdofs_tag = 0;
00890     std::vector< int > valsDOFs;
00891     Tag gdsTag;
00892     rval = mb->tag_get_handle( "GLOBAL_DOFS", gdsTag );
00893 
00894     if( meshCells.size() > 0 )
00895     {
00896         oneCell = meshCells[0];  // it is possible we do not have any cells, even after migration
00897         rval    = mb->tag_get_data( orgSendProcTag, &oneCell, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sending processor value" );
00898         if( gdsTag )
00899         {
00900             DataType dtype;
00901             rval = mb->tag_get_data_type( gdsTag, dtype );
00902             if( MB_SUCCESS == rval && MB_TYPE_INTEGER == dtype )
00903             {
00904                 // find the values on first cell
00905                 int lenTag = 0;
00906                 rval       = mb->tag_get_length( gdsTag, lenTag );
00907                 if( MB_SUCCESS == rval && lenTag > 0 )
00908                 {
00909                     valsDOFs.resize( lenTag );
00910                     rval = mb->tag_get_data( gdsTag, &oneCell, 1, &valsDOFs[0] );
00911                     if( MB_SUCCESS == rval && valsDOFs[0] > 0 )
00912                     {
00913                         // first value positive means we really need to transport this data during
00914                         // coverage
00915                         size_gdofs_tag = lenTag;
00916                     }
00917                 }
00918             }
00919         }
00920     }
00921 
00922     // another collective call, to see if the mesh is migrated and if the GLOBAL_DOFS tag need to be
00923     // transferred over to the coverage mesh it is possible that there is no initial mesh source
00924     // mesh on the task, so we do not know that info from the tag but TupleList needs to be sized
00925     // uniformly for all tasks do a collective MPI_MAX to see if it is migrated and if we have
00926     // (collectively) a GLOBAL_DOFS task
00927 
00928     int local_int_array[2], global_int_array[2];
00929     local_int_array[0] = orig_sender;
00930     local_int_array[1] = size_gdofs_tag;
00931     // now reduce over all processors
00932     int mpi_err =
00933         MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
00934     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00935     orig_sender    = global_int_array[0];
00936     size_gdofs_tag = global_int_array[1];
00937 #ifdef VERBOSE
00938     std::cout << "proc: " << my_rank << " size_gdofs_tag:" << size_gdofs_tag << "\n";
00939 #endif
00940     valsDOFs.resize( size_gdofs_tag );
00941 
00942     // finally, we have correct global info, we can decide if mesh was migrated and if we have
00943     // global dofs tag that need to be sent with coverage info
00944     int migrated_mesh = 0;
00945     if( orig_sender != -1 ) migrated_mesh = 1;  //
00946     // if size_gdofs_tag>0, we are sure valsDOFs got resized to what we need
00947 
00948     // get all mesh verts1
00949     Range mesh_verts;
00950     rval = mb->get_connectivity( meshCells, mesh_verts );MB_CHK_SET_ERR( rval, "can't get  mesh vertices" );
00951     int num_mesh_verts = (int)mesh_verts.size();
00952 
00953     // now see the mesh points positions; to what boxes should we send them?
00954     std::vector< double > coords_mesh( 3 * num_mesh_verts );
00955     rval = mb->get_coords( mesh_verts, &coords_mesh[0] );MB_CHK_SET_ERR( rval, "can't get mesh points position" );
00956 
00957     // decide gnomonic plane for each vertex, as in the compute boxes
00958     std::vector< int > gnplane;
00959     gnplane.resize( num_mesh_verts );
00960     for( int i = 0; i < num_mesh_verts; i++ )
00961     {
00962         CartVect pos( &coords_mesh[3 * i] );
00963         int pl;
00964         IntxUtils::decide_gnomonic_plane( pos, pl );
00965         gnplane[i] = pl;
00966     }
00967 
00968     std::vector< int > gids( num_mesh_verts );
00969     rval = mb->tag_get_data( gid, mesh_verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get vertices gids" );
00970 
00971     // ranges to send to each processor; will hold vertices and elements (quads/ polygons)
00972     // will look if the box of the mesh cell covers bounding box(es) (within tolerances)
00973     std::map< int, Range > Rto;
00974     int numprocs = parcomm->proc_config().proc_size();
00975 
00976     for( Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
00977     {
00978         EntityHandle q = *eit;
00979         const EntityHandle* conn;
00980         int num_nodes;
00981         rval = mb->get_connectivity( q, conn, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity on cell" );
00982 
00983         // first decide what planes need to consider
00984         std::set< int > planes;  // if this list contains more than 3 planes, we have a very bad mesh!!!
00985         std::vector< double > elco( 3 * num_nodes );
00986         for( int i = 0; i < num_nodes; i++ )
00987         {
00988             EntityHandle v = conn[i];
00989             int index      = mesh_verts.index( v );
00990             planes.insert( gnplane[index] );
00991             for( int j = 0; j < 3; j++ )
00992             {
00993                 elco[3 * i + j] = coords_mesh[3 * index + j];  // extract from coords
00994             }
00995         }
00996         // now loop over all planes that need to be considered for this element
00997         for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
00998         {
00999             int pl         = *st;  // gnomonic plane considered
01000             double qmin[2] = { DBL_MAX, DBL_MAX };
01001             double qmax[2] = { -DBL_MAX, -DBL_MAX };
01002             for( int i = 0; i < num_nodes; i++ )
01003             {
01004                 CartVect dp( &elco[3 * i] );  // uses constructor for CartVect that takes a
01005                                               // pointer to double
01006                 // gnomonic projection
01007                 double c2[2];
01008                 IntxUtils::gnomonic_projection( dp, Rsrc, pl, c2[0], c2[1] );  // 2 coordinates
01009                 for( int j = 0; j < 2; j++ )
01010                 {
01011                     if( qmin[j] > c2[j] ) qmin[j] = c2[j];
01012                     if( qmax[j] < c2[j] ) qmax[j] = c2[j];
01013                 }
01014             }
01015             // now decide if processor p should be interested in this cell, by looking at plane pl
01016             // 2d box this is one of the few size n loops;
01017             for( int p = 0; p < numprocs; p++ )  // each cell q can be sent to more than one processor
01018             {
01019                 double procMin1 = allBoxes[24 * p + 4 * ( pl - 1 )];  // these were determined before
01020                 //
01021                 if( procMin1 >= DBL_MAX )  // the processor has no targets on this plane
01022                     continue;
01023                 double procMin2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
01024                 double procMax1 = allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
01025                 double procMax2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
01026                 // test overlap of 2d boxes
01027                 if( procMin1 > qmax[0] + box_error || procMin2 > qmax[1] + box_error ) continue;  //
01028                 if( qmin[0] > procMax1 + box_error || qmin[1] > procMax2 + box_error ) continue;
01029                 // good to be inserted
01030                 Rto[p].insert( q );
01031             }
01032         }
01033     }
01034 
01035     // here, we will use crystal router to send each cell to designated tasks (mesh migration)
01036 
01037     // a better implementation would be to use pcomm send / recv entities; a good test case
01038     // pcomm send / receives uses point to point communication, not global gather / scatter
01039 
01040     // now, build TLv and TLq  (tuple list for vertices and cells, separately sent)
01041     size_t numq = 0;
01042     size_t numv = 0;
01043 
01044     // merge the list of vertices to be sent
01045     for( int p = 0; p < numprocs; p++ )
01046     {
01047         if( p == (int)my_rank ) continue;  // do not "send" it to current task, because it is already here
01048         Range& range_to_P = Rto[p];
01049         // add the vertices to it
01050         if( range_to_P.empty() ) continue;  // nothing to send to proc p
01051 #ifdef VERBOSE
01052         std::cout << " proc : " << my_rank << " to proc " << p << " send " << range_to_P.size() << " cells "
01053                   << " psize: " << range_to_P.psize() << "\n";
01054 #endif
01055         Range vertsToP;
01056         rval = mb->get_connectivity( range_to_P, vertsToP );MB_CHK_SET_ERR( rval, "can't get connectivity" );
01057         numq = numq + range_to_P.size();
01058         numv = numv + vertsToP.size();
01059         range_to_P.merge( vertsToP );
01060     }
01061 
01062     TupleList TLv;  // send vertices with a different tuple list
01063     TupleList TLq;
01064     TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, 3 real coordinates
01065     TLv.enableWriteAccess();
01066 
01067     // add also GLOBAL_DOFS info, if found on the mesh cell; it should be found only on HOMME cells!
01068     int sizeTuple =
01069         2 + max_edges_1 + migrated_mesh + size_gdofs_tag;  // max edges could be up to MAXEDGES :) for polygons
01070     TLq.initialize( sizeTuple, 0, 0, 0,
01071                     numq );  // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v), plus
01072                              // original sender if set (migrated mesh case)
01073     // we will not send the entity handle, global ID should be more than enough
01074     // we will not need more than 2B vertices
01075     // if we need more than 2B, we will need to use a different marker anyway (GLOBAL ID is not
01076     // enough then)
01077 
01078     TLq.enableWriteAccess();
01079 #ifdef VERBOSE
01080     std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
01081 #endif
01082 
01083     for( int to_proc = 0; to_proc < numprocs; to_proc++ )
01084     {
01085         if( to_proc == (int)my_rank ) continue;
01086         Range& range_to_P = Rto[to_proc];
01087         Range V           = range_to_P.subset_by_type( MBVERTEX );
01088 
01089         for( Range::iterator it = V.begin(); it != V.end(); ++it )
01090         {
01091             EntityHandle v = *it;
01092             int index      = mesh_verts.index( v );  //
01093             assert( -1 != index );
01094             int n                = TLv.get_n();             // current size of tuple list
01095             TLv.vi_wr[2 * n]     = to_proc;                 // send to processor
01096             TLv.vi_wr[2 * n + 1] = gids[index];             // global id needs index in the second_mesh_verts range
01097             TLv.vr_wr[3 * n]     = coords_mesh[3 * index];  // departure position, of the node local_verts[i]
01098             TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
01099             TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
01100             TLv.inc_n();  // increment tuple list size
01101         }
01102         // also, prep the 2d cells for sending ...
01103         Range Q = range_to_P.subset_by_dimension( 2 );
01104         for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
01105         {
01106             EntityHandle q = *it;  // this is a second mesh cell (or src, lagrange set)
01107             int global_id;
01108             rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_SET_ERR( rval, "can't get gid for polygon" );
01109             int n                    = TLq.get_n();  // current size
01110             TLq.vi_wr[sizeTuple * n] = to_proc;      //
01111             TLq.vi_wr[sizeTuple * n + 1] =
01112                 global_id;  // global id of element, used to identify it for debug purposes only
01113             const EntityHandle* conn4;
01114             int num_nodes;  // could be up to MAXEDGES; max_edges?;
01115             rval = mb->get_connectivity( q, conn4, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity for cell" );
01116             if( num_nodes > max_edges_1 )
01117             {
01118                 mb->list_entities( &q, 1 );MB_CHK_SET_ERR( MB_FAILURE, "too many nodes in a cell (" << num_nodes << "," << max_edges_1 << ")" );
01119             }
01120             for( int i = 0; i < num_nodes; i++ )
01121             {
01122                 EntityHandle v = conn4[i];
01123                 int index      = mesh_verts.index( v );
01124                 assert( -1 != index );
01125                 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
01126             }
01127             for( int k = num_nodes; k < max_edges_1; k++ )
01128             {
01129                 TLq.vi_wr[sizeTuple * n + 2 + k] =
01130                     0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
01131             }
01132             int currentIndexIntTuple = 2 + max_edges_1;
01133             // is the mesh migrated before or not?
01134             if( migrated_mesh )
01135             {
01136                 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" );
01137                 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender;  // should be different than -1
01138                 currentIndexIntTuple++;
01139             }
01140             // GLOBAL_DOFS info, if available
01141             if( size_gdofs_tag )
01142             {
01143                 rval = mb->tag_get_data( gdsTag, &q, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't get gdofs data in HOMME" );
01144                 for( int i = 0; i < size_gdofs_tag; i++ )
01145                 {
01146                     TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
01147                         valsDOFs[i];  // should be different than 0 or -1
01148                 }
01149             }
01150 
01151             TLq.inc_n();  // increment tuple list size
01152         }
01153     }  // end for loop over total number of processors
01154 
01155     // now we are done populating the tuples; route them to the appropriate processors
01156     // this does the communication magic
01157     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01158     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
01159 
01160     // the first mesh elements are in localEnts; we do not need them at all
01161 
01162     // maps from global ids to new vertex and cell handles, that are added
01163 
01164     std::map< int, EntityHandle > globalID_to_vertex_handle;
01165     // we already have some vertices from second mesh set; they are already in the processor, even
01166     // before receiving other verts from neighbors this is an inverse map from gid to vertex handle,
01167     // which is local here, we do not want to duplicate vertices their identifier is the global ID!!
01168     // it must be unique per mesh ! (I mean, second mesh); gid gor first mesh is not needed here
01169     int k = 0;
01170     for( Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
01171     {
01172         globalID_to_vertex_handle[gids[k]] = *vit;
01173     }
01174     /*std::map<int, EntityHandle> globalID_to_eh;*/  // do we need this one?
01175     globalID_to_eh.clear();                          // we do not really need it, but we keep it for debugging mostly
01176 
01177     // now, look at every TLv, and see if we have to create a vertex there or not
01178     int n = TLv.get_n();  // the size of the points received
01179     for( int i = 0; i < n; i++ )
01180     {
01181         int globalId = TLv.vi_rd[2 * i + 1];
01182         if( globalID_to_vertex_handle.find( globalId ) ==
01183             globalID_to_vertex_handle.end() )  // we do not have locally this vertex (yet)
01184                                                // so we have to create it, and add to the inverse map
01185         {
01186             EntityHandle new_vert;
01187             double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
01188             rval             = mb->create_vertex( dp_pos, new_vert );MB_CHK_SET_ERR( rval, "can't create new vertex " );
01189             globalID_to_vertex_handle[globalId] = new_vert;  // now add it to the map
01190             // set the GLOBAL ID tag on the new vertex
01191             rval = mb->tag_set_data( gid, &new_vert, 1, &globalId );MB_CHK_SET_ERR( rval, "can't set global ID tag on new vertex " );
01192         }
01193     }
01194 
01195     // now, all necessary vertices should be created
01196     // look in the local list of 2d cells for this proc, and add all those cells to covering set
01197     // also
01198 
01199     Range& local  = Rto[my_rank];
01200     Range local_q = local.subset_by_dimension( 2 );
01201 
01202     for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
01203     {
01204         EntityHandle q = *it;  // these are from lagr cells, local
01205         int gid_el;
01206         rval = mb->tag_get_data( gid, &q, 1, &gid_el );MB_CHK_SET_ERR( rval, "can't get global id of cell " );
01207         assert( gid_el >= 0 );
01208         globalID_to_eh[gid_el] = q;  // do we need this? yes, now we do; parent tags are now using it heavily
01209         rval                   = mb->tag_set_data( sendProcTag, &q, 1, &my_rank );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
01210     }
01211 
01212     // now look at all elements received through; we do not want to duplicate them
01213     n = TLq.get_n();  // number of elements received by this processor
01214     // a cell should be received from one proc only; so why are we so worried about duplicated
01215     // elements? a vertex can be received from multiple sources, that is fine
01216 
01217     for( int i = 0; i < n; i++ )
01218     {
01219         int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
01220         // int from_proc=TLq.vi_rd[sizeTuple * i ]; // we do not need from_proc anymore
01221 
01222         // do we already have a quad with this global ID, represented? no way !
01223         // if (globalID_to_eh.find(globalIdEl) == globalID_to_eh.end())
01224         //{
01225         // construct the conn triangle , quad or polygon
01226         EntityHandle new_conn[MAXEDGES];  // we should use std::vector with max_edges_1
01227         int nnodes = -1;
01228         for( int j = 0; j < max_edges_1; j++ )
01229         {
01230             int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
01231             if( vgid == 0 )
01232                 new_conn[j] = 0;  // this can actually happen for polygon mesh (when we have less
01233                                   // number of vertices than max_edges)
01234             else
01235             {
01236                 assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
01237                 new_conn[j] = globalID_to_vertex_handle[vgid];
01238                 nnodes      = j + 1;  // nodes are at the beginning, and are variable number
01239             }
01240         }
01241         EntityHandle new_element;
01242         //
01243         EntityType entType = MBQUAD;
01244         if( nnodes > 4 ) entType = MBPOLYGON;
01245         if( nnodes < 4 ) entType = MBTRI;
01246         rval = mb->create_element( entType, new_conn, nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element for second mesh " );
01247 
01248         globalID_to_eh[globalIdEl] = new_element;
01249         local_q.insert( new_element );
01250         rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set gid for cell " );
01251         int currentIndexIntTuple = 2 + max_edges_1;
01252         if( migrated_mesh )
01253         {
01254             orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
01255             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" );
01256             currentIndexIntTuple++;  // add one more
01257         }
01258         // check if we need to retrieve and set GLOBAL_DOFS data
01259         if( size_gdofs_tag )
01260         {
01261             for( int j = 0; j < size_gdofs_tag; j++ )
01262             {
01263                 valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
01264             }
01265             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" );
01266         }
01267         // store also the processor this coverage element came from
01268         int from_proc = TLq.vi_rd[sizeTuple * i];
01269         rval          = mb->tag_set_data( sendProcTag, &new_element, 1, &from_proc );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
01270     }
01271 
01272     // now, create a new set, covering_set
01273     rval = mb->add_entities( covering_set, local_q );MB_CHK_SET_ERR( rval, "can't add entities to new mesh set " );
01274 #ifdef VERBOSE
01275     std::cout << " proc " << my_rank << " add " << local_q.size() << " cells to covering set \n";
01276 #endif
01277     return MB_SUCCESS;
01278 }
01279 
01280 #endif  // MOAB_HAVE_MPI
01281 //#undef VERBOSE
01282 #undef CHECK_CONVEXITY
01283 
01284 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines