moab
|
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