X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh_intersection.cxx;h=b233f1dcb7ad980485d086f9edabbe54671b3a25;hb=3e9ebc98a6bc02022900a8b2295e1a775d78ece0;hp=8fdf018821e362b7ef3fef7985181fd6fb364e79;hpb=14f30db13a9749dd4a47f5c14944c7331e7390d1;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 8fdf01882..b233f1dcb 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -20,6 +20,7 @@ #include "MEDCouplingUMesh.hxx" #include "MEDCoupling1GTUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "CellModel.hxx" #include "VolSurfUser.txx" #include "InterpolationUtils.hxx" @@ -33,6 +34,9 @@ #include "InterpKernelGeo2DEdgeLin.hxx" #include "InterpKernelGeo2DEdgeArcCircle.hxx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" +#include "TranslationRotationMatrix.hxx" +#include "VectorUtils.hxx" +#include "MEDCouplingSkyLineArray.hxx" #include #include @@ -40,6 +44,7 @@ #include #include #include +#include using namespace MEDCoupling; @@ -277,7 +282,7 @@ namespace MEDCoupling } /** - * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI). + * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI). */ void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo, const int *desc1Bg, const int *desc1End, const std::vector >& intesctEdges1, @@ -332,7 +337,7 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg // This initializes posBaseElt. if(nbOfTurn==0) { - for(unsigned i=1;i& edges, const std::vector< MCAuto >& edgesPtr); std::size_t size() const { return _pool.size(); } int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const; - void setMeshAt(int pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs); + void setMeshAt(std::size_t pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs); const std::vector& getConnOf(int pos) const { return get(pos)._edges; } const std::vector< MCAuto >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; } MCAuto getZeMesh() const { return _ze_mesh; } @@ -895,7 +900,7 @@ int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) co return zeMesh->getCellContainingPoint(barys->begin(),eps); } -void VectorOfCellInfo::setMeshAt(int pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs) +void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs) { get(pos);//to check pos bool isFast(pos==0 && _pool.size()==1); @@ -907,7 +912,7 @@ void VectorOfCellInfo::setMeshAt(int pos, const MCAuto& mesh, _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back())); // std::vector pool(_pool.size()-1+sz); - for(int i=0;i >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, std::vector& addCoo, std::map& mergedNodes) { static const int SPACEDIM=2; - INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin()); // Build BB tree of all edges in the tool mesh (second mesh) - MCAuto bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree()); + MCAuto bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps)); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells()); intersectEdge1.resize(nDescCell1); @@ -1264,7 +1269,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin()); int offset2(offset1+m2->getNumberOfNodes()); int offset3(offset2+((int)addCoords.size())/2); - MCAuto bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree()); + MCAuto bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps)); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); // Here a BBTree on 2D-cells, not on segments: BBTree myTree(bbox2,0,0,m2->getNumberOfCells(),eps); @@ -1395,7 +1400,7 @@ void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair& p=cut3DSurf[desc[offset+j]]; @@ -1426,7 +1431,7 @@ void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::paircheckFullyDefined(); m2->checkFullyDefined(); + INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!"); @@ -1612,7 +1619,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 /*! * Partitions the first given 2D mesh using the second given 1D mesh as a tool. * Thus the final result contains the aggregation of nodes of \a mesh2D, then nodes of \a mesh1D, then new nodes that are the result of the intersection - * and finaly, in case of quadratic polygon the centers of edges new nodes. + * and finally, in case of quadratic polygon the centers of edges new nodes. * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input. * * \param [in] mesh2D - the 2D mesh (spacedim=meshdim=2) to be intersected using \a mesh1D tool. The mesh must be so that each point in the space covered by \a mesh2D @@ -1642,8 +1649,8 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, // Step 1: compute all edge intersections (new nodes) std::vector< std::vector > intersectEdge1, colinear2, subDiv2; std::vector addCoo,addCoordsQuadratic; // coordinates of newly created nodes - INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); // // Build desc connectivity DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New()); @@ -1697,7 +1704,9 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MCAuto baryRet1(ret1NonCol->computeCellCenterOfMass()); MCAuto elts,eltsIndex; mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex); - MCAuto eltsIndex2(eltsIndex->deltaShiftIndex()); + MCAuto eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1); + if (eltsIndex->getNumberOfTuples() > 1) + eltsIndex2 = eltsIndex->deltaShiftIndex(); MCAuto eltsIndex3(eltsIndex2->findIdsEqual(1)); if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells()) throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !"); @@ -1770,7 +1779,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MCAuto splitOfOneCell(BuildMesh2DCutFrom(eps,*it,m1Desc,partOfMesh1CuttingCur2DCell,dd1->begin()+dd2->getIJ(*it,0),dd1->begin()+dd2->getIJ((*it)+1,0),intersectEdge1,ret2->getNumberOfTuples(),partOfRet3)); ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true); outMesh2DSplit.push_back(splitOfOneCell); - for(int i=0;igetNumberOfCells();i++) + for(std::size_t i=0;igetNumberOfCells();i++) ret2->pushBackSilent(*it); } // @@ -1783,7 +1792,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MCAuto ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp)); // To finish - filter ret3 - std::numeric_limits::max() -> -1 - negate values must be resolved. ret3->rearrange(1); - MCAuto edgesToDealWith(ret3->findIdsStricltyNegative()); + MCAuto edgesToDealWith(ret3->findIdsStrictlyNegative()); for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++) { int old2DCellId(-ret3->getIJ(*it,0)-1); @@ -1831,14 +1840,14 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) MCAuto desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New()); MCAuto mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1)); const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin()); - MCAuto bboxArr(mDesc->getBoundingBoxForBBTree()); + MCAuto bboxArr(mDesc->getBoundingBoxForBBTree(eps)); const double *bbox(bboxArr->begin()),*coords(getCoords()->begin()); int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells()); std::vector< std::vector > intersectEdge(nDescCell),overlapEdge(nDescCell); std::vector addCoo; BBTree myTree(bbox,0,0,nDescCell,-eps); - INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); for(int i=0;i candidates; @@ -1965,8 +1974,8 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) checkConsistencyLight(); if(getSpaceDimension()!=2 || getMeshDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !"); - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; - INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; + INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin()); MCAuto newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0); @@ -1994,4 +2003,461 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) return ret.retn(); } +///@cond INTERNAL +/** + * c, cI describe a wire mesh in 3D space, in local numbering + * startNode, endNode in global numbering + *\return true if the segment is indeed split + */ +bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode, + const int * c, const int * cI, const int *idsBg, const int *endBg, + std::vector & pointIds, std::vector & hitSegs) +{ + using namespace std; + + const int SPACEDIM=3; + typedef pair PairDI; + set< PairDI > x; + for (const int * it = idsBg; it != endBg; ++it) + { + assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2); + int start = c[cI[*it]+1], end = c[cI[*it]+2]; + x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords + x.insert(make_pair(coo[end*SPACEDIM], end)); + } + + vector xx(x.begin(), x.end()); + sort(xx.begin(),xx.end()); + pointIds.reserve(xx.size()); + + // Keep what is inside [startNode, endNode]: + int go = 0; + for (vector::const_iterator it=xx.begin(); it != xx.end(); ++it) + { + const int idx = (*it).second; + if (!go) + { + if (idx == startNode) go = 1; + if (idx == endNode) go = 2; + if (go) pointIds.push_back(idx); + continue; + } + pointIds.push_back(idx); + if (idx == endNode || idx == startNode) + break; + } + +// vector pointIds2(pointIds.size()+2); +// copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1); +// pointIds2[0] = startNode; +// pointIds2[pointIds2.size()-1] = endNode; + + if (go == 2) + reverse(pointIds.begin(), pointIds.end()); + + // Now identify smaller segments that are not sub-divided - those won't need any further treatment: + for (const int * it = idsBg; it != endBg; ++it) + { + int start = c[cI[*it]+1], end = c[cI[*it]+2]; + vector::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start); + if (itStart == pointIds.end()) continue; + vector::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end); + if (itEnd == pointIds.end()) continue; + if (abs(distance(itEnd, itStart)) != 1) continue; + hitSegs.push_back(*it); // segment is undivided. + } + + return (pointIds.size() > 2); // something else apart start and end node +} + +void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode, + const std::vector& insidePoints, std::vector& modifiedFace) +{ + using namespace std; + int dst = distance(sIdxConn, sIdxConnE); + modifiedFace.reserve(dst + insidePoints.size()-2); + modifiedFace.resize(dst); + copy(sIdxConn, sIdxConnE, modifiedFace.data()); + + vector::iterator shortEnd = modifiedFace.begin()+dst; + vector::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode); + if (startPos == shortEnd) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!"); + vector::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode); + if (endPos == shortEnd) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!"); + int d = distance(startPos, endPos); + if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ... + modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted. + else + modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend()); +} + +///@endcond + + +/*! + * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty). + * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) ! + * This method performs a conformization of \b this. + * + * Only polyhedron cells are supported. You can call convertAllToPoly() + * + * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too. + * This method expects that all nodes in \a this are not closer than \a eps. + * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method. + * + * \param [in] eps the relative error to detect merged edges. + * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array + * that the user is expected to deal with. + * + * \throw If \a this is not coherent. + * \throw If \a this has not spaceDim equal to 3. + * \throw If \a this has not meshDim equal to 3. + * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly() + */ +DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) +{ + using namespace std; + + static const int SPACEDIM=3; + checkConsistencyLight(); + if(getSpaceDimension()!=3 || getMeshDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!"); + if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first."); + + MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); + const double * coo(_coords->begin()); + MCAuto ret(DataArrayInt::New()); + + { + /************************* + * STEP 1 -- faces + *************************/ + MCAuto descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + MCAuto mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI)); + const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); + MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity())); + + // Build BBTree + MCAuto bboxArr(mDesc->getBoundingBoxForBBTree(eps)); + const double *bbox(bboxArr->begin()); getCoords()->begin(); + int nDescCell(mDesc->getNumberOfCells()); + BBTree myTree(bbox,0,0,nDescCell,-eps); + // Surfaces - handle biggest first + MCAuto surfF = mDesc->getMeasureField(true); + DataArrayDouble * surfs = surfF->getArray(); + // Normal field + MCAuto normalsF = mDesc->buildOrthogonalField(); + DataArrayDouble * normals = normalsF->getArray(); + const double * normalsP = normals->getConstPointer(); + + // Sort faces by decreasing surface: + vector< pair > S; + for(std::size_t i=0;i < surfs->getNumberOfTuples();i++) + { + pair p = make_pair(surfs->begin()[i], i); + S.push_back(p); + } + sort(S.rbegin(),S.rend()); // reverse sort + vector hit(nDescCell); + fill(hit.begin(), hit.end(), false); + vector hitPoly; // the final result: which 3D cells have been modified. + + for( vector >::const_iterator it = S.begin(); it != S.end(); it++) + { + int faceIdx = (*it).second; + if (hit[faceIdx]) continue; + + vector candidates, cands2; + myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates); + // Keep only candidates whose normal matches the normal of current face + for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) + { + bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps); + if (*it2 != faceIdx && col) + cands2.push_back(*it2); + } + if (!cands2.size()) + continue; + + // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later + INTERP_KERNEL::TranslationRotationMatrix rotation; + INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]), + coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]), + coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation); + + MCAuto mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called + MCAuto mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called + double * cooPartRef(mPartRef->_coords->getPointer()); + double * cooPartCand(mPartCand->_coords->getPointer()); + for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartCand+SPACEDIM*ii); + + // Localize faces in 2D thanks to barycenters + MCAuto baryPart = mPartCand->computeCellCenterOfMass(); + vector compo; compo.push_back(2); + MCAuto baryPartZ = baryPart->keepSelectedComponents(compo); + MCAuto idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps); + if (!idsGoodPlane->getNumberOfTuples()) + continue; + + baryPart = baryPart->selectByTupleId(*idsGoodPlane); + + compo[0] = 0; compo.push_back(1); + MCAuto baryPartXY = baryPart->keepSelectedComponents(compo); + mPartRef->changeSpaceDimension(2,0.0); + MCAuto cc(DataArrayInt::New()), ccI(DataArrayInt::New()); + mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI); + + if (!cc->getNumberOfTuples()) + continue; + MCAuto dsi = ccI->deltaShiftIndex(); + + { + MCAuto tmp = dsi->findIdsInRange(0, 2); + if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) + { + ostringstream oss; + oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!"; + throw INTERP_KERNEL::Exception(oss.str()); + } + } + + MCAuto ids = dsi->findIdsEqual(1); + // Boundary face: + if (!ids->getNumberOfTuples()) + continue; + + double checkSurf=0.0; + const int * idsGoodPlaneP(idsGoodPlane->begin()); + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + { + int faceIdx2 = cands2[idsGoodPlaneP[*ii]]; + hit[faceIdx2] = true; + checkSurf += surfs->begin()[faceIdx2]; + } + if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps) + { + ostringstream oss; + oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!"; + throw INTERP_KERNEL::Exception(oss.str()); + } + + // For all polyhedrons using this face, replace connectivity: + vector polyIndices, packsIds, facePack; + for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++) + polyIndices.push_back(revDescP[ii]); + ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size()); + + // Current face connectivity + const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1; + const int * sIdxConnE = cDesc + cIDesc[faceIdx+1]; + connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); + // Deletion of old faces + int jj=0; + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj) + { + if (packsIds[jj] == -1) + // The below should never happen - if a face is used several times, with a different layout of the nodes + // it means that is is already conform, so it is *not* hit by the algorithm. The algorithm only hits + // faces which are actually used only once, by a single cell. This is different for edges below. + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error."); + else + connSla->deletePack(*it2, packsIds[jj]); + } + // Insertion of new faces: + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + { + // Build pack from the face to insert: + int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)]; + int facePack2Sz; + const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type! + // Insert it in all hit polyhedrons: + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2) + connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type + } + } + } // end step1 + + // Set back modified connectivity + MCAuto cAuto; cAuto.takeRef(_nodal_connec); + MCAuto cIAuto; cIAuto.takeRef(_nodal_connec_index); + connSla->convertToPolyhedronConn(cAuto, cIAuto); + + { + /************************ + * STEP 2 -- edges + ************************/ + // Now we have a face-conform mesh. + + // Recompute descending + MCAuto desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + // Rebuild desc connectivity with orientation this time!! + MCAuto mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI)); + const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer()); + const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); + MCAuto ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy()); + MCAuto cDeepC(mDesc->getNodalConnectivity()->deepCopy()); + MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC)); + MCAuto desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New()); + MCAuto mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2); +// std::cout << "writing!\n"; +// mDesc->writeVTK("/tmp/toto_desc_confInter.vtu"); +// mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu"); + const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); + const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); + MCAuto bboxArr(mDesc2->getBoundingBoxForBBTree(eps)); + const double *bbox2(bboxArr->begin()); + int nDesc2Cell=mDesc2->getNumberOfCells(); + BBTree myTree2(bbox2,0,0,nDesc2Cell,-eps); + + // Edges - handle longest first + MCAuto lenF = mDesc2->getMeasureField(true); + DataArrayDouble * lens = lenF->getArray(); + + // Sort edges by decreasing length: + vector > S; + for(std::size_t i=0;i < lens->getNumberOfTuples();i++) + { + pair p = make_pair(lens->getIJ(i, 0), i); + S.push_back(p); + } + sort(S.rbegin(),S.rend()); // reverse sort + + vector hit(nDesc2Cell); + fill(hit.begin(), hit.end(), false); + + for( vector >::const_iterator it = S.begin(); it != S.end(); it++) + { + int eIdx = (*it).second; + if (hit[eIdx]) + continue; + + vector candidates, cands2; + myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates); + // Keep only candidates colinear with current edge + double vCurr[3]; + unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; + for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? + vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3]; + for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) + { + double vOther[3]; + unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; + for (int i3=0; i3 < 3; i3++) + vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3]; + bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps); + // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need + // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine()) + if (col) + cands2.push_back(*it2); + } + if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above + continue; + + // Now rotate edges to bring them on Ox + int startNode = cDesc2[cIDesc2[eIdx]+1]; + int endNode = cDesc2[cIDesc2[eIdx]+2]; + INTERP_KERNEL::TranslationRotationMatrix rotation; + INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation); + MCAuto mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called + MCAuto mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called + MCAuto nodeMap(mPartCand->zipCoordsTraducer()); + int nbElemsNotM1; + { + MCAuto tmp(nodeMap->findIdsNotEqual(-1)); + nbElemsNotM1 = tmp->getNbOfElems(); + } + MCAuto nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1); + double * cooPartRef(mPartRef->_coords->getPointer()); + double * cooPartCand(mPartCand->_coords->getPointer()); + for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartCand+SPACEDIM*ii); + + + // Eliminate all edges for which y or z is not null + MCAuto baryPart = mPartCand->computeCellCenterOfMass(); + vector compo; compo.push_back(1); + MCAuto baryPartY = baryPart->keepSelectedComponents(compo); + compo[0] = 2; + MCAuto baryPartZ = baryPart->keepSelectedComponents(compo); + MCAuto idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps); + MCAuto idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps); + MCAuto idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2); + if (!idsGoodLine->getNumberOfTuples()) + continue; + + // Now the ordering along the Ox axis: + std::vector insidePoints, hitSegs; + bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode], + mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(), + idsGoodLine->begin(), idsGoodLine->end(), + /*out*/insidePoints, hitSegs); + // Optim: smaller segments completely included in eIdx and not split won't need any further treatment: + for (vector::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its) + hit[cands2[*its]] = true; + + if (!isSplit) // current segment remains in one piece + continue; + + // Get original node IDs in global coords array + for (std::vector::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit) + *iit = nodeMapInv->begin()[*iit]; + + vector polyIndices, packsIds, facePack; + // For each face implying this edge + for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++) + { + int faceIdx = revDescP2[ii]; + // each cell where this face is involved connectivity will be modified: + ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]); + + // Current face connectivity + const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1; + const int * sIdxConnE = cDesc + cIDesc[faceIdx+1]; + + vector modifiedFace; + ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace); + modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON); + connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size()); + } + } + + // Rebuild 3D connectivity from descending: + MCAuto newConn(MEDCouplingSkyLineArray::New()); + MCAuto superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0); + MCAuto idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0); + MCAuto vals(DataArrayInt::New()); vals->alloc(0); + newConn->set3(superIdx, idx, vals); + for(std::size_t ii = 0; ii < getNumberOfCells(); ii++) + for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++) + { + int sz, faceIdx = abs(descP[jj])-1; + bool orient = descP[jj]>0; + const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz); + if (orient) + newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type + else + { + vector rev(sz-1); + for (int kk=0; kkpushBackPack(ii, rev.data(), rev.data()+sz-1); + } + } + // And finally: + newConn->convertToPolyhedronConn(cAuto, cIAuto); + } // end step2 + + ret = ret->buildUniqueNotSorted(); + return ret.retn(); +} +