From: abn Date: Thu, 16 Dec 2021 14:30:25 +0000 (+0100) Subject: WIP warped faces ... still not working properly ... X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=eaffffee3b032a692e07f8a78e8dcc3f1a30dccc;p=tools%2Fmedcoupling.git WIP warped faces ... still not working properly ... --- diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx index 5deca19a6..5cafb7e47 100755 --- a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx @@ -581,7 +581,7 @@ void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double const double *st2=coo+3*st; vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2]; double pos=-(vec4[0]*vec2[0]+vec4[1]*vec2[1]+vec4[2]*vec2[2])/((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])); - if(pos>eps && pos<1-eps) + if(pos>eps && pos<1-eps) // intersecting and creating a new point sufficiently far away from the extremities { mcIdType nNode=ToIdType(addCoo.size())/3; vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2]; @@ -1612,65 +1612,115 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size()); for(mcIdType i=0;i res; + std::vector res; // list of points that will be used to reconstruct intersection polygon + std::vector resIsSeg; // same size as res, true if the point is part of a segment + mcIdType segCnt = 0; mcIdType offset=descIndx[i]; mcIdType nbOfSeg=descIndx[i+1]-offset; - for(mcIdType j=0;j -1 ; // whether a full segment will be added when scanning the last edge + for(mcIdType j=0;j-1) // ... just once, hence using one of the newly created point - res.push_back(status); - else // twice: the segment is inside the plane. Take two existing points { - res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]); - res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]); + // if we pushed a seg just before, this intersection point is very close to one extremity + // "just before" can also mean the last edge in the face, so we need to look this too ... + if (lastPushed != 2 && !lastWillBeSeg) + { res.push_back(status); resIsSeg.push_back(false); + lastPushed = 1; } + else + lastPushed = 0; + } + else // status==-1 : two intersec: the segment is inside the plane. Take two existing points. + { + if (lastPushed == 1) // Remove last added point as it is very close to extremity + { res.pop_back(); resIsSeg.pop_back(); } + res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]); resIsSeg.push_back(true); + res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]); resIsSeg.push_back(true); + lastPushed = 2; + segCnt++; } } + else + lastPushed = 0; } + // At the end of the loop, res might have been filled several times even for a single face. - switch(res.size()) - { - case 2: // regular case, two intersection points - { - cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1]; - break; - } - case 1: // 1 and 0 means colinear or surface cell touching plane by exact corner points - case 0: - { - std::set s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]); - std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector >(res)); - if(res.size()==2) - { - cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1]; - } - else // only one touching point (or none) - nothing recorded - { - cut3DSurf[i].first=-1; cut3DSurf[i].second=-1; - } - break; - } - default: - {// case when the face is completly included in the plane - if(ToIdType(res.size())==2*nbOfSeg) - { - cut3DSurf[i].first=-2; cut3DSurf[i].second=i; - } - else { - cut3DSurf[i].first = -3; cut3DSurf[i].second=i; - // Fills in warpMap: it is important for the rest of the algo that the nodes are filled in the proper order - // i.e. in the same order as the edges around the face. Luckily for us, we have this for free: - // the for loop above on the edges in the right order because the descending connectivity of a mesh already ensures - // the proper ordering of the segments in desc/descI!! - std::vector& nodesW = warpMap[i]; - nodesW.resize(res.size()); - std::copy(res.begin(), res.end(), nodesW.begin()); + if(res.size() == 2) + { // regular case, two intersection points + cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1]; + } + else if (res.size() <= 1 && segCnt == 0) + { // one segment included in the plane, or just one point touching + std::set s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]); + std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(), + std::back_insert_iterator< std::vector >(res)); + if(res.size()==2) + { cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1]; } + else // only one touching point (or none) - nothing recorded + { cut3DSurf[i].first=-1; cut3DSurf[i].second=-1; } + } + else if(segCnt == nbOfSeg) + { // exact inclusion of a (clean plane) face into the plane + cut3DSurf[i].first=-2; cut3DSurf[i].second=i; + if(ToIdType(res.size()) != 2*nbOfSeg) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyForSplitFrom3DCurve: woopsy, unexpected situation..."); + } + else if (res.size() == 4 || segCnt >=2) // warped case!! Only support warped QUAD4 for now + { + cut3DSurf[i].first = -3; cut3DSurf[i].second=i; + // Rebuild the proper connectivity of the central polygon that will be inserted for this face + + // With QUAD4, the coplanar segments are necessarily consecutive (even if warped) + // Fills in warpMap: it is important for the rest of the algo that the nodes are filled in the proper order + // i.e. in the same order as the edges around the face. Luckily for us, we have this for free: + // the for loop above on the edges in the right order because the descending connectivity of a mesh already ensures + // the proper ordering of the segments in desc/descI!! + std::vector& nodesW = warpMap[i]; + bool inSegs = false; + std::vector segs; + for (mcIdType kk = 0; kk <= ToIdType(res.size()); kk++) // yes one slot more to cover last iteration + { + if (inSegs && (kk==ToIdType(res.size()) || !resIsSeg[kk])) + { // Unstack properly previously recorded set of (consecutive) segments to re-order points correctly + if (segs[1]==segs[2] || segs[1]==segs[3]) + { nodesW.push_back(segs[0]); + nodesW.push_back(segs[1]); + } + else { + nodesW.push_back(segs[1]); + nodesW.push_back(segs[0]); + } + if(nodesW[1] == segs[2]) nodesW.push_back(segs[3]); + else nodesW.push_back(segs[2]); + + for (mcIdType jj = 2; jj < ToIdType(segs.size()/2); jj++) + if (segs[jj*2] != segs[(jj-1)*2+1]) + nodesW.push_back(segs[jj*2+1]); + else + nodesW.push_back(segs[jj*2]); + inSegs = false; + } + if (kk==ToIdType(res.size())) break; // last iteration, we should stop here + if (!resIsSeg[kk]) // regular point (not part of a segment) + nodesW.push_back(res[kk]); + else + { inSegs = true; segs.push_back(res[kk]); } } - } - } + } + else // face close to a warp face might lead to funny things - discard + { + cut3DSurf[i].first = -4; cut3DSurf[i].second=i; + } } } @@ -1696,43 +1746,36 @@ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair > m; + std::map > m; // map giving for each vertex of the reconstructed 2D cells, the set set of points + // connected to it - said differently this describes the edges of the polygon to build mcIdType offset=descIndx[i]; mcIdType nbOfFaces=descIndx[i+1]-offset; mcIdType start=-1; mcIdType end=-1; + const std::vector *vecWarp = nullptr; for(int j=0;j& p=cut3DSurf[desc[offset+j]]; if(p.first!=-1 && p.second!=-1) { - if(p.first!=-2 && p.first!=-3) // not a coplanar polygon, not a warped face + if(p.first > -2) // not a coplanar polygon (-2), not a warped face (-3), not an ignored case (-4) { start=p.first; end=p.second; m[p.first].insert(p.second); m[p.second].insert(p.first); } else if (p.first == -3) - { // Warped face ... for now we only handle warped QUADs ... -// const std::vector & vec = warpMap.at(offset+j); -// if (vec.size() != 4) -// throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf: only warped QUAD4 are properly handled! Looks like we are hitting a more complex case ... sorry."); -// start = vec.front(); end = vec.back(); -// mcIdType k=0, sz=vec.size(); -// // Handle this as if we had two different faces: insert (0,1), (1,0) in the map -// // and then (2,3) and (3,2). Nothing more (not (1,2) ...) -// m[vec[0]].insert(vec[1]); -// m[vec[1]].insert(vec[0]); -// m[vec[2]].insert(vec[3]); -// m[vec[3]].insert(vec[2]); -// for (auto it = vec.begin(); it != vec.end(); it++, k++) -// { -// mcIdType id1=*it, id2=vec[(k+1)%sz]; -// m[id1].insert(id2); -// m[id2].insert(id1); -// } + { // Warped face ... for now we only handle warped QUADs, and only one warped face per 3D cell + if (vecWarp) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf: more than one warped face in a cell!! Not supported."); + vecWarp = &warpMap.at(p.second); } - else // face is coplanar to the plane, copy it entirely. TODO: here we do this + else if (p.first == -4) + { + // Do nothing ... face close to another warped face + } + else // p.first == -2 + // face is coplanar to the plane, copy it entirely. TODO: here we do this // by extracting its descending connectivity so that the reconstruction algo (10 lines down below) works the same. // this could probably be done more efficiently! (convertToPoly and copy connectivity???) { @@ -1753,33 +1796,73 @@ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON)); - mcIdType prev=end; - while(end!=start) + // Classical (non-warped) case: + if (!vecWarp) { - std::map >::const_iterator it=m.find(start); - const std::set& s=(*it).second; - std::set s2; s2.insert(prev); - std::set s3; - std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin())); - if(s3.size()==1) + std::vector conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON)); + mcIdType prev=end; + while(end!=start) + { + std::map >::const_iterator it=m.find(start); + const std::set& s=(*it).second; + std::set s2; s2.insert(prev); + std::set s3; + std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin())); + if(s3.size()==1) + { + mcIdType val=*s3.begin(); + conn.push_back(start); + prev=start; + start=val; + } + else + start=end; + } + conn.push_back(end); + if(conn.size()>3) { - mcIdType val=*s3.begin(); - conn.push_back(start); - prev=start; - start=val; + nodalRes->insertAtTheEnd(conn.begin(),conn.end()); + nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples()); + cellIds->pushBackSilent(i); } - else - start=end; } - conn.push_back(end); - if(conn.size()>3) + else // Cell with a warped face { - nodalRes->insertAtTheEnd(conn.begin(),conn.end()); + // When a plane intersects a warped face, we have (in the worst case): + // - two (disjoint) triangles + // - connected by a single quad + // We start with the two triangles by identifying the points with a set containing two items + std::set consumed; + mcIdType kk=0; + for(auto it=m.begin(); it != m.end(); it++, kk++) + { + const std::set& s = (*it).second; + const mcIdType nodId = (*it).first; + if (consumed.find(nodId) != consumed.end()) continue; + if(s.size() == 2) + { + const mcIdType nod1=*s.begin(), nod2=*s.rbegin(); + const auto &s1 = m.at(nod1), &s2 = m.at(nod2); + if (s1.size() != 1 || s2.size() != 1 ) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf: unexpected situation when building intersection with warped face ..."); + nodalRes->pushBackSilent(ToIdType(INTERP_KERNEL::NORM_POLYGON)); + nodalRes->pushBackSilent(nod1); consumed.insert(nod1); + nodalRes->pushBackSilent(nodId); consumed.insert(nodId); + nodalRes->pushBackSilent(nod2); consumed.insert(nod2); + nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples()); + cellIds->pushBackSilent(i); + } + } + // Finally insert the QUAD4 (but as a polygon) corresponding to the central piece + if (vecWarp->size() < 4) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf: should not happen?"); + std::vector connQuad(vecWarp->size()+1,ToIdType(INTERP_KERNEL::NORM_POLYGON)); + std::copy(vecWarp->begin(), vecWarp->end(), connQuad.begin()+1); // +1 to skip type + nodalRes->insertAtTheEnd(connQuad.begin(), connQuad.end()); nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples()); cellIds->pushBackSilent(i); } - } + } // for(mcIdType i=0;i