From: Anida Khizar Date: Mon, 9 May 2022 07:42:14 +0000 (+0200) Subject: Conformize meshes aligned on Oxyz to treat non-partition-like non-conformities X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=ab4076ca6e1947da2d9161eb1fcc548871cb68d0;p=tools%2Fmedcoupling.git Conformize meshes aligned on Oxyz to treat non-partition-like non-conformities --- diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 38989739e..541d40083 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -245,6 +245,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayIdType *colinearize2D(double eps); MEDCOUPLING_EXPORT DataArrayIdType *colinearizeKeepingConform2D(double eps); MEDCOUPLING_EXPORT DataArrayIdType *conformize3D(double eps); + MEDCOUPLING_EXPORT DataArrayIdType *conformize3DIJK(double eps); MEDCOUPLING_EXPORT mcIdType split2DCells(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI, const DataArrayIdType *midOpt=0, const DataArrayIdType *midOptI=0); MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da); MEDCOUPLING_EXPORT static MCAuto Build1DMeshFromCoords(DataArrayDouble *da); @@ -358,6 +359,8 @@ namespace MEDCoupling DataArrayIdType *internalColinearize2D(double eps, bool stayConform); template void renumberNodesInConnT(const MAPCLS& newNodeNumbersO2N); + void conformize3DEdges(const double * coo, double eps, MCAuto& c, MCAuto& cI, MCAuto& ret); + public: MEDCOUPLING_EXPORT static DataArrayIdType *ComputeRangesFromTypeDistribution(const std::vector& code); MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=25; diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx index 7bbe37f8e..a2ce44bc3 100755 --- a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx @@ -42,6 +42,7 @@ #include "InterpKernelGeo2DEdgeArcCircle.hxx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" #include "MEDCouplingUMesh_internal.hxx" +#include "TranslationRotationMatrix.hxx" #include #include @@ -1842,3 +1843,177 @@ void MEDCouplingUMesh::attractSeg3MidPtsAroundNodesUnderground(double ratio, con } } } + + +void MEDCouplingUMesh::conformize3DEdges(const double * coo, double eps, MCAuto& c, MCAuto& cI, MCAuto& ret) +{ + + static const int SPACEDIM=3; + // Now we have a face-conform mesh. + + // Recompute descending + MCAuto desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New()); + // Rebuild desc connectivity with orientation this time!! + MCAuto mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI)); + const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer()); + const mcIdType *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(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::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 mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); + const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); + MCAuto bboxArr(mDesc2->getBoundingBoxForBBTree(eps)); + const double *bbox2(bboxArr->begin()); + mcIdType 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: + std::vector > S; + for(mcIdType i=0;i < lens->getNumberOfTuples();i++) + { + std::pair p = std::make_pair(lens->getIJ(i, 0), i); + S.push_back(p); + } + sort(S.rbegin(),S.rend()); // reverse sort + + std::vector hit(nDesc2Cell); + fill(hit.begin(), hit.end(), false); + + for( std::vector >::const_iterator it = S.begin(); it != S.end(); it++) + { + mcIdType eIdx = (*it).second; + if (hit[eIdx]) + continue; + + std::vector candidates, cands2; + myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates); + // Keep only candidates colinear with current edge + double vCurr[3]; + mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; + for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? + vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3]; + for(std::vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) + { + double vOther[3]; + mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; + for (mcIdType 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 + mcIdType startNode = cDesc2[cIDesc2[eIdx]+1]; + mcIdType 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()); + mcIdType 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 (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (mcIdType 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(); + std::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 (std::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]; + + std::vector polyIndices, packsIds, facePack; + // For each face implying this edge + for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++) + { + mcIdType 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 mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1; + const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1]; + + std::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(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0); + MCAuto idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0); + MCAuto vals(DataArrayIdType::New()); vals->alloc(0); + newConn->set3(superIdx, idx, vals); + mcIdType nbCells=getNumberOfCells(); + for(mcIdType ii = 0; ii < nbCells; ii++) + for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++) + { + mcIdType sz, faceIdx = abs(descP[jj])-1; + bool orient = descP[jj]>0; + const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz); + if (orient) + newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type + else + { + std::vector rev(sz-1); + for (mcIdType kk=0; kkpushBackPack(ii, rev.data(), rev.data()+sz-1); + } + } + // And finally: + newConn->convertToPolyhedronConn(c, cI); + +} + + + + diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index e7cbaf8f3..4f5e96377 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -1129,6 +1129,7 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh // Build circular permutation to shift consecutive edges together renumb->iota(i+1); renumb->applyModulus(nbCellsInSplitMesh1D); + splitMesh1D->renumberCells(renumbP, false); cSplitPtr = splitMesh1D->getNodalConnectivity()->begin(); ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin(); @@ -1176,6 +1177,7 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh for(mcIdType mm=0;mm cIAuto; cIAuto.takeRef(_nodal_connec_index); connSla->convertToPolyhedronConn(cAuto, cIAuto); + + /************************ + * STEP 2 -- edges + ************************/ + conformize3DEdges(coo, eps, cAuto, cIAuto, ret); + + + ret = ret->buildUniqueNotSorted(); + return ret.retn(); +} + + +/*! + * \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 for meshes aligned on Oxyz. + * + * 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 DataArrayIdType * - 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. + * \throw If the mesh is not aligned on Oxyz. + * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly() + */ +DataArrayIdType *MEDCouplingUMesh::conformize3DIJK(double eps) +{ + using namespace std; + + static const int SPACEDIM=3; + checkConsistencyLight(); + if(getSpaceDimension()!=3 || getMeshDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK : 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::conformize3DIJK : This method only works for polyhedrons! Call convertAllToPoly first."); + + MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); + const double * coo(_coords->begin()); + MCAuto ret(DataArrayIdType::New()); ret->alloc(0,1); + { - /************************ - * STEP 2 -- edges - ************************/ - // Now we have a face-conform mesh. - - // Recompute descending - MCAuto desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New()); - // Rebuild desc connectivity with orientation this time!! - MCAuto mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI)); + /************************* + * STEP 1 -- faces + *************************/ + MCAuto descDNU(DataArrayIdType::New()),descIDNU(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New()); + MCAuto mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI)); const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); - const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer()); const mcIdType *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(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::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 mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); - const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); - MCAuto bboxArr(mDesc2->getBoundingBoxForBBTree(eps)); - const double *bbox2(bboxArr->begin()); - mcIdType 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(mcIdType i=0;i < lens->getNumberOfTuples();i++) + MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity())); + + // Build BBTree + MCAuto bboxArr(mDesc->getBoundingBoxForBBTree(eps)); + const double *bbox(bboxArr->begin()); getCoords()->begin(); + mcIdType 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(mcIdType i=0;i < surfs->getNumberOfTuples();i++) { - pair p = make_pair(lens->getIJ(i, 0), i); + pair p = make_pair(surfs->begin()[i], i); S.push_back(p); } sort(S.rbegin(),S.rend()); // reverse sort - - vector hit(nDesc2Cell); + 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++) { - mcIdType eIdx = (*it).second; - if (hit[eIdx]) - continue; + mcIdType faceIdx = (*it).second; + if (hit[faceIdx]) continue; vector candidates, cands2; - myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates); - // Keep only candidates colinear with current edge - double vCurr[3]; - mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; - for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? - vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3]; + 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++) { - double vOther[3]; - mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; - for (mcIdType 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) + bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps); + if (*it2 != faceIdx && col) cands2.push_back(*it2); } - if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above + if (!cands2.size()) continue; - // Now rotate edges to bring them on Ox - mcIdType startNode = cDesc2[cIDesc2[eIdx]+1]; - mcIdType endNode = cDesc2[cIDesc2[eIdx]+2]; + 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 + + // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later 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()); - mcIdType nbElemsNotM1; - { - MCAuto tmp(nodeMap->findIdsNotEqual(-1)); - nbElemsNotM1 = tmp->getNbOfElems(); - } - MCAuto nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1); + INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]), + coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]), + coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation); + double * cooPartRef(mPartRef->_coords->getPointer()); double * cooPartCand(mPartCand->_coords->getPointer()); for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) @@ -2616,82 +2641,161 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) rotation.transform_vector(cooPartCand+SPACEDIM*ii); - // Eliminate all edges for which y or z is not null + // Localize faces in 2D thanks to barycenters MCAuto baryPart = mPartCand->computeCellCenterOfMass(); - vector compo; compo.push_back(1); - MCAuto baryPartY = baryPart->keepSelectedComponents(compo); - compo[0] = 2; + vector compo; compo.push_back(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()) + MCAuto idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps); + if (!idsGoodPlane->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; + baryPart = baryPart->selectByTupleId(*idsGoodPlane); - // Get original node IDs in global coords array - for (std::vector::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit) - *iit = nodeMapInv->begin()[*iit]; + compo[0] = 0; compo.push_back(1); + MCAuto baryPartXY = baryPart->keepSelectedComponents(compo); + mPartRef->changeSpaceDimension(2,0.0); + MCAuto cc(DataArrayIdType::New()), ccI(DataArrayIdType::New()); + mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI); - vector polyIndices, packsIds, facePack; - // For each face implying this edge - for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++) - { - mcIdType 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 mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1; - const mcIdType * 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()); - } - } + if (!cc->getNumberOfTuples()) + continue; + MCAuto dsi = ccI->deltaShiftIndex(); - // Rebuild 3D connectivity from descending: - MCAuto newConn(MEDCouplingSkyLineArray::New()); - MCAuto superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0); - MCAuto idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0); - MCAuto vals(DataArrayIdType::New()); vals->alloc(0); - newConn->set3(superIdx, idx, vals); - mcIdType nbCells=getNumberOfCells(); - for(mcIdType ii = 0; ii < nbCells; ii++) - for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++) { - mcIdType sz, faceIdx = abs(descP[jj])-1; - bool orient = descP[jj]>0; - const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz); - if (orient) - newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type - else + MCAuto tmp = dsi->findIdsInRange(0, 2); + if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) { - vector rev(sz-1); - for (mcIdType kk=0; kkpushBackPack(ii, rev.data(), rev.data()+sz-1); + ostringstream oss; + oss << "MEDCouplingUMesh::conformize3DIJK: Non expected non-conformity. Only simple non-conformities are handled. Face #" << faceIdx << " violates this condition!"; + throw INTERP_KERNEL::Exception(oss.str()); } } - // And finally: - newConn->convertToPolyhedronConn(cAuto, cIAuto); - } // end step2 + + MCAuto ids = dsi->findIdsEqual(1); + // Boundary face: + if (!ids->getNumberOfTuples()) + continue; + + const mcIdType * idsGoodPlaneP(idsGoodPlane->begin()); + vector goFaces; //Faces that intersect current face + for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++) + { + mcIdType faceIdx2 = cands2[idsGoodPlaneP[*ii]]; + goFaces.push_back(faceIdx2); + hit[faceIdx2] = true; + } + + MCAuto meshGoodPlane(mDesc->buildPartOfMySelf(&goFaces[0], &goFaces[0]+goFaces.size(), false)); // false=zipCoords is called + + double faceIdxNormal = normalsP[faceIdx]; + int direction = -1; + for(int d=0; d<3; d++) + { + double diff = fabs(normalsP[faceIdx*3+d]) - 1; + if(fabs(diff) < eps) + direction = d; + } + if(direction == -1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK: Not a mesh aligned on Oxyz !"); + + // Projecting the plane in a 2D space + const DataArrayDouble *cooGoodPlane(meshGoodPlane->getCoords()); + double extraDim = cooGoodPlane->begin()[direction]; + map > projection = {{0, make_pair(1,2)}, {1, make_pair(0,2)}, {2, make_pair(0,1)}}; + int dir1 = projection.at(direction).first, dir2 = projection.at(direction).second; + vector compo2D; compo2D.push_back(dir1); compo2D.push_back(dir2); + MCAuto coo2D = cooGoodPlane->keepSelectedComponents(compo2D); + meshGoodPlane->setCoords(coo2D); + + // Intersection between the plane and the current face + const double vec[3] = {0.,0.,1.}; + mPartRef->changeSpaceDimension(3,0.0); meshGoodPlane->changeSpaceDimension(3,0.0); + mPartRef->orientCorrectly2DCells(vec, false); meshGoodPlane->orientCorrectly2DCells(vec, false); + mPartRef->changeSpaceDimension(2,0.0); meshGoodPlane->changeSpaceDimension(2,0.0); + DataArrayIdType *cellIdInPlane(0),*cellIdInCurrentFace(0); + MCAuto mInter=MEDCouplingUMesh::Intersect2DMeshes(meshGoodPlane,mPartRef,eps,cellIdInPlane,cellIdInCurrentFace); + const DataArrayDouble *cooInter(mInter->getCoords()); + const double * coo_InterD(cooInter->begin()); + + goFaces.push_back(faceIdx); + mcIdType index = 0; + for( vector::const_iterator face=goFaces.begin(); face!=goFaces.end(); ++face, ++index ) + { + // For all polyhedrons using this face, replace connectivity: + vector polyIndices, packsIds, facePack; + for (mcIdType ii=revDescIP[*face]; ii < revDescIP[*face+1]; ii++) + polyIndices.push_back(revDescP[ii]); + ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size()); + + // Current face connectivity + const mcIdType * sIdxConn = cDesc + cIDesc[*face] + 1; + const mcIdType * sIdxConnE = cDesc + cIDesc[*face+1]; + connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); + + // Deletion of old faces + mcIdType jj=0; + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj) + { + if (packsIds[jj] == -1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK: Could not find face in connectivity! Internal error."); + else + connSla->deletePack(*it2, packsIds[jj]); + } + + // Insertion of new faces coming from the intersection + MCAuto newFaces= *face==goFaces.back() ? cellIdInCurrentFace->findIdsEqual(0) : cellIdInPlane->findIdsEqual(index); + if(newFaces->getNumberOfTuples()) + { + for (const mcIdType *conformFace = newFaces->begin(); conformFace != newFaces->end(); conformFace++) + { + MCAuto connSlaDescIntersectMesh(MEDCouplingSkyLineArray::New(mInter->getNodalConnectivityIndex(), mInter->getNodalConnectivity())); + mcIdType facePack2Sz; + const mcIdType * facePack2 = connSlaDescIntersectMesh->getSimplePackSafePtr(*conformFace, facePack2Sz); // contains the type! + + vector origNodes; + for(mcIdType nodeI=1; nodeI useArray(cooIn3DMesh_,false,DeallocType::C_DEALLOC,1,3); + MCAuto nodeIndexIn3DMesh = _coords->findClosestTupleId(cooIn3DMesh); + origNodes.push_back(nodeIndexIn3DMesh->begin()[0]); + } + + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2) + connSla->pushBackPack(*it2, origNodes.data(), origNodes.data()+origNodes.size()); // without the type + } + } + else + { + mcIdType facePack2Sz; + const mcIdType * facePack2 = connSlaDesc->getSimplePackSafePtr(*face, 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 + } + } + cellIdInCurrentFace->decrRef(); + cellIdInPlane->decrRef(); + } + } // 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 + ************************/ + conformize3DEdges(coo, eps, cAuto, cIAuto, ret); ret = ret->buildUniqueNotSorted(); return ret.retn(); } + + diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 50c0ab5ba..4c41a5b1e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -2075,6 +2075,7 @@ namespace MEDCoupling //tools DataArrayIdType *conformize2D(double eps); DataArrayIdType *conformize3D(double eps); + DataArrayIdType *conformize3DIJK(double eps); DataArrayIdType *colinearize2D(double eps); DataArrayIdType *colinearizeKeepingConform2D(double eps); void shiftNodeNumbersInConn(int delta); diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 3887e6f72..d0aa693b3 100644 --- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py @@ -1505,5 +1505,40 @@ class MEDCouplingIntersectTest(unittest.TestCase): self.assertEqual(set([18]), set(ret.getValues())) pass + def testSwig2Conformize3DIJK(self): + eps = 1.0e-8 + mesh = MEDCouplingUMesh('merge', 3) + coo = DataArrayDouble([(0.0, 0.0, -0.5), (1.0, 0.0, -0.5), (0.0, 1.0, -0.5), (1.0, 1.0, -0.5), (0.0, 0.0, 0.5), (1.0, 0.0, 0.5), (0.0, 1.0, 0.5), (1.0, 1.0, 0.5), (0.0, 0.0, 1.5), (1.0, 0.0, 1.5), (0.0, 1.0, 1.5), (1.0, 1.0, 1.5), (-1.0, 0.0, 0.0), (0.0, 0.0, 0.0), (-1.0, 1.0, 0.0), (0.0, 1.0, 0.0), (-1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-1.0, 1.0, 1.0), (0.0, 1.0, 1.0)]) + mesh.setCoords(coo) + c = DataArrayInt([31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3, -1, 3, 7, 5, 1, 31, 5, 4, 6, 7, -1, 9, 11, 10, 8, -1, 5, 9, 8, 4, -1, 4, 8, 10, 6, -1, 6, 10, 11, 7, -1, 7, 11, 9, 5, 31, 13, 12, 14, 15, -1, 17, 19, 18, 16, -1, 13, 17, 16, 12, -1, 12, 16, 18, 14, -1, 14, 18, 19, 15, -1, 15, 19, 17, 13]) + cI = DataArrayInt([0, 30, 60, 90]) + mesh.setConnectivity(c, cI) + + ret = mesh.conformize3DIJK(eps) + + mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity() + mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity() + c0, cI0 = mesh.getNodalConnectivity().getValues(), mesh.getNodalConnectivityIndex().getValues() + c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues() + c2, cI2 = mretDesc2.getNodalConnectivity().getValues(), mretDesc2.getNodalConnectivityIndex().getValues() + + cRef0 = [31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 13, 0, -1, 2, 15, 6, 7, 3, -1, 3, 7, 5, 1, -1, 15, 6, 4, 13, -1, 13, 0, 2, 15, 31, 4, 6, 7, 5, -1, 9, 11, 10, 8, -1, 5, 9, 8, 17, 4, -1, 6, 19, 10, 11, 7, -1, 7, 11, 9, 5, -1, 17, 4, 6, 19, -1, 19, 10, 8, 17, 31, 13, 12, 14, 15, -1, 17, 19, 18, 16, -1, 13, 4, 17, 16, 12, -1, 12, 16, 18, 14, -1, 14, 18, 19, 6, 15, -1, 15, 6, 4, 13, -1, 17, 4, 6, 19] + cIRef0 = [0, 37, 74, 111] + cRef = [5, 1, 0, 2, 3, 5, 5, 7, 6, 4, 5, 1, 5, 4, 13, 0, 5, 2, 15, 6, 7, 3, 5, 3, 7, 5, 1, 5, 15, 6, 4, 13, 5, 13, 0, 2, 15, 5, 9, 11, 10, 8, 5, 5, 9, 8, 17, 4, 5, 6, 19, 10, 11, 7, 5, 7, 11, 9, 5, 5, 17, 4, 6, 19, 5, 19, 10, 8, 17, 5, 13, 12, 14, 15, 5, 17, 19, 18, 16, 5, 13, 4, 17, 16, 12, 5, 12, 16, 18, 14, 5, 14, 18, 19, 6, 15] + cIRef = [0, 5, 10, 16, 22, 27, 32, 37, 42, 48, 54, 59, 64, 69, 74, 79, 85, 90, 96] + cRef2 = [1, 1, 0, 1, 0, 2, 1, 2, 3, 1, 3, 1, 1, 5, 7, 1, 7, 6, 1, 6, 4, 1, 4, 5, 1, 1, 5, 1, 4, 13, 1, 13, 0, 1, 2, 15, 1, 15, 6, 1, 7, 3, 1, 13, 15, 1, 9, 11, 1, 11, 10, 1, 10, 8, 1, 8, 9, 1, 5, 9, 1, 8, 17, 1, 17, 4, 1, 6, 19, 1, 19, 10, 1, 11, 7, 1, 19, 17, 1, 13, 12, 1, 12, 14, 1, 14, 15, 1, 19, 18, 1, 18, 16, 1, 16, 17, 1, 16, 12, 1, 18, 14] + cIRef2 = [0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78, 81, 84, 87, 90, 93, 96, 99, 102] + self.assertEqual(18, mretDesc.getNumberOfCells()) + self.assertEqual(34, mretDesc2.getNumberOfCells()) + self.assertEqual(cRef0, c0) + self.assertEqual(cIRef0, cI0) + self.assertEqual(cRef, c) + self.assertEqual(cIRef, cI) + self.assertEqual(cRef2, c2) + self.assertEqual(cIRef2, cI2) + self.assertEqual(set([0,1,2]), set(ret.getValues())) + pass + + if __name__ == '__main__': unittest.main()