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