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