MOAB: Mesh Oriented datABase  (version 5.4.1)
IntxRllCssphere.cpp
Go to the documentation of this file.
00001 /*
00002  * IntxRllCssphere.cpp
00003  *
00004  *  Created on: Aug 8, 2014
00005  *      Author: iulian
00006  */
00007 
00008 #include "moab/IntxMesh/IntxRllCssphere.hpp"
00009 #include "moab/GeomUtil.hpp"
00010 #include "moab/IntxMesh/IntxUtils.hpp"
00011 
00012 namespace moab
00013 {
00014 
00015 IntxRllCssphere::IntxRllCssphere( Interface* mbimpl ) : Intx2Mesh( mbimpl ), R( 0.0 ), plane( 0 ) {}
00016 
00017 IntxRllCssphere::~IntxRllCssphere() {}
00018 
00019 /*
00020  * return also the area for robustness verification
00021  */
00022 double IntxRllCssphere::setup_tgt_cell( EntityHandle tgt, int& nsTgt )
00023 {
00024     // get coordinates of the tgt quad, to decide the gnomonic plane
00025     double cellArea = 0;
00026 
00027     int num_nodes;
00028     ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );
00029 
00030     if( MB_SUCCESS != rval ) return 1;
00031     nsTgt = num_nodes;
00032     // these edges will never be polygons, only quads or triangles
00033 
00034     // CartVect coords[4];
00035     rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );
00036     if( MB_SUCCESS != rval ) return 1;
00037     CartVect middle = tgtCoords[0];
00038     for( int i = 1; i < nsTgt; i++ )
00039         middle += tgtCoords[i];
00040     middle = 1. / nsTgt * middle;
00041 
00042     IntxUtils::decide_gnomonic_plane( middle, plane );  // output the plane
00043     for( int j = 0; j < nsTgt; j++ )
00044     {
00045         // populate coords in the plane for intersection
00046         // they should be oriented correctly, positively
00047         int rc = IntxUtils::gnomonic_projection( tgtCoords[j], R, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );
00048         if( rc != 0 ) return 1;
00049     }
00050 
00051     for( int j = 1; j < nsTgt - 1; j++ )
00052         cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
00053 
00054     // take tgt coords in order and compute area in plane
00055     return cellArea;
00056 }
00057 
00058 /* the elements are convex for sure, then do a gnomonic projection of both,
00059  *  compute intersection in the plane, then go back to the sphere for the points
00060  *  */
00061 ErrorCode IntxRllCssphere::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
00062                                                                 EntityHandle src,
00063                                                                 double* P,
00064                                                                 int& nP,
00065                                                                 double& area,
00066                                                                 int markb[MAXEDGES],
00067                                                                 int markr[MAXEDGES],
00068                                                                 int& nsSrc,
00069                                                                 int& nsTgt,
00070                                                                 bool check_boxes_first )
00071 {
00072     // the area will be used from now on, to see how well we fill the tgt cell with polygons
00073     // the points will be at most 40; they will describe a convex patch, after the points will be
00074     // ordered and collapsed (eliminate doubles)
00075 
00076     // CartVect srccoords[4];
00077     int num_nodes  = 0;
00078     ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
00079 
00080     nsSrc = num_nodes;
00081     rval  = mb->get_coords( srcConn, nsSrc, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
00082 
00083     // determine the type of edge: const lat or not?
00084     // just look at the consecutive z coordinates for the edge
00085     for( int i = 0; i < nsSrc; i++ )
00086     {
00087         int nexti = ( i + 1 ) % nsSrc;
00088         if( fabs( srcCoords[i][2] - srcCoords[nexti][2] ) < 1.e-6 )
00089             srcEdgeType[i] = 1;
00090         else
00091             srcEdgeType[i] = 0;
00092     }
00093     area = 0.;
00094     nP   = 0;  // number of intersection points we are marking the boundary of src!
00095     if( check_boxes_first )
00096     {
00097         // look at the boxes formed with vertices; if they are far away, return false early
00098         // make sure the tgt is setup already
00099         setup_tgt_cell( tgt, nsTgt );  // we do not need area here
00100         if( !GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsSrc, box_error ) )
00101             return MB_SUCCESS;  // no error, but no intersection, decide early to get out
00102     }
00103 #ifdef ENABLE_DEBUG
00104     if( dbg_1 )
00105     {
00106         std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
00107         for( int j = 0; j < nsTgt; j++ )
00108         {
00109             std::cout << tgtCoords[j] << "\n";
00110         }
00111         std::cout << "src " << mb->id_from_handle( src ) << "\n";
00112         for( int j = 0; j < nsSrc; j++ )
00113         {
00114             std::cout << srcCoords[j] << "\n";
00115         }
00116         mb->list_entities( &tgt, 1 );
00117         mb->list_entities( &src, 1 );
00118     }
00119 #endif
00120     for( int j = 0; j < nsSrc; j++ )
00121     {
00122         rval = IntxUtils::gnomonic_projection( srcCoords[j], R, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
00123     }
00124 #ifdef ENABLE_DEBUG
00125     if( dbg_1 )
00126     {
00127         std::cout << "gnomonic plane: " << plane << "\n";
00128         std::cout << " tgt                                src\n";
00129         for( int j = 0; j < nsTgt; j++ )
00130         {
00131             std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
00132         }
00133         for( int j = 0; j < nsSrc; j++ )
00134         {
00135             std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
00136         }
00137     }
00138 #endif
00139     rval = IntxUtils::EdgeIntxRllCs( srcCoords2D, srcCoords, srcEdgeType, nsSrc, tgtCoords2D, tgtCoords, nsTgt, markb,
00140                                      markr, plane, R, P, nP );MB_CHK_ERR( rval );
00141 
00142     int side[MAXEDGES] = { 0 };  // this refers to what side? src or tgt?// more tolerant here with epsilon_area
00143     int extraPoints    = IntxUtils::borderPointsOfXinY2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, &( P[2 * nP] ), side,
00144                                                          2 * epsilon_area );
00145     if( extraPoints >= 1 )
00146     {
00147         for( int k = 0; k < nsSrc; k++ )
00148         {
00149             if( side[k] )
00150             {
00151                 // this means that vertex k of src is inside convex tgt; mark edges k-1 and k in
00152                 // src,
00153                 //   as being "intersected" by tgt; (even though they might not be intersected by
00154                 //   other edges, the fact that their apex is inside, is good enough)
00155                 markb[k] = 1;
00156                 markb[( k + nsSrc - 1 ) % nsSrc] =
00157                     1;  // it is the previous edge, actually, but instead of doing -1, it is
00158                 // better to do modulo +3 (modulo 4)
00159                 // null side b for next call
00160                 side[k] = 0;
00161             }
00162         }
00163     }
00164     nP += extraPoints;
00165 
00166     extraPoints =
00167         IntxUtils::borderPointsOfCSinRLL( tgtCoords, tgtCoords2D, nsTgt, srcCoords, nsSrc, srcEdgeType, &( P[2 * nP] ),
00168                                           side,
00169                                           100 * epsilon_area );  // we need to compare with 0 a volume from 3 vector
00170                                                                  // product; // lots of round off errors at stake
00171     if( extraPoints >= 1 )
00172     {
00173         for( int k = 0; k < nsTgt; k++ )
00174         {
00175             if( side[k] )
00176             {
00177                 // this is to mark that tgt edges k-1 and k are intersecting src
00178                 markr[k] = 1;
00179                 markr[( k + nsTgt - 1 ) % nsTgt] =
00180                     1;  // it is the previous edge, actually, but instead of doing -1, it is
00181                 // better to do modulo +3 (modulo 4)
00182                 // null side b for next call
00183             }
00184         }
00185     }
00186     nP += extraPoints;
00187 
00188     // now sort and orient the points in P, such that they are forming a convex polygon
00189     // this will be the foundation of our new mesh
00190     // this works if the polygons are convex
00191     IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 );  // nP should be at most 8 in the end ?
00192     // if there are more than 3 points, some area will be positive
00193 
00194     if( nP >= 3 )
00195     {
00196         for( int k = 1; k < nP - 1; k++ )
00197             area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
00198     }
00199 
00200     return MB_SUCCESS;  // no error
00201 }
00202 
00203 // this method will also construct the triangles/polygons in the new mesh
00204 // if we accept planar polygons, we just save them
00205 // also, we could just create new vertices every time, and merge only in the end;
00206 // could be too expensive, and the tolerance for merging could be an
00207 // interesting topic
00208 ErrorCode IntxRllCssphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP )
00209 {
00210     // first of all, check against tgt and src vertices
00211     //
00212 #ifdef ENABLE_DEBUG
00213     if( dbg_1 )
00214     {
00215         std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
00216                   << "\n";
00217         for( int n = 0; n < nP; n++ )
00218             std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
00219     }
00220 #endif
00221 
00222     // get the edges for the tgt triangle; the extra points will be on those edges, saved as
00223     // lists (unordered)
00224 
00225     // first get the list of edges adjacent to the tgt cell
00226     // use the neighTgtEdgeTag
00227     EntityHandle adjTgtEdges[MAXEDGES];
00228     ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge tgt tag" );
00229     // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
00230     // potatoes
00231 
00232     // these will be in the new mesh, mbOut
00233     // some of them will be handles to the initial vertices from src or tgt meshes (lagr or euler)
00234 
00235     EntityHandle* foundIds = new EntityHandle[nP];
00236     for( int i = 0; i < nP; i++ )
00237     {
00238         double* pp = &iP[2 * i];  // iP+2*i
00239         // project the point back on the sphere
00240         CartVect pos;
00241         IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], R, plane, pos );
00242         int found = 0;
00243         // first, are they on vertices from tgt or src?
00244         // priority is the tgt mesh (mb2?)
00245         int j                = 0;
00246         EntityHandle outNode = (EntityHandle)0;
00247         for( j = 0; j < nsTgt && !found; j++ )
00248         {
00249             // int node = tgtTri.v[j];
00250             double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
00251             if( d2 < epsilon_1 )
00252             {
00253 
00254                 foundIds[i] = tgtConn[j];  // no new node
00255                 found       = 1;
00256 #ifdef ENABLE_DEBUG
00257                 if( dbg_1 )
00258                     std::cout << "  tgt node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
00259                               << " 2d coords:" << tgtCoords2D[2 * j] << "  " << tgtCoords2D[2 * j + 1] << " d2: " << d2
00260                               << " \n";
00261 #endif
00262             }
00263         }
00264 
00265         for( j = 0; j < nsSrc && !found; j++ )
00266         {
00267             // int node = srcTri.v[j];
00268             double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
00269             if( d2 < epsilon_1 )
00270             {
00271                 // suspect is srcConn[j] corresponding in mbOut
00272 
00273                 foundIds[i] = srcConn[j];  // no new node
00274                 found       = 1;
00275 #ifdef ENABLE_DEBUG
00276                 if( dbg_1 )
00277                     std::cout << "  src node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2 << " \n";
00278 #endif
00279             }
00280         }
00281         if( !found )
00282         {
00283             // find the edge it belongs, first, on the tgt element
00284             //
00285             for( j = 0; j < nsTgt; j++ )
00286             {
00287                 int j1      = ( j + 1 ) % nsTgt;
00288                 double area = IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp );
00289 #ifdef ENABLE_DEBUG
00290                 if( dbg_1 )
00291                     std::cout << "   edge " << j << ": " << mb->id_from_handle( adjTgtEdges[j] ) << " " << tgtConn[j]
00292                               << " " << tgtConn[j1] << "  area : " << area << "\n";
00293 #endif
00294                 if( fabs( area ) < epsilon_1 / 2 )
00295                 {
00296                     // found the edge; now find if there is a point in the list here
00297                     // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
00298                     int indx = TgtEdges.index( adjTgtEdges[j] );
00299                     // CID 181167 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
00300                     if( indx < 0 )
00301                     {
00302                         std::cerr << " error in adjacent tgt edge: " << mb->id_from_handle( adjTgtEdges[j] ) << "\n";
00303                         delete[] foundIds;
00304                         return MB_FAILURE;
00305                     }
00306                     std::vector< EntityHandle >* expts = extraNodesVec[indx];
00307                     // if the points pp is between extra points, then just give that id
00308                     // if not, create a new point, (check the id)
00309                     // get the coordinates of the extra points so far
00310                     int nbExtraNodesSoFar = expts->size();
00311                     if( nbExtraNodesSoFar > 0 )
00312                     {
00313                         CartVect* coords1 = new CartVect[nbExtraNodesSoFar];
00314                         mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
00315                         // std::list<int>::iterator it;
00316                         for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
00317                         {
00318                             // int pnt = *it;
00319                             double d2 = ( pos - coords1[k] ).length_squared();
00320                             if( d2 < epsilon_1 )
00321                             {
00322                                 found       = 1;
00323                                 foundIds[i] = ( *expts )[k];
00324 #ifdef ENABLE_DEBUG
00325                                 if( dbg_1 ) std::cout << " found node:" << foundIds[i] << std::endl;
00326 #endif
00327                             }
00328                         }
00329                         delete[] coords1;
00330                     }
00331                     if( !found )
00332                     {
00333                         // create a new point in 2d (at the intersection)
00334                         // foundIds[i] = m_num2dPoints;
00335                         // expts.push_back(m_num2dPoints);
00336                         // need to create a new node in mbOut
00337                         // this will be on the edge, and it will be added to the local list
00338                         mb->create_vertex( pos.array(), outNode );
00339                         ( *expts ).push_back( outNode );
00340                         foundIds[i] = outNode;
00341                         found       = 1;
00342 #ifdef ENABLE_DEBUG
00343                         if( dbg_1 ) std::cout << " new node: " << outNode << std::endl;
00344 #endif
00345                     }
00346                 }
00347             }
00348         }
00349         if( !found )
00350         {
00351             std::cout << " tgt quad: ";
00352             for( int j1 = 0; j1 < nsTgt; j1++ )
00353             {
00354                 std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
00355             }
00356             std::cout << " a point pp is not on a tgt quad " << *pp << " " << pp[1] << " tgt quad "
00357                       << mb->id_from_handle( tgt ) << " \n";
00358             delete[] foundIds;
00359             return MB_FAILURE;
00360         }
00361     }
00362 #ifdef ENABLE_DEBUG
00363     if( dbg_1 )
00364     {
00365         std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
00366         for( int i1 = 0; i1 < nP; i1++ )
00367             std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
00368     }
00369 #endif
00370     // first, find out if we have nodes collapsed; shrink them
00371     // we may have to reduce nP
00372     // it is possible that some nodes are collapsed after intersection only
00373     // nodes will always be in order (convex intersection)
00374     correct_polygon( foundIds, nP );
00375     // now we can build the triangles, from P array, with foundIds
00376     // we will put them in the out set
00377     if( nP >= 3 )
00378     {
00379         EntityHandle polyNew;
00380         mb->create_element( MBPOLYGON, foundIds, nP, polyNew );
00381         mb->add_entities( outSet, &polyNew, 1 );
00382 
00383         // tag it with the index ids from tgt and src sets
00384         int id = rs1.index( src );  // index starts from 0
00385         mb->tag_set_data( srcParentTag, &polyNew, 1, &id );
00386         id = rs2.index( tgt );
00387         mb->tag_set_data( tgtParentTag, &polyNew, 1, &id );
00388 
00389         counting++;
00390         mb->tag_set_data( countTag, &polyNew, 1, &counting );
00391 
00392 #ifdef ENABLE_DEBUG
00393         if( dbg_1 )
00394         {
00395 
00396             std::cout << "Counting: " << counting << "\n";
00397             std::cout << " polygon " << mb->id_from_handle( polyNew ) << "  nodes: " << nP << " :";
00398             for( int i1 = 0; i1 < nP; i1++ )
00399                 std::cout << " " << mb->id_from_handle( foundIds[i1] );
00400             std::cout << " plane: " << plane << "\n";
00401             std::vector< CartVect > posi( nP );
00402             mb->get_coords( foundIds, nP, &( posi[0][0] ) );
00403             for( int i1 = 0; i1 < nP; i1++ )
00404                 std::cout << foundIds[i1] << " " << posi[i1] << "\n";
00405 
00406             std::stringstream fff;
00407             fff << "file0" << counting << ".vtk";
00408             mb->write_mesh( fff.str().c_str(), &outSet, 1 );
00409         }
00410 #endif
00411     }
00412     // disable_debug();
00413     delete[] foundIds;
00414     foundIds = NULL;
00415     return MB_SUCCESS;
00416 }
00417 
00418 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines