From: abn Date: Fri, 10 Dec 2021 19:53:16 +0000 (+0100) Subject: WIP handling warped cells in buildSlice3D X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=fe8575b0e9b9d265db18f7f6a7fa3f2a195f6ba0;p=tools%2Fmedcoupling.git WIP handling warped cells in buildSlice3D --- diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 0762ba266..908f3ac49 100755 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -3846,12 +3846,14 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const dou cut3DCurve[*it]=-1; mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve); std::vector< std::pair > cut3DSurf(mDesc2->getNumberOfCells()); + std::map > 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 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 ret=MEDCouplingUMesh::New("Slice3D",2); @@ -3907,9 +3909,10 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3DSurf(const double *origin, const mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve); mcIdType ncellsSub=subMesh->getNumberOfCells(); std::vector< std::pair > cut3DSurf(ncellsSub); + std::map > 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 conn(DataArrayIdType::New()),connI(DataArrayIdType::New()),cellIds2(DataArrayIdType::New()); connI->pushBackSilent(0); conn->alloc(0,1); const mcIdType *nodal=subMesh->getNodalConnectivity()->getConstPointer(); @@ -3975,14 +3978,15 @@ MCAuto MEDCouplingUMesh::clipSingle3DCellByPlane(const double sameNbNodes=(mDesc1->getNumberOfNodes()==oldNbNodes); } std::vector< std::pair > cut3DSurf(mDesc2->getNumberOfCells()); + std::map > 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 conn(DataArrayIdType::New()),connI(DataArrayIdType::New()); connI->pushBackSilent(0); conn->alloc(0,1); { MCAuto 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 !"); } diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 2d5f408d7..c12504d63 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -340,9 +340,12 @@ namespace MEDCoupling std::vector& addCoordsQuadratic, std::vector& cr, std::vector& crI, std::vector& cNb1, std::vector& cNb2); static void AssemblyForSplitFrom3DCurve(const std::vector& cut3DCurve, std::vector& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf, const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve, - const mcIdType *desc, const mcIdType *descIndx, std::vector< std::pair >& cut3DSurf); + const mcIdType *desc, const mcIdType *descIndx, std::vector< std::pair >& cut3DSurf, + std::map > & warpMap); void assemblyForSplitFrom3DSurf(const std::vector< std::pair >& cut3DSurf, - const mcIdType *desc, const mcIdType *descIndx, DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const; + const mcIdType *desc, const mcIdType *descIndx, + const std::map > &warpMap, + DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const; void buildSubCellsFromCut(const std::vector< std::pair >& cut3DSurf, const mcIdType *desc, const mcIdType *descIndx, const double *coords, double eps, std::vector >& 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); diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx index 7bbe37f8e..5deca19a6 100755 --- a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx @@ -1605,7 +1605,8 @@ void MEDCouplingUMesh::AppendExtrudedCell(const mcIdType *connBg, const mcIdType void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& cut3DCurve, std::vector& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf, const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve, const mcIdType *desc, const mcIdType *descIndx, - std::vector< std::pair >& cut3DSurf) + std::vector< std::pair >& cut3DSurf, + std::map > & warpMap) { std::set nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end()); mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size()); @@ -1614,29 +1615,30 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& std::vector res; mcIdType offset=descIndx[i]; mcIdType nbOfSeg=descIndx[i+1]-offset; - for(mcIdType j=0;j-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 s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]); @@ -1645,20 +1647,28 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& { 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& nodesW = warpMap[i]; + nodesW.resize(res.size()); + std::copy(res.begin(), res.end(), nodesW.begin()); + } } } } @@ -1676,6 +1686,7 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& */ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair >& cut3DSurf, const mcIdType *desc, const mcIdType *descIndx, + const std::map > &warpMap, DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const { checkFullyDefined(); @@ -1695,13 +1706,35 @@ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair& 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 & 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; @@ -1719,6 +1752,7 @@ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON)); mcIdType prev=end; while(end!=start)