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];
mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
for(mcIdType i=0;i<nbOf3DSurfCell;i++)
{
- std::vector<mcIdType> res;
+ std::vector<mcIdType> res; // list of points that will be used to reconstruct intersection polygon
+ std::vector<bool> 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<nbOfSeg;j++) // for all segments of a face of the initial 3D mesh
+ // Loop: for all segments of a face of the initial 3D mesh, record either a (single) point of intersection with the given segment
+ // or a full segment if it is coplanar.
+ // We also should deal with limiting cases here: if we just added a full segment, the next edge can not have a single
+ // intersecting point (if it is so, it means that this point is actually very close to one of the extremity, and we should
+ // discard it).
+ unsigned lastPushed=-1; // 1 means single point, 2 means seg, 0 means nothing
+ bool lastWillBeSeg = cut3DCurve[nbOfSeg-1] > -1 ; // whether a full segment will be added when scanning the last edge
+ for(mcIdType j=0;j<nbOfSeg;j++)
{
mcIdType edgeId=desc[offset+j];
mcIdType status=cut3DCurve[edgeId];
if(status!=-2) // segment is being intersected (BUT a segment might be colinear with a -2)
{
if(status>-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<mcIdType> 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<mcIdType> >(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<mcIdType>& 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<mcIdType> 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<mcIdType> >(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<mcIdType>& nodesW = warpMap[i];
+ bool inSegs = false;
+ std::vector<mcIdType> 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;
+ }
}
}
mcIdType nbOfCells=getNumberOfCells();
for(mcIdType i=0;i<nbOfCells;i++)
{
- std::map<mcIdType, std::set<mcIdType> > m;
+ std::map<mcIdType, std::set<mcIdType> > 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<mcIdType> *vecWarp = nullptr;
for(int j=0;j<nbOfFaces;j++)
{
const std::pair<mcIdType,mcIdType>& 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<mcIdType> & 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???)
{
if(m.empty())
continue;
// Now rebuild new 2D cell connectivity corresponding to the currently intersected 3D cell:
- std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
- mcIdType prev=end;
- while(end!=start)
+ // Classical (non-warped) case:
+ if (!vecWarp)
{
- std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
- const std::set<mcIdType>& s=(*it).second;
- std::set<mcIdType> s2; s2.insert(prev);
- std::set<mcIdType> s3;
- std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
- if(s3.size()==1)
+ std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
+ mcIdType prev=end;
+ while(end!=start)
+ {
+ std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
+ const std::set<mcIdType>& s=(*it).second;
+ std::set<mcIdType> s2; s2.insert(prev);
+ std::set<mcIdType> 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<mcIdType> consumed;
+ mcIdType kk=0;
+ for(auto it=m.begin(); it != m.end(); it++, kk++)
+ {
+ const std::set<mcIdType>& 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<mcIdType> 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<nbOfCells ...
}
