From 4c2430b2dfd804d83146d71b98452075ee7c5e96 Mon Sep 17 00:00:00 2001 From: abn Date: Fri, 17 Mar 2017 16:43:00 +0100 Subject: [PATCH] Working on edges ... not satisfactory yet ! --- .../TranslationRotationMatrix.hxx | 24 +- src/MEDCoupling/MEDCouplingSkyLineArray.cxx | 104 ++- src/MEDCoupling/MEDCouplingSkyLineArray.hxx | 9 +- src/MEDCoupling/MEDCouplingUMesh.hxx | 5 + .../MEDCouplingUMesh_intersection.cxx | 614 ++++++++++-------- .../MEDCouplingBasicsTest5.py | 34 +- src/MEDCoupling_Swig/MEDCouplingCommon.i | 20 +- .../MEDCouplingIntersectTest.py | 97 ++- 8 files changed, 582 insertions(+), 325 deletions(-) diff --git a/src/INTERP_KERNEL/TranslationRotationMatrix.hxx b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx index 94d7fd33f..90a0f1fa7 100644 --- a/src/INTERP_KERNEL/TranslationRotationMatrix.hxx +++ b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx @@ -124,9 +124,7 @@ namespace INTERP_KERNEL */ static void Rotate3DTriangle(const double* PP1,const double*PP2,const double*PP3, TranslationRotationMatrix& rotation_matrix) { - //initializes - rotation_matrix.translate(PP1); - + double P1w[3]; double P2w[3]; double P3w[3]; P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2]; @@ -135,10 +133,14 @@ namespace INTERP_KERNEL // translating to set P1 at the origin for (int i=0; i<3; i++) { - P2w[i]-=PP1[i]; - P3w[i]-=PP1[i]; + P1w[i] = -PP1[i]; + P2w[i] -= PP1[i]; + P3w[i] -= PP1[i]; } + //initializes matrix + rotation_matrix.translate(P1w); + // rotating to set P2 on the Oxy plane TranslationRotationMatrix A; A.rotate_x(P2w); @@ -162,15 +164,19 @@ namespace INTERP_KERNEL */ static void Rotate3DBipoint(const double* PP1,const double*PP2, TranslationRotationMatrix& rotation_matrix) { - //initializes - rotation_matrix.translate(PP1); - + double P1w[3]; double P2w[3]; P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2]; // translating to set P1 at the origin for (int i=0; i<3; i++) - P2w[i]-=PP1[i]; + { + P1w[i] = -PP1[i]; + P2w[i] -= PP1[i]; + } + + //initializes + rotation_matrix.translate(P1w); // rotating to set P2 on the Oxy plane TranslationRotationMatrix A; diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx index b2b6d5be1..ffb964e1f 100644 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx @@ -221,6 +221,38 @@ void MEDCouplingSkyLineArray::checkSuperIndex(const std::string& func) const } } +void MEDCouplingSkyLineArray::validSuperIndex(const std::string& func, int superIndex) const +{ + if(superIndex < 0 || superIndex >= _super_index->getNbOfElems()) + { + std::ostringstream oss; + oss << "MEDCouplingSkyLineArray::" << func << ": invalid super index!"; + throw INTERP_KERNEL::Exception(oss.str()); + } +} + +void MEDCouplingSkyLineArray::validIndex(const std::string& func, int idx) const +{ + if(idx < 0 || idx >= _index->getNbOfElems()) + { + std::ostringstream oss; + oss << "MEDCouplingSkyLineArray::" << func << ": invalid index!"; + throw INTERP_KERNEL::Exception(oss.str()); + } +} + +void MEDCouplingSkyLineArray::validSuperIndexAndIndex(const std::string& func, int superIndex, int index) const +{ + validSuperIndex(func, superIndex); + int idx = _super_index->begin()[superIndex] + index; + if(idx < 0 || idx >= _index->getNbOfElems()) + { + std::ostringstream oss; + oss << "MEDCouplingSkyLineArray::" << func << ": invalid index!"; + throw INTERP_KERNEL::Exception(oss.str()); + } +} + std::string MEDCouplingSkyLineArray::simpleRepr() const { std::ostringstream oss; @@ -281,7 +313,7 @@ std::string MEDCouplingSkyLineArray::simpleRepr() const /** * For a 2- or 3-level SkyLine array, return a copy of the absolute pack with given identifier. */ -void MEDCouplingSkyLineArray::getPackSafe(const int absolutePackId, std::vector & pack) const +void MEDCouplingSkyLineArray::getSimplePackSafe(const int absolutePackId, std::vector & pack) const { if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems()) throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!"); @@ -294,7 +326,7 @@ void MEDCouplingSkyLineArray::getPackSafe(const int absolutePackId, std::vector< /** * Same as getPackSafe, but directly returns a pointer to the internal data with the size of the pack. */ -const int * MEDCouplingSkyLineArray::getPackSafePtr(const int absolutePackId, int & packSize) const +const int * MEDCouplingSkyLineArray::getSimplePackSafePtr(const int absolutePackId, int & packSize) const { if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems()) throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!"); @@ -355,6 +387,7 @@ void MEDCouplingSkyLineArray::deletePack(const int superIdx, const int idx) using namespace std; checkSuperIndex("deletePack"); + validSuperIndexAndIndex("deletePack", superIdx, idx); int * vP = _values->getPointer(); int * siP(_super_index->getPointer()), *iP(_index->getPointer()); @@ -383,6 +416,7 @@ void MEDCouplingSkyLineArray::pushBackPack(const int superIdx, const int * packB using namespace std; checkSuperIndex("pushBackPack"); + validSuperIndex("pushBackPack", superIdx); int *siP(_super_index->getPointer()), *iP(_index->getPointer()); const int sz(distance(packBg, packEnd)); @@ -408,33 +442,73 @@ void MEDCouplingSkyLineArray::pushBackPack(const int superIdx, const int * packB (siP[ii])++; } +/** + * Replace pack with absolute index 'idx' with the provided new pack. Function can be used either + * for 2-level SkyLine or 3-level SkyLine. + */ +void MEDCouplingSkyLineArray::replaceSimplePack(const int idx, const int * packBg, const int * packEnd) +{ + using namespace std; + + validIndex("replaceSimplePack", idx); + + int * iP(_index->getPointer()); + int newSz = distance(packBg, packEnd); + const int start = iP[idx], end = iP[idx+1]; + + // _values + int initValSz = _values->getNbOfElems(); + int deltaSz = newSz-(end-start); // can be negative + if (deltaSz) + { + if (deltaSz > 0) + _values->reAlloc(initValSz+deltaSz); + int *vP(_values->getPointer()); + copy(vP+end, vP+initValSz, vP+end+deltaSz); + if (deltaSz < 0) + _values->reAlloc(initValSz+deltaSz); + } + + // copy new pack + copy(packBg, packEnd, _values->getPointer()+start); + + // _index + for(int ii = idx+1; ii < _index->getNbOfElems(); ii++) + iP[ii] += deltaSz; +} /** * Replace pack with super index 'superIdx' and index 'idx' with the provided new pack. + * Function can be used only for 3-level SkyLine. */ void MEDCouplingSkyLineArray::replacePack(const int superIdx, const int idx, const int * packBg, const int * packEnd) { using namespace std; checkSuperIndex("replacePack"); + validSuperIndexAndIndex("replacePack", superIdx, idx); - int * vP = _values->getPointer(); int * siP(_super_index->getPointer()), *iP(_index->getPointer()); int newSz = distance(packBg, packEnd); const int start = iP[siP[superIdx]+idx], end = iP[siP[superIdx]+idx+1]; - // _values - copy(vP+end, vP+_values->getNbOfElems(), vP+start); - _values->reAlloc(_values->getNbOfElems() - (end-start)); - // _index - int nt = _index->getNbOfElems(); - copy(iP+siP[superIdx]+idx+1, iP+nt, iP+siP[superIdx]+idx); - _index->reAlloc(nt-1); iP = _index->getPointer(); // better not forget this ... - for(int ii = siP[superIdx]+idx; ii < nt-1; ii++) - iP[ii] -= (end-start); + // _values + int initValSz = _values->getNbOfElems(); + int deltaSz = newSz-(end-start); // can be negative + if (deltaSz) + { + if (deltaSz > 0) + _values->reAlloc(initValSz+deltaSz); + int *vP(_values->getPointer()); + copy(vP+end, vP+initValSz, vP+end+deltaSz); + if (deltaSz < 0) + _values->reAlloc(initValSz+deltaSz); + } - // _super_index - for(int ii = superIdx+1; ii < _super_index->getNbOfElems(); ii++) - (siP[ii])--; + // copy new pack + copy(packBg, packEnd, _values->getPointer()+start); + // _index + for(int ii = siP[superIdx]+idx+1; ii < _index->getNbOfElems(); ii++) + iP[ii] += deltaSz; } diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx index d6d601ad1..bf5d74cab 100644 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx @@ -84,14 +84,15 @@ namespace MEDCoupling std::string simpleRepr() const; - void getPackSafe(const int absolutePackId, std::vector & pack) const; - const int * getPackSafePtr(const int absolutePackId, int & packSize) const; + void getSimplePackSafe(const int absolutePackId, std::vector & pack) const; + const int * getSimplePackSafePtr(const int absolutePackId, int & packSize) const; void findPackIds(const std::vector & superPackIndices, const int *packBg, const int *packEnd, std::vector& out) const; void deletePack(const int superIdx, const int idx); void pushBackPack(const int superIdx, const int * packBg, const int * packEnd); + void replaceSimplePack(const int idx, const int * packBg, const int * packEnd); void replacePack(const int superIdx, const int idx, const int * packBg, const int * packEnd); void convertToPolyhedronConn( MCAuto& c, MCAuto& cI) const; @@ -101,10 +102,14 @@ namespace MEDCoupling ~MEDCouplingSkyLineArray(); void checkSuperIndex(const std::string& func) const; + void validSuperIndex(const std::string& func, int superIndex) const; + void validIndex(const std::string& func, int index) const; + void validSuperIndexAndIndex(const std::string& func, int superIndex, int index) const; MCAuto _super_index; MCAuto _index; MCAuto _values; }; + } # endif diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 6192f0565..ccdae4a34 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -347,6 +347,11 @@ namespace MEDCoupling int split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI); static bool Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords); static void ComputeAllTypesInternal(std::set& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex); + static bool OrderPointsAlongLine(const double * coo, int startNode, int endNode, + const int * c, const int * cI, const int *idsBg, const int *endBg, + std::vector & pointIds); + static void ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode, + const std::vector& insidePoints, std::vector& modifiedFace); public: MEDCOUPLING_EXPORT static DataArrayInt *ComputeRangesFromTypeDistribution(const std::vector& code); MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=24; diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index e2ad0792a..940c9199b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -44,6 +44,7 @@ #include #include #include +#include using namespace MEDCoupling; @@ -2000,14 +2001,75 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) ///@cond INTERNAL /** - * + * c, cI describe a wire mesh in 3D space. + *\return true if the segment is indeed split */ -bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, const double * startNode, const double *endNode, +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) { + 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; + continue; + } + if (idx == endNode || idx == startNode) + break; + pointIds.push_back(idx); + } + if (go == 2) + reverse(pointIds.begin(), pointIds.end()); + return (pointIds.size() != 0); +} + +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()); + 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 > 0 || d == (dst-1)) + modifiedFace.insert(++startPos, insidePoints.begin(), insidePoints.end()); + else + modifiedFace.insert(++endPos, insidePoints.rbegin(), insidePoints.rend()); } + ///@endcond @@ -2041,290 +2103,320 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) 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 desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); - MCAuto mDesc(buildDescendingConnectivity(desc,descI,revDesc,revDescI)); - const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + int *c(getNodalConnectivity()->getPointer()),*cI(getNodalConnectivityIndex()->getPointer()); MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); - const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); - MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity())); const double * coo(_coords->begin()); + MCAuto ret(DataArrayInt::New()); + std::cout << "BEGIN\n"; + std::cout << connSla->simpleRepr(); - /************************* - * STEP 1 -- faces - *************************/ - // Build BBTree - MCAuto bboxArr(mDesc->getBoundingBoxForBBTree()); - const double *bbox(bboxArr->begin()),*coords(getCoords()->begin()); - int nCell=getNumberOfCells(), 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> S; - for(int i=0;i < surfs->getNumberOfTuples();i++){ - pair p = make_pair(surfs->getIJ(i, 0), 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. + { + /************************* + * 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()); + const double *bbox(bboxArr->begin()),*coords(getCoords()->begin()); + int nCell=getNumberOfCells(), 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> S; + for(int 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 sIdx = (*it).second; - if (hit[sIdx]) continue; + for( vector>::const_iterator it = S.begin(); it != S.end(); it++) + { + int sIdx = (*it).second; + if (hit[sIdx]) continue; + + vector candidates, cands2; + myTree.getIntersectingElems(bbox+sIdx*2*SPACEDIM,candidates); + // Keep only candidates whose normal matches the normal of current face + for(vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + { + bool col = INTERP_KERNEL::isColinear3D(normalsP + sIdx*SPACEDIM, normalsP + *(it)*SPACEDIM, eps); + if (*it != sIdx && col) + cands2.push_back(*it); + } + 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[sIdx]+1]), + coo+SPACEDIM*(cDesc[cIDesc[sIdx]+2]), + coo+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation); + + MCAuto mPartRef(mDesc->buildPartOfMySelfSlice(sIdx, sIdx+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 (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartCand+SPACEDIM*ii); + double ze_z = cooPartRef[2]; + + // 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(ze_z-eps, ze_z+eps); + if (!idsGoodPlane->getNumberOfTuples()) + continue; + + 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(); - vector candidates, cands2; - myTree.getIntersectingElems(bbox+sIdx*2*SPACEDIM,candidates); - // Keep only candidates whose normal matches the normal of current face - for(vector::const_iterator it=candidates.begin();it!=candidates.end();it++) { - bool col = INTERP_KERNEL::isColinear3D(normalsP + sIdx*SPACEDIM, normalsP + *(it)*SPACEDIM, eps); - if (*it != sIdx && col) - cands2.push_back(*it); + MCAuto tmp = dsi->findIdsInRange(0, 2); + if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled."); } - 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[sIdx]+1]), - coo+SPACEDIM*(cDesc[cIDesc[sIdx]+2]), - coo+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation); - - MCAuto mPartRef(mDesc->buildPartOfMySelfSlice(sIdx, sIdx+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 (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) - rotation.transform_vector(cooPartRef+SPACEDIM*ii); - for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) - rotation.transform_vector(cooPartCand+SPACEDIM*ii); - double ze_z = cooPartRef[2]; - - // 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(ze_z-eps, ze_z+eps); - if (!idsGoodPlane->getNumberOfTuples()) - continue; - - 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()) - throw INTERP_KERNEL::Exception("Non expected non-conformity. Only simple (=partition-like) non-conformities are handled."); + MCAuto ids = dsi->findIdsEqual(1); + // Boundary face: + if (!ids->getNumberOfTuples()) + continue; + + double checkSurf=0.0; + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + { + int faceIdx = cands2[idsGoodPlane->begin()[*ii]]; + hit[faceIdx] = true; + checkSurf += surfs->begin()[faceIdx]; + } + if (fabs(checkSurf - surfs->begin()[sIdx]) > eps) + { + ostringstream oss; + oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << sIdx << " 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[sIdx]; ii < revDescIP[sIdx+1]; ii++) + polyIndices.push_back(revDescP[ii]); + ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size()); + + // Current face connectivity + const int * sIdxConn = cDesc + cIDesc[sIdx] + 1; + const int * sIdxConnE = cDesc + cIDesc[sIdx+1]; + connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); + // Deletion of old faces + int jj=0; + for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it, ++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(*it, packsIds[jj]); + } + // Insertion of new faces: + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + { + // Build pack from the face to insert: + int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)]; + int facePack2Sz; +// vector toto; +// connSlaDesc->getSimplePackSafe(faceIdx, toto); + const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx, facePack2Sz); // contains the type! + // Insert it in all hit polyhedrons: + for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it) + connSla->pushBackPack(*it, facePack2+1, facePack2+facePack2Sz); // without the type + } } - // TODO check matching total surface! - - MCAuto ids = dsi->findIdsEqual(1); - // Boundary face: - if (!ids->getNumberOfTuples()) - continue; - - for (const int * ii = ids->begin(); ii != ids->end(); ii++) - hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true; - - // For all polyhedrons using this face, replace connectivity: - vector polyIndices, packsIds, facePack; - for (int ii=*(revDescIP+sIdx); ii < *(revDescIP+sIdx+1); ii++) - polyIndices.push_back(*(revDescP+ii)); - // Current face connectivity - const int * sIdxConn = cDesc + cIDesc[sIdx] + 1; - const int * sIdxConnE = cDesc + cIDesc[sIdx+1]; - connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); - // Deletion of old faces - int jj=0; - for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it, ++jj) - connSla->deletePack(*it, packsIds[jj]); - // Insertion of new faces: - for (const int * ii = ids->begin(); ii != ids->end(); ii++) - { - // Build pack from the face to insert: - int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)]; - int facePack2Sz; - vector toto; - connSlaDesc->getPackSafe(faceIdx, toto); - const int * facePack2 = connSlaDesc->getPackSafePtr(faceIdx, facePack2Sz); // contains the type! - // Insert it in all hit polyhedrons: - for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it) - connSla->pushBackPack(*it, facePack2+1, facePack2+facePack2Sz); // without the type - } - } + } // 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. - mDesc = buildDescendingConnectivity(desc,descI,revDesc,revDescI); - MCAuto desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New()); - MCAuto mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2); - const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); - bboxArr = mDesc2->getBoundingBoxForBBTree(); - const double *bbox2(bboxArr->begin()); - int nDesc2Cell=mDesc2->getNumberOfCells(); - BBTree myTree(bbox2,0,0,nDesc2Cell,-eps); - - // Edges - handle longest first - MCAuto lenF = mDesc2->getMeasureField(true); - DataArrayDouble * lens = lenF->getArray(); - - // Sort edges by decreasing length: - vector> S2; - for(int i=0;i < lens->getNumberOfTuples();i++){ - pair p = make_pair(lens->getIJ(i, 0), i); - S2.push_back(p); - } - sort(S.rbegin(),S.rend()); // reverse sort - hit.resize(nDesc2Cell); - fill(hit.begin(), hit.end(), false); + std::cout << "BEFORE STEP2\n"; + std::cout << connSla->simpleRepr(); - for( vector>::const_iterator it = S.begin(); it != S.end(); it++) - { - int eIdx = (*it).second; - if (hit[eIdx]) continue; - - vector candidates, cands2; - myTree.getIntersectingElems(bbox+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 it=candidates.begin();it!=candidates.end();it++) - { - double vOther[3]; - unsigned start2 = cDesc2[cIDesc2[*it]+1], end2 = cDesc2[cIDesc2[*it]+2]; - for (int i3=0; i3 < 3; i3++) - vCurr[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3]; - bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps); - if (*it != eIdx && col) - cands2.push_back(*it); - } - if (!cands2.size()) - continue; - - // Now rotate edges to bring them on Ox - const double * startNode = coo+SPACEDIM*(cDesc[cIDesc[eIdx]+1]); - const double * endNode = coo+SPACEDIM*(cDesc[cIDesc[eIdx]+2]); - INTERP_KERNEL::TranslationRotationMatrix rotation; - INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(startNode, 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()); - double * cooPartRef(mPartRef->_coords->getPointer()); - double * cooPartCand(mPartCand->_coords->getPointer()); - for (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) - rotation.transform_vector(cooPartRef+SPACEDIM*ii); - for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) - rotation.transform_vector(cooPartCand+SPACEDIM*ii); - double ze_y = cooPartRef[1], ze_z = cooPartRef[2]; - - // 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; - bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), startNode, endNode, - mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(), - idsGoodLine->begin(), idsGoodLine->end(), - /*out*/insidePoints); - if (!isSplit) // current segment remains in one piece - continue; - - 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(); + { + /************************ + * STEP 2 -- edges + /************************/ + // Now we have a face-conform mesh. + + // Recompute descending + MCAuto descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + // Rebuild desc connectivity with orientation this time!! + MCAuto mDesc(buildDescendingConnectivity2(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())); + std::cout << "writing VTK!\n"; + mDesc->writeVTK("/tmp/toto_desc_confInter.vtu", true); + MCAuto desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New()); + MCAuto mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2); + std::cout << "writing VTK2!\n"; + mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu", true); + const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); + const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); + MCAuto bboxArr(mDesc2->getBoundingBoxForBBTree()); + 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(int 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 + for( vector>::const_iterator it = S.begin(); it != S.end(); it++) { - MCAuto tmp = dsi->findIdsInRange(0, 2); - if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) - throw INTERP_KERNEL::Exception("Non expected non-conformity. Only simple (=partition-like) non-conformities are handled."); - } - // TODO check matching total surface! - - MCAuto ids = dsi->findIdsEqual(1); - // Boundary face: - if (!ids->getNumberOfTuples()) - continue; - - for (const int * ii = ids->begin(); ii != ids->end(); ii++) - hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true; - - // For all polyhedrons using this face, replace connectivity: - vector polyIndices, packsIds, facePack; - for (int ii=*(revDescIP+eIdx); ii < *(revDescIP+eIdx+1); ii++) - polyIndices.push_back(*(revDescP+ii)); - // Current face connectivity - const int * sIdxConn = cDesc + cIDesc[eIdx] + 1; - const int * sIdxConnE = cDesc + cIDesc[eIdx+1]; - connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); - // Deletion of old faces - int jj=0; - for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it, ++jj) - connSla->deletePack(*it, packsIds[jj]); - // Insertion of new faces: - for (const int * ii = ids->begin(); ii != ids->end(); ii++) + int eIdx = (*it).second; + + 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 it=candidates.begin();it!=candidates.end();it++) + { + double vOther[3]; + unsigned start2 = cDesc2[cIDesc2[*it]+1], end2 = cDesc2[cIDesc2[*it]+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(*it); + } + 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; { - // Build pack from the face to insert: - int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)]; - int facePack2Sz; - vector toto; - connSlaDesc->getPackSafe(faceIdx, toto); - const int * facePack2 = connSlaDesc->getPackSafePtr(faceIdx, facePack2Sz); // contains the type! - // Insert it in all hit polyhedrons: - for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it) - connSla->pushBackPack(*it, facePack2+1, facePack2+facePack2Sz); // without the type + MCAuto tmp(nodeMap->findIdsNotEqual(-1)); + nbElemsNotM1 = tmp->getNbOfElems(); } - } - + MCAuto nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1); + double * cooPartRef(mPartRef->_coords->getPointer()); + double * cooPartCand(mPartCand->_coords->getPointer()); + for (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartCand+SPACEDIM*ii); + double ze_y = cooPartRef[1], ze_z = cooPartRef[2]; + + // 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; + bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode], + mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(), + idsGoodLine->begin(), idsGoodLine->end(), + /*out*/insidePoints); + std::cout << "edge " << eIdx << "\n"; + for (int kk =0; kk < insidePoints.size(); kk++) + std::cout << " " << nodeMapInv->begin()[insidePoints[kk]] << "\n"; + + 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 = abs(revDescP2[ii]) - 1; +// +// //connSlaDesc +// +// // Current face connectivity - find it in the original 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); +// std::cout << "modifiedFace WAS " << eIdx << "\n"; +// for (const int * kk =sIdxConn; kk < sIdxConnE; kk++) +// std::cout << " " << *kk << "\n"; +// std::cout << "modifiedFace IS " << eIdx << "\n"; +// for (int kk =0; kk < modifiedFace.size(); kk++) +// std::cout << " " << modifiedFace[kk] << "\n"; +// +// int k3 = 0; +// std::cout << "BEFORE\n"; +// std::cout << connSla->simpleRepr(); +// connSlaDesc->replacePack(faceIdx, modifiedFace); +// std::cout << "AFTER\n"; +// std::cout << connSla->simpleRepr(); +// } + } + } // Set back modified connectivity - MCAuto cAuto; cAuto.takeRef(_nodal_connec); - MCAuto cIAuto; cIAuto.takeRef(_nodal_connec_index); connSla->convertToPolyhedronConn(cAuto, cIAuto); - - - MCAuto ret(DataArrayInt::New()); + ret = ret->buildUniqueNotSorted(); return ret.retn(); } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index 7f4615e91..0dd4550a3 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -1502,7 +1502,7 @@ class MEDCouplingBasicsTest5(unittest.TestCase): pass def testSwig2Hexa8HavingFacesWarped1(self): - """ This test is bases on a "error" of interpolation detected. After investigation cell #3 of src is warped that leads to the fact that when trg is + """ This test is bases on a "error" of interpolation detected. After investigation cell #3 of src is warped that leads to the fact that when trg is intersected with src the sum of intersection volume is greater than the volume of the trg cell. A test that can be done is to split the cell #3 of src into tetrohedrons and by summing all the volumes it does not fit the volume computed of cell#3 unsplitted (expect for GENERAL_24). @@ -3638,7 +3638,7 @@ class MEDCouplingBasicsTest5(unittest.TestCase): sla0.set3( superi.deepCopy(), index.deepCopy(), value.deepCopy() ) self.assertTrue( superi.isEqual( sla0.getSuperIndexArray() )) - pack = sla0.getPackSafe(2) + pack = sla0.getSimplePackSafe(2) self.assertEqual([9,10,1,12], pack) ids = sla0.findPackIds([0,1], [9,10,1,12]) self.assertEqual([-1,1], ids) @@ -3691,8 +3691,36 @@ class MEDCouplingBasicsTest5(unittest.TestCase): self.assertEqual(siRef, si.getValues()) self.assertEqual(idxRef, idx.getValues()) self.assertEqual(cRef, val.getValues()) - pass + idxRef2 = [0,4,8,12,16,20,23,26,29,32,36,40,42,46,50,53,56, 60, 64, 68, 72, 76, 80 ] + cRef2 = [1,0,2,3, 5,7,6,4, 1,5,4,0, 0,4,6,2, 2,6,7,3, 3,7,8, 7,5,8, 5,1,8, 1,3,8, + 9,1,3,10, 11,12,7,5, 300,300, 3,7,12,10, 10,12,11,9, 3,7,8, 7,5,8, + 11,5,7,12, 14,16,15,13, 11,14,13,5, 5,13,15,7, 7,15,16,12, 12,16,14,11] + sla0.replacePack(1,2, [300,300]) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual(siRef, si.getValues()) + self.assertEqual(idxRef2, idx.getValues()) + self.assertEqual(cRef2, val.getValues()) + + sla0.replacePack(1,2, [9,11,5,1]) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual(siRef, si.getValues()) + self.assertEqual(idxRef, idx.getValues()) + self.assertEqual(cRef, val.getValues()) + + sla0.replaceSimplePack(11, [300,300]) # 11 is the abs index of pack (superIdx=1,idx=2) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual(siRef, si.getValues()) + self.assertEqual(idxRef2, idx.getValues()) + self.assertEqual(cRef2, val.getValues()) + + sla0.replaceSimplePack(11, [9,11,5,1]) # 11 is the abs index of pack (superIdx=1,idx=2) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual(siRef, si.getValues()) + self.assertEqual(idxRef, idx.getValues()) + self.assertEqual(cRef, val.getValues()) + + pass def testMEDCouplingUMeshgenerateGraph(self): # cartesian mesh 3x3 diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 5e50a2422..a1a557ecc 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -1193,7 +1193,7 @@ namespace MEDCoupling DataArrayInt* getIndexArray() const; DataArrayInt* getValuesArray() const; - void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception); + void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception); %extend { @@ -1222,10 +1222,10 @@ namespace MEDCoupling return self->simpleRepr(); } - PyObject *getPackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception) + PyObject *getSimplePackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception) { std::vector ret; - self->getPackSafe(absolutePackId,ret); + self->getSimplePackSafe(absolutePackId,ret); return convertIntArrToPyList2(ret); } @@ -1247,6 +1247,20 @@ namespace MEDCoupling self->pushBackPack(i,vpack.data(), vpack.data()+vpack.size()); } + void replaceSimplePack(const int idx, PyObject *pack) throw(INTERP_KERNEL::Exception) + { + std::vector vpack; + convertPyToNewIntArr3(pack,vpack); + self->replaceSimplePack(idx, vpack.data(), vpack.data()+vpack.size()); + } + + void replacePack(const int superIdx, const int idx, PyObject *pack) throw(INTERP_KERNEL::Exception) + { + std::vector vpack; + convertPyToNewIntArr3(pack,vpack); + self->replacePack(superIdx, idx, vpack.data(), vpack.data()+vpack.size()); + } + PyObject *convertToPolyhedronConn() const throw(INTERP_KERNEL::Exception) { MCAuto d0=DataArrayInt::New(); diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 4f7ef5b81..0e2469a68 100644 --- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py @@ -316,7 +316,7 @@ class MEDCouplingIntersectTest(unittest.TestCase): def testSwig2Intersect2DMeshesQuadra1(self): import cmath def createDiagCircle(lX, lY, R, cells=[0,1]): - """ A circle in a square box, cut along the diagonal. + """ A circle in a square box, cut along the diagonal. """ c = [] for i in range(8): @@ -917,41 +917,74 @@ class MEDCouplingIntersectTest(unittest.TestCase): def testSwig2Conformize3D1(self): """ Simple test where no edge merge is required, only face merging (first part of the algo) """ - m = MEDCouplingCMesh("cube") - dac = DataArrayDouble([0.0, 1.0]) - m.setCoordsAt(0, dac); m.setCoordsAt(1, dac); m.setCoordsAt(2, dac) - m = m.buildUnstructured() - m.convertToPolyTypes(range(m.getNumberOfCells())) - coo = m.getCoords().getValues() - coo.extend([1.0, 0.5, 0.5]) - m.setCoords(DataArrayDouble(coo, len(coo)/3, 3)) - #c, cI = m.getNodalConnectivity().getValues(), m.getNodalConnectivityIndex().getValues() - ## Initial: - #c = [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] - c = [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, 8, - -1, 7, 5, 8, - -1, 5, 1, 8, - -1, 1, 3, 8] - cI = [0, len(c)] - m.setConnectivity(DataArrayInt(c), DataArrayInt(cI)) - m2 = MEDCouplingCMesh("cubes") - m2.setCoordsAt(0, DataArrayDouble([1.0, 2.0])) - m2.setCoordsAt(1, dac) - m2.setCoordsAt(2, DataArrayDouble([0.0, 1.0, 2.0])) - m2 = m2.buildUnstructured(); m2.convertToPolyTypes(range(m2.getNumberOfCells())) - mret = MEDCouplingUMesh.MergeUMeshes([m, m2]) - mret.mergeNodes(1.0e-8) - mret.writeVTK("/tmp/toto.vtu") - mretDesc, _, _, _, _ = mret.buildDescendingConnectivity() + mesh = MEDCouplingUMesh('merge', 3) + coo = DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0,0,1),(1,0,1),(0,1,1),(1,1,1),(1,0.5,0.5), + (2,0,0),(2,1,0),(2,0,1),(2,1,1),(1,0,2),(2,0,2),(1,1,2),(2,1,2)]) + 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, 8, -1, 7, 5, 8, -1, 5, 1, 8, -1, 1, 3, 8, 31, 9, 1, 3, 10, -1, 11, 12, 7, 5, -1, + 9, 11, 5, 1, -1, 1, 5, 7, 3, -1, 3, 7, 12, 10, -1, 10, 12, 11, 9, 31, 11, 5, 7, 12, -1, + 14, 16, 15, 13, -1, 11, 14, 13, 5, -1, 5, 13, 15, 7, -1, 7, 15, 16, 12, -1, 12, 16, 14, 11]) + cI = DataArrayInt([0, 41, 71, 101]) + mesh.setConnectivity(c, cI) + ret = mesh.conformize3D(1.0e-8) + + mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity() + mretDesc.writeVTK("/tmp/toto_conf_desc.vtu") + c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues() + cRef = [5, 1, 0, 2, 3, 5, 5, 7, 6, 4, 5, 1, 5, 4, 0, 5, 0, 4, 6, 2, 5, 2, 6, 7, 3, 5, 3, 7, 8, 5, + 7, 5, 8, 5, 5, 1, 8, 5, 1, 3, 8, 5, 9, 1, 3, 10, 5, 11, 12, 7, 5, 5, 9, 11, 5, 1, 5, 3, + 7, 12, 10, 5, 10, 12, 11, 9, 5, 14, 16, 15, 13, 5, 11, 14, 13, 5, 5, 5, 13, 15, 7, 5, 7, + 15, 16, 12, 5, 12, 16, 14, 11] + cIRef = [0, 5, 10, 15, 20, 25, 29, 33, 37, 41, 46, 51, 56, 61, 66, 71, 76, 81, 86, 91] + self.assertEqual(19, mretDesc.getNumberOfCells()) + self.assertEqual(cRef, c) + self.assertEqual(cIRef, cI) + self.assertEqual([1], ret.getValues()) + pass + + def testSwig2Conformize3D2(self): + """ More advanced test where edge merge is required. """ + mesh = MEDCouplingUMesh('merge', 3) + coo = DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0,0,1),(1,0,1),(0,1,1),(1,1,1),(0,0,2),(1,0,2),(0,1,2), + (1,1,2),(1,0.6,0),(1,0.6,1),(1,0.3,1),(1,0.3,2),(2,0,0),(2,1,0),(2,0,1),(2,1,1),(2,0,2),(2,1,2)]) + mesh.setCoords(coo) + c = DataArrayInt([31, 1, 0, 2, 3, 12, -1, 5, 13, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3, -1, 3, 7, 13, 12, -1, + 5, 1, 12, 13, 31, 5, 4, 6, 7, 14, -1, 9, 15, 11, 10, 8, -1, 5, 9, 8, 4, -1, 4, 8, 10, 6, -1, 6, 10, 11, 7, -1, + 7, 11, 15, 14, -1, 9, 5, 14, 15, 31, 16, 1, 3, 17, -1, 18, 19, 7, 5, -1, 16, 18, 5, 1, -1, 1, 5, 7, 3, -1, + 3, 7, 19, 17, -1, 17, 19, 18, 16, 31, 18, 5, 7, 19, -1, 20, 21, 11, 9, -1, 18, 20, 9, 5, -1, 5, 9, 11, 7, -1, + 7, 11, 21, 19, -1, 19, 21, 20, 18]) + cI = DataArrayInt([0, 37, 74, 104, 134]) + mesh.setConnectivity(c, cI) + + mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity() + mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity() mretDesc.writeVTK("/tmp/toto_desc.vtu") + mretDesc2.writeVTK("/tmp/toto_desc2.vtu") - mret.conformize3D(1.0e-8) + ret = mesh.conformize3D(1.0e-8) - mret.writeVTK("/tmp/toto_conf.vtu") - mretDesc, _, _, _, _ = mret.buildDescendingConnectivity() + mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity() + mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity() mretDesc.writeVTK("/tmp/toto_conf_desc.vtu") - return mret + mretDesc2.writeVTK("/tmp/toto_conf_desc2.vtu") + c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues() + c2, cI2 = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues() + cRef = [] + cIRef = [] + cRef2 = [] + cIRef2 = [] + print mretDesc2.getNumberOfCells() + return + self.assertEqual(22, mretDesc.getNumberOfCells()) + self.assertEqual(40, mretDesc2.getNumberOfCells()) + self.assertEqual(cRef, c) + self.assertEqual(cIRef, cI) + self.assertEqual(cRef2, c2) + self.assertEqual(cIRef2, cI2) + self.assertEqual([0,1,2,3], ret.getValues()) + pass + if __name__ == '__main__': unittest.main() -- 2.39.2