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