cut3DCurve[*it]=-1;
mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
std::vector< std::pair<mcIdType,mcIdType> > cut3DSurf(mDesc2->getNumberOfCells());
+ std::map<mcIdType, std::vector<mcIdType> > warpMap; // dealing with warp cells ...
AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,mDesc2->getNodalConnectivity()->getConstPointer(),mDesc2->getNodalConnectivityIndex()->getConstPointer(),
mDesc1->getNodalConnectivity()->getConstPointer(),mDesc1->getNodalConnectivityIndex()->getConstPointer(),
- desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf);
+ desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf,
+ warpMap);
MCAuto<DataArrayIdType> conn(DataArrayIdType::New()),connI(DataArrayIdType::New()),cellIds2(DataArrayIdType::New());
connI->pushBackSilent(0); conn->alloc(0,1); cellIds2->alloc(0,1);
- subMesh->assemblyForSplitFrom3DSurf(cut3DSurf,desc2->getConstPointer(),descIndx2->getConstPointer(),conn,connI,cellIds2);
+ subMesh->assemblyForSplitFrom3DSurf(cut3DSurf,desc2->getConstPointer(),descIndx2->getConstPointer(),warpMap, conn,connI,cellIds2);
if(cellIds2->empty())
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane !");
MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Slice3D",2);
mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
mcIdType ncellsSub=subMesh->getNumberOfCells();
std::vector< std::pair<mcIdType,mcIdType> > cut3DSurf(ncellsSub);
+ std::map<mcIdType, std::vector<mcIdType> > dnu;
AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,subMesh->getNodalConnectivity()->getConstPointer(),subMesh->getNodalConnectivityIndex()->getConstPointer(),
mDesc1->getNodalConnectivity()->getConstPointer(),mDesc1->getNodalConnectivityIndex()->getConstPointer(),
- desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf);
+ desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf, dnu);
MCAuto<DataArrayIdType> conn(DataArrayIdType::New()),connI(DataArrayIdType::New()),cellIds2(DataArrayIdType::New()); connI->pushBackSilent(0);
conn->alloc(0,1);
const mcIdType *nodal=subMesh->getNodalConnectivity()->getConstPointer();
sameNbNodes=(mDesc1->getNumberOfNodes()==oldNbNodes);
}
std::vector< std::pair<mcIdType,mcIdType> > cut3DSurf(mDesc2->getNumberOfCells());
+ std::map<mcIdType, std::vector<mcIdType> > dnu;
AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,mDesc2->getNodalConnectivity()->begin(),mDesc2->getNodalConnectivityIndex()->begin(),
mDesc1->getNodalConnectivity()->begin(),mDesc1->getNodalConnectivityIndex()->begin(),
- desc1->begin(),descIndx1->begin(),cut3DSurf);
+ desc1->begin(),descIndx1->begin(),cut3DSurf, dnu);
MCAuto<DataArrayIdType> conn(DataArrayIdType::New()),connI(DataArrayIdType::New());
connI->pushBackSilent(0); conn->alloc(0,1);
{
MCAuto<DataArrayIdType> cellIds2(DataArrayIdType::New()); cellIds2->alloc(0,1);
- assemblyForSplitFrom3DSurf(cut3DSurf,desc2->begin(),descIndx2->begin(),conn,connI,cellIds2);
+ assemblyForSplitFrom3DSurf(cut3DSurf,desc2->begin(),descIndx2->begin(),dnu, conn,connI,cellIds2);
if(cellIds2->empty())
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane !");
}
std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& cr, std::vector<mcIdType>& crI, std::vector<mcIdType>& cNb1, std::vector<mcIdType>& cNb2);
static void AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
- const mcIdType *desc, const mcIdType *descIndx, std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf);
+ const mcIdType *desc, const mcIdType *descIndx, std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
+ std::map<mcIdType, std::vector<mcIdType> > & warpMap);
void assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
- const mcIdType *desc, const mcIdType *descIndx, DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const;
+ const mcIdType *desc, const mcIdType *descIndx,
+ const std::map<mcIdType, std::vector<mcIdType> > &warpMap,
+ DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const;
void buildSubCellsFromCut(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf, const mcIdType *desc, const mcIdType *descIndx, const double *coords, double eps, std::vector<std::vector<mcIdType> >& res) const;
void split2DCellsLinear(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI);
mcIdType split2DCellsQuadratic(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI, const DataArrayIdType *mid, const DataArrayIdType *midI);
void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
const mcIdType *desc, const mcIdType *descIndx,
- std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf)
+ std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
+ std::map<mcIdType, std::vector<mcIdType> > & warpMap)
{
std::set<mcIdType> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
std::vector<mcIdType> res;
mcIdType offset=descIndx[i];
mcIdType nbOfSeg=descIndx[i+1]-offset;
- for(mcIdType j=0;j<nbOfSeg;j++)
+ for(mcIdType j=0;j<nbOfSeg;j++) // for all segments of a face of the initial 3D mesh
{
mcIdType edgeId=desc[offset+j];
mcIdType status=cut3DCurve[edgeId];
- if(status!=-2)
+ if(status!=-2) // segment is being intersected (BUT a segment might be colinear with a -2)
{
- if(status>-1)
+ if(status>-1) // ... just once, hence using one of the newly created point
res.push_back(status);
- else
+ 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]);
}
}
}
+ // At the end of the loop, res might have been filled several times even for a single face.
switch(res.size())
{
- case 2:
+ case 2: // regular case, two intersection points
{
cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
break;
}
- case 1:
+ 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]);
{
cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
}
- else
+ else // only one touching point (or none) - nothing recorded
{
cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
}
break;
}
default:
- {// case when plane is on a multi colinear edge of a polyhedron
+ {// 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
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
+ 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());
+ }
}
}
}
*/
void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
const mcIdType *desc, const mcIdType *descIndx,
+ const std::map<mcIdType, std::vector<mcIdType> > &warpMap,
DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const
{
checkFullyDefined();
const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
if(p.first!=-1 && p.second!=-1)
{
- if(p.first!=-2)
+ if(p.first!=-2 && p.first!=-3) // not a coplanar polygon, not a warped face
{
start=p.first; end=p.second;
m[p.first].insert(p.second);
m[p.second].insert(p.first);
}
- else
+ 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);
+// }
+ }
+ else // 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???)
{
const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
}
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)
