From 0ba3453939dda0697b09ed7e728a01d4f33e3ce2 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Fri, 11 Aug 2023 09:55:19 +0200 Subject: [PATCH] [EDF27988] : Implementation of MEDCouplingUMesh.explodeMeshTo for MEDFileUMesh.reduceToCells --- src/MEDCoupling/MEDCouplingUMesh.cxx | 38 +++++++++++ src/MEDCoupling/MEDCouplingUMesh.hxx | 2 + src/MEDCoupling/MEDCouplingUMesh_internal.cxx | 68 +++++++++++++++++++ .../MEDCouplingBasicsTest7.py | 65 +++++++++++++++++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 13 ++++ src/MEDLoader/MEDFileMesh.cxx | 2 +- src/MEDLoader/MEDFileMesh.hxx | 1 + src/MEDLoader/Swig/CMakeLists.txt | 2 +- src/MEDLoader/Swig/MEDLoaderFinalize.i | 3 + src/MEDLoader/Swig/MEDLoaderFinalize.py | 66 ++++++++++++++++++ src/MEDLoader/Swig/MEDLoaderTest3.py | 59 ++++++++++++++++ 11 files changed, 314 insertions(+), 5 deletions(-) create mode 100644 src/MEDLoader/Swig/MEDLoaderFinalize.py diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 0bd908c17..c06cb2838 100755 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -697,6 +697,44 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity(DataArrayIdType return buildDescendingConnectivityGen(desc,descIndx,revDesc,revDescIndx,MEDCouplingFastNbrer); } +/*! + * This method is a generalization of buildDescendingConnectivity method. This method return mesh containing subcells of this at level specified by \a targetDeltaLevel + * and return descending and reverse descending correspondances to this. + * + * \param [in] targetDeltaLevel target delta level compared to \a this mesh dimension. This parameter is expected to be lower than zero. + * + * \throw If targetDeltaLevel is greater or equal to zero + * \throw If targetDeltaLevel is lower than -meshDim + * \sa MEDCouplingUMesh::buildDescendingConnectivity, MEDCouplingUMesh::explode3DMeshTo1D + */ +MCAuto MEDCouplingUMesh::explodeMeshTo(int targetDeltaLevel, MCAuto& desc, MCAuto& descIndx, MCAuto& revDesc, MCAuto& revDescIndx) const +{ + this->checkConsistencyLight(); + if( targetDeltaLevel >= 0 ) + THROW_IK_EXCEPTION("Input parameter targetDeltaLevel is expected to be lower than zero !"); + if( targetDeltaLevel == -1 ) + { + desc = DataArrayIdType::New(); descIndx = DataArrayIdType::New(); revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New(); + MCAuto ret( this->buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx) ); + return ret; + } + if( targetDeltaLevel == -2 && this->getMeshDimension() == 3 ) + { + desc = DataArrayIdType::New(); descIndx = DataArrayIdType::New(); revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New(); + MCAuto ret( this->explode3DMeshTo1D(desc,descIndx,revDesc,revDescIndx) ); + return ret; + } + if( targetDeltaLevel == -this->getMeshDimension() ) + { + MCAuto ret = MEDCouplingUMesh::Build0DMeshFromCoords( const_cast( this->getCoords() ) ); + MEDCouplingUMesh::DeleteCellTypeInIndexedArray(getNodalConnectivity(),getNodalConnectivityIndex(),desc,descIndx); + revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New(); + this->getReverseNodalConnectivity(revDesc,revDescIndx); + return ret; + } + THROW_IK_EXCEPTION("Not valid input targetDeltaLevel regarding mesh dimension"); +} + /*! * \a this has to have a mesh dimension equal to 3. If it is not the case an INTERP_KERNEL::Exception will be thrown. * This behaves exactly as MEDCouplingUMesh::buildDescendingConnectivity does except that this method compute directly the transition from mesh dimension 3 to sub edges (dimension 1) diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 22524c028..af74fd7be 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -123,6 +123,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT bool areCellsIncludedInPolicy7(const MEDCouplingUMesh *other, DataArrayIdType *& arr) const; MEDCOUPLING_EXPORT void getReverseNodalConnectivity(DataArrayIdType *revNodal, DataArrayIdType *revNodalIndx) const; MEDCOUPLING_EXPORT MCAuto explodeIntoEdges(MCAuto& desc, MCAuto& descIndex, MCAuto& revDesc, MCAuto& revDescIndx) const; + MEDCOUPLING_EXPORT MCAuto explodeMeshTo(int targetDeltaLevel, MCAuto& desc, MCAuto& descIndx, MCAuto& revDesc, MCAuto& revDescIndx) const; MEDCOUPLING_EXPORT MEDCouplingUMesh *explode3DMeshTo1D(DataArrayIdType *desc, DataArrayIdType *descIndx, DataArrayIdType *revDesc, DataArrayIdType *revDescIndx) const; MEDCOUPLING_EXPORT MEDCouplingUMesh *buildDescendingConnectivity(DataArrayIdType *desc, DataArrayIdType *descIndx, DataArrayIdType *revDesc, DataArrayIdType *revDescIndx) const; MEDCOUPLING_EXPORT MEDCouplingUMesh *buildDescendingConnectivity2(DataArrayIdType *desc, DataArrayIdType *descIndx, DataArrayIdType *revDesc, DataArrayIdType *revDescIndx) const; @@ -320,6 +321,7 @@ namespace MEDCoupling MCAuto& elts, MCAuto& eltsIndex, std::function sensibilityTo2DQuadraticLinearCellsFunc) const; /// @cond INTERNAL + static void DeleteCellTypeInIndexedArray(const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, MCAuto& arrOut, MCAuto& arrIndxOut); static MEDCouplingUMesh *MergeUMeshesLL(const std::vector& a); typedef mcIdType (*DimM1DescNbrer)(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2); template diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx index b69dd2270..4acfeb560 100755 --- a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx @@ -1313,6 +1313,74 @@ DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCoupling return ret.retn(); } +/*! + * \param [in] arrIn array part of the indexed array pair to be treated + * \param [in] arrIndxIn array index part of the indexed array pair to be treated + * \param [out] arrOut array part of the indexed array pair containing result + * \param [out] arrIndxOut array index part of the indexed array pair containing result + */ +void MEDCouplingUMesh::DeleteCellTypeInIndexedArray(const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, MCAuto& arrOut, MCAuto& arrIndxOut) +{ + if( !arrIn || !arrIn->isAllocated() ) + THROW_IK_EXCEPTION("input array is null or not allocated !"); + if( !arrIndxIn || !arrIndxIn->isAllocated() ) + THROW_IK_EXCEPTION("input indexed array is null or not allocated !"); + arrIn->checkNbOfComps(1,"input array"); arrIndxIn->checkNbOfComps(1,"input indexed array"); + mcIdType arrNbTuples(arrIn->getNumberOfTuples()),arrIndxNbTuples(arrIndxIn->getNumberOfTuples()); + if( arrIndxNbTuples < 1 ) + THROW_IK_EXCEPTION("Indexed array is supposed to have length >= 1 !"); + if( arrNbTuples < arrIndxNbTuples ) + THROW_IK_EXCEPTION("Number of tuples of input array is to low compared to indexed array one !"); + const mcIdType *inArrPtr(arrIn->begin()),*inArrIndxPtr(arrIndxIn->begin()); + if( *inArrIndxPtr != 0 ) + THROW_IK_EXCEPTION("First value of indexed array must be 0 !"); + arrOut = DataArrayIdType::New(); arrOut->alloc(arrNbTuples-arrIndxNbTuples+1,1); + arrIndxOut = DataArrayIdType::New(); arrIndxOut->alloc(arrIndxNbTuples,1); + bool presenceOfPolyh = false; + { + mcIdType *outArrPtr(arrOut->getPointer()),*outArrIndxPtr(arrIndxOut->getPointer()); + *outArrIndxPtr++ = 0; + for( mcIdType i = 0 ; i < arrIndxNbTuples - 1 ; ++i ) + { + mcIdType startPos(*inArrIndxPtr++); + mcIdType endPos(*inArrIndxPtr); + if(inArrPtr[startPos] == INTERP_KERNEL::NORM_POLYHED) + { + presenceOfPolyh = true; + break; + } + outArrPtr = std::copy(inArrPtr+startPos+1,inArrPtr+endPos,outArrPtr); + *outArrIndxPtr++ = endPos - i - 1; + } + } + if(!presenceOfPolyh) + return ; + // + { + arrOut = DataArrayIdType::New(); arrOut->alloc(0,1); + inArrIndxPtr = arrIndxIn->begin(); + mcIdType *outArrIndxPtr = arrIndxOut->getPointer(); + *outArrIndxPtr++ = 0; + for( mcIdType i = 0 ; i < arrIndxNbTuples - 1 ; ++i,++outArrIndxPtr ) + { + mcIdType startPos(*inArrIndxPtr++); + mcIdType endPos(*inArrIndxPtr); + if(inArrPtr[startPos] != INTERP_KERNEL::NORM_POLYHED) + { + arrOut->insertAtTheEnd(inArrPtr+startPos+1,inArrPtr+endPos); + outArrIndxPtr[0] = outArrIndxPtr[-1] + (endPos-startPos-1); + } + else + { + std::set s(inArrPtr+startPos+1,inArrPtr+endPos); + s.erase( -1 ); + arrOut->insertAtTheEnd(s.begin(),s.end()); + outArrIndxPtr[0] = outArrIndxPtr[-1] + s.size(); + } + } + } +} + MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector& a) { if(a.empty()) diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index 4f3a298f1..2a75dce56 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -449,9 +449,9 @@ class MEDCouplingBasicsTest7(unittest.TestCase): # non regression test in python wrapping rg=DataArrayInt64([0,10,29,56,75,102,121,148,167,194,213,240,259,286,305,332,351,378,397,424,443,470,489,516]) a,b,c=DataArrayInt64([75]).splitByValueRange(rg) - assert(a.isEqual(DataArrayInt64([4]))) - assert(b.isEqual(DataArrayInt64([0]))) - assert(c.isEqual(DataArrayInt64([4]))) + self.assertTrue(a.isEqual(DataArrayInt64([4]))) + self.assertTrue(b.isEqual(DataArrayInt64([0]))) + self.assertTrue(c.isEqual(DataArrayInt64([4]))) pass def testDAIBuildExplicitArrByRanges1(self): @@ -1247,6 +1247,65 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertTrue( disc.getMeasureField(tri,True).getArray().isEqual( tri.getMeasureField(True).getArray(), 1e-10) ) pass + def testUMeshExplodeMeshTo(self): + """ + [EDF27988] : implementation of reduceToCells implies implementation of MEDCouplingUMesh.explodeMeshTo + """ + arr = DataArrayDouble(5) ; arr.iota() + m = MEDCouplingCMesh() ; m.setCoords(arr,arr,arr) + m = m.buildUnstructured() + m1 = m[::2] ; m2 = m[1::2] + m1.simplexize(PLANAR_FACE_5) + m = MEDCouplingUMesh.MergeUMeshesOnSameCoords([m1,m2]) + mE1 = m.explodeMeshTo(-1) + ref1 = m.buildDescendingConnectivity() + mE2 = m.explodeMeshTo(-2) + ref2 = m.explode3DMeshTo1D() + mE3 = m.explodeMeshTo(-3) + self.assertTrue( len(mE1) ==5 ) + self.assertTrue( mE1[0].getNodalConnectivity().isEqual(ref1[0].getNodalConnectivity()) ) + self.assertTrue( mE1[0].getNodalConnectivityIndex().isEqual(ref1[0].getNodalConnectivityIndex()) ) + self.assertTrue( mE1[0].getCoords().getHiddenCppPointer() == m.getCoords().getHiddenCppPointer() ) + for i in range(1,5): + self.assertTrue( mE1[i].isEqual(ref1[i]) ) + # + self.assertTrue( len(mE2) ==5 ) + self.assertTrue( mE2[0].getNodalConnectivity().isEqual(ref2[0].getNodalConnectivity()) ) + self.assertTrue( mE2[0].getNodalConnectivityIndex().isEqual(ref2[0].getNodalConnectivityIndex()) ) + self.assertTrue( mE2[0].getCoords().getHiddenCppPointer() == m.getCoords().getHiddenCppPointer() ) + for i in range(1,5): + self.assertTrue( mE2[i].isEqual(ref2[i]) ) + # + self.assertTrue( mE3[0].getMeshDimension() == 0 ) + self.assertTrue( mE3[0].getNumberOfCells() == mE3[0].getNumberOfNodes() ) + a,b = m.getReverseNodalConnectivity() + self.assertTrue( mE3[3].isEqual(a) and mE3[4].isEqual(b) ) + ref3_2 = (m.getNodalConnectivityIndex().deltaShiftIndex()-1) ; ref3_2.computeOffsetsFull() + self.assertTrue( ref3_2.isEqual(mE3[2]) ) + tmp = m.getNodalConnectivityIndex().deepCopy() ; tmp.popBackSilent() ; tmp = tmp.buildComplement( len(m.getNodalConnectivity()) ) ; ref3_1 = m.getNodalConnectivity()[tmp] + self.assertTrue( ref3_1.isEqual(mE3[1]) ) + # + cellsInPolyh = [37,160] + polyh = m[cellsInPolyh] + polyh.convertAllToPoly() + m[cellsInPolyh] = polyh + pE3 = m.explodeMeshTo(-3) + self.assertTrue( pE3[0].getMeshDimension() == 0 ) + self.assertTrue( pE3[0].getNumberOfCells() == pE3[0].getNumberOfNodes() ) + a,b = m.getReverseNodalConnectivity() + self.assertTrue( pE3[3].isEqual(a) and pE3[4].isEqual(b) ) + self.assertTrue( pE3[2].isEqual(mE3[2]) ) # indexed arrays are the same + + ref_a,ref_b = DataArrayInt.ExtractFromIndexedArrays( DataArrayInt(cellsInPolyh).buildComplement(m.getNumberOfCells()), mE3[1], mE3[2] ) + a,b = DataArrayInt.ExtractFromIndexedArrays( DataArrayInt(cellsInPolyh).buildComplement(m.getNumberOfCells()), pE3[1], pE3[2] ) + self.assertTrue( ref_a.isEqual(a) ) + self.assertTrue( ref_b.isEqual(b) ) + for cell in cellsInPolyh: + ref_c,ref_d = DataArrayInt.ExtractFromIndexedArrays( cell, mE3[1], mE3[2] ) ; ref_c.sort() + c,d = DataArrayInt.ExtractFromIndexedArrays( cell, pE3[1], pE3[2] ) + self.assertTrue( ref_c.isEqual(c) ) + self.assertTrue( ref_d.isEqual(d) ) + pass if __name__ == '__main__': diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 25280d253..448d3e769 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -2734,6 +2734,19 @@ namespace MEDCoupling return ret; } + PyObject *explodeMeshTo(int targetDeltaLevel) const + { + MCAuto desc,descIndx,revDesc,revDescIndx; + MCAuto m=self->explodeMeshTo(targetDeltaLevel,desc,descIndx,revDesc,revDescIndx); + PyObject *ret=PyTuple_New(5); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(m.retn()),SWIGTYPE_p_MEDCoupling__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(desc.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(descIndx.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,3,SWIG_NewPointerObj(SWIG_as_voidptr(revDesc.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,4,SWIG_NewPointerObj(SWIG_as_voidptr(revDescIndx.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + return ret; + } + PyObject *explodeIntoEdges() const { MCAuto desc,descIndex,revDesc,revDescIndx; diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index d00f5ba7a..4de31cc17 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -4534,7 +4534,7 @@ DataArrayIdType *MEDFileUMesh::deduceNodeSubPartFromCellSubPart(const std::map >& extractDef) const { diff --git a/src/MEDLoader/MEDFileMesh.hxx b/src/MEDLoader/MEDFileMesh.hxx index 993625e66..2cdbaa0f7 100644 --- a/src/MEDLoader/MEDFileMesh.hxx +++ b/src/MEDLoader/MEDFileMesh.hxx @@ -368,6 +368,7 @@ MCAuto& coords, MCAuto& partCoords, MCAuto >& extractDef) const; MEDLOADER_EXPORT MEDFileUMesh *extractPart(const std::map >& extractDef) const; + // MCAuto extractPart(int level, DataArrayIdType *keepCells, bool removeOrphanNodes = true) const; // implemented in python MEDLOADER_EXPORT MEDFileUMesh *buildExtrudedMesh(const MEDCouplingUMesh *m1D, int policy) const; MEDLOADER_EXPORT MEDFileUMesh *linearToQuadratic(int conversionType=0, double eps=1e-12) const; MEDLOADER_EXPORT MEDFileUMesh *quadraticToLinear(double eps=1e-12) const; diff --git a/src/MEDLoader/Swig/CMakeLists.txt b/src/MEDLoader/Swig/CMakeLists.txt index d718bafb0..0cb6bc26d 100644 --- a/src/MEDLoader/Swig/CMakeLists.txt +++ b/src/MEDLoader/Swig/CMakeLists.txt @@ -96,7 +96,7 @@ INSTALL(FILES MEDLoader.i MEDLoaderTypemaps.i MEDLoaderCommon.i MEDLoaderFinaliz SALOME_INSTALL_SCRIPTS(${CMAKE_CURRENT_BINARY_DIR}/MEDLoader.py ${MEDCOUPLING_INSTALL_PYTHON} EXTRA_DPYS "${SWIG_MODULE_MEDLoader_REAL_NAME}") -INSTALL(FILES MEDLoaderSplitter.py DESTINATION ${MEDCOUPLING_INSTALL_PYTHON}) +INSTALL(FILES MEDLoaderSplitter.py MEDLoaderFinalize.py DESTINATION ${MEDCOUPLING_INSTALL_PYTHON}) INSTALL(FILES ${ALL_TESTS} MEDLoaderDataForTest.py MEDLoaderTest1.py MEDLoaderTest2.py MEDLoaderTest3.py DESTINATION ${MEDCOUPLING_INSTALL_SCRIPT_SCRIPTS}) diff --git a/src/MEDLoader/Swig/MEDLoaderFinalize.i b/src/MEDLoader/Swig/MEDLoaderFinalize.i index c8fab8b86..84eecd379 100644 --- a/src/MEDLoader/Swig/MEDLoaderFinalize.i +++ b/src/MEDLoader/Swig/MEDLoaderFinalize.i @@ -19,6 +19,9 @@ // Author : Anthony Geay (EDF R&D) %pythoncode %{ +import MEDLoaderFinalize +MEDFileUMesh.reduceToCells = MEDLoaderFinalize.MEDFileUMeshReduceToCells +del MEDLoaderFinalize MEDFileMeshesIterator.__next__ = MEDFileMeshesIterator.next MEDFileAnyTypeFieldMultiTSIterator.__next__ = MEDFileAnyTypeFieldMultiTSIterator.next MEDFileFieldsIterator.__next__ = MEDFileFieldsIterator.next diff --git a/src/MEDLoader/Swig/MEDLoaderFinalize.py b/src/MEDLoader/Swig/MEDLoaderFinalize.py new file mode 100644 index 000000000..c5689dbbd --- /dev/null +++ b/src/MEDLoader/Swig/MEDLoaderFinalize.py @@ -0,0 +1,66 @@ +# -*- coding: iso-8859-1 -*- +# Copyright (C) 2023 CEA, EDF +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# + +def MEDFileUMeshReduceToCells(self, level, keepCells, removeOrphanNodes=True): + """ + Method returning a new MEDFileUMesh, restriction of self to level and keepCell cells at this level. + This method also + + :param level: Specifies the top level of the returned MEDFileUMesh expected + :param keepCells: A DataArrayInt specifying cell ids at level level of self + :param removeOrphanNodes: Specifies if orphan nodes should be removed at the end + + see also MEDFileUMesh.extractPart + """ + import MEDLoader as ml + subLevs = [l for l in self.getNonEmptyLevels() if l<=level] + subMeshes = [self[lev] for lev in subLevs] + allFamilyFields = [self.getFamilyFieldAtLevel(lev) for lev in subLevs] + allRefMesh = subMeshes[0] + refMesh = allRefMesh[keepCells] + + mmOut = ml.MEDFileUMesh() + # level 0 + mmOut[0] = refMesh + mmOut.setFamilyFieldArr(0,allFamilyFields[0][keepCells]) + + # subLevels + for curLev,meshLev,famFieldLev in zip(subLevs[1:],subMeshes[1:],allFamilyFields[1:]): + allMeshLev,d,di, rd,rdi = allRefMesh.explodeMeshTo( curLev-level ) + a,b = allMeshLev.areCellsIncludedIn(meshLev,2) + if not a: + raise RuntimeError("Error in mesh {}") + dlev,dlevi = ml.DataArrayInt.ExtractFromIndexedArrays( keepCells, d,di ) + dlev2 = dlev.buildUniqueNotSorted() + cellsToKeepLev = ml.DataArrayInt.BuildIntersection([dlev2,b]) + cellsToKeepLev = b.indicesOfSubPart(cellsToKeepLev) + cellsToKeepLev.sort() + mmOut[curLev] = meshLev[cellsToKeepLev] + mmOut.setFamilyFieldArr(curLev,famFieldLev[cellsToKeepLev]) + + allFamNodes = mmOut.getFamilyFieldAtLevel(1) + if allFamNodes: + mmOut.setFamilyFieldArr(1,allFamNodes[:]) + + if removeOrphanNodes: + mmOut.zipCoords() + + mmOut.copyFamGrpMapsFrom(self) + return mmOut diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index d8afdb960..80a1b87ab 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -7132,6 +7132,65 @@ class MEDLoaderTest3(unittest.TestCase): pass + def testUMeshReduceToCells(self): + """ + [EDF27988] : test of MEDFileUMesh.reduceToCells method + """ + arr = DataArrayDouble(4) ; arr.iota() + arrZ = DataArrayDouble(2) ; arrZ.iota() + m = MEDCouplingCMesh() ; m.setCoords(arr,arr,arrZ) + m = m.buildUnstructured() + mm = MEDFileUMesh() + mm[0] = m + m1 = m.buildDescendingConnectivity()[0] + bas1 = DataArrayInt([0,6,11,16,21,25,29,34,38]) + haut1 = DataArrayInt([1,7,12,17,22,26,30,35,39]) + lat1 = DataArrayInt([8,9,23,36]) + allcellsInMM = DataArrayInt.Aggregate([bas1,haut1,lat1]) + allcellsInMM.sort() + m1InMM = m1[allcellsInMM] + mm[-1] = m1InMM + bas1 = allcellsInMM.findIdForEach(bas1) ; bas1.setName("bas1") + haut1 = allcellsInMM.findIdForEach(haut1) ; haut1.setName("haut1") + lat1 = allcellsInMM.findIdForEach(lat1) ; lat1.setName("lat1") + m2 = m1InMM.buildDescendingConnectivity()[0] + lat2 = DataArrayInt( [ 14, 15, 16, 17, 34, 35, 50, 51] ) + allCellsInMM2 = DataArrayInt( [4, 5, 6, 7, 11, 12, 13, 14, 15, 16, 17, 21, 22, 23, 27, 28, 29, 32, 33, 34, 35, 38, 39, 43, 44, 45, 48, 49, 50, 51, 54, 55] ) + lat2 = allCellsInMM2.findIdForEach(lat2) ; lat2.setName("lat2") + m2InMM = m2[allCellsInMM2] + mm[-2] = m2InMM + mm.setName("mesh") + ## + bas1_id = -1 ; haut1_id = -2 ; lat1_id = -3 ; lat2_id = -4 + fam0 = DataArrayInt( m.getNumberOfCells() ) ; fam0.iota() + mm.setFamilyFieldArr(0,fam0) + fam1 = DataArrayInt( m1InMM.getNumberOfCells() ) ; fam1[:] = 0 ; fam1[bas1] = bas1_id ; fam1[haut1] = haut1_id ; fam1[lat1] = lat1_id + mm.setFamilyFieldArr(-1,fam1) + fam2 = DataArrayInt( m2InMM.getNumberOfCells() ) ; fam2[:] = 0 ; fam2[lat2] = lat2_id + mm.setFamilyFieldArr(-2,fam2) + for grp,famid in [(bas1,bas1_id),(haut1,haut1_id),(lat1,lat1_id),(lat2,lat2_id)]: + mm.addFamily(grp.getName(),famid) + mm.setFamiliesIdsOnGroup(grp.getName(),[famid]) + # ze heart of the test + mmr = mm.reduceToCells(0,DataArrayInt([4])) + # assert section + tmp = m[4] ; tmp.zipCoords() + self.assertTrue( mmr[0].isEqual(tmp,1e-12) ) + tmp = MEDCoupling1SGTUMesh( mmr[-1] ) + self.assertTrue( tmp.getCellModelEnum() == NORM_QUAD4 ) + self.assertTrue( tmp.getNodalConnectivity().isEqual( DataArrayInt([0, 4, 5, 1, 1, 0, 2, 3, 5, 7, 6, 4, 2, 6, 7, 3]) ) ) + tmp = MEDCoupling1SGTUMesh( mmr[-2] ) + self.assertTrue( tmp.getCellModelEnum() == NORM_SEG2 ) + self.assertTrue( tmp.getNodalConnectivity().isEqual( DataArrayInt([5, 4, 0, 4, 5, 1, 4, 6, 5, 7, 7, 6, 2, 6, 7, 3]) ) ) + # + self.assertTrue( mmr.getFamilyFieldAtLevel(0).isEqual(DataArrayInt([4])) ) + self.assertTrue( mmr.getFamilyFieldAtLevel(-1).isEqual(DataArrayInt([-3, -1, -2, -3])) ) + self.assertTrue( mmr.getFamilyFieldAtLevel(-2).isEqual(DataArrayInt([0, -4, -4, 0, 0, 0, -4, -4])) ) + # + self.assertTrue( set( mm.getGroupsNames() ) == set( mmr.getGroupsNames() ) ) + for grp in mm.getGroupsNames(): + self.assertTrue( set(mm.getFamiliesIdsOnGroup(grp)) == set(mmr.getFamiliesIdsOnGroup(grp)) ) + pass if __name__ == "__main__": -- 2.39.2