From ed6274a0347cfb6dc4dfa519b0f1c0072749f010 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Thu, 19 Sep 2024 16:20:02 +0200 Subject: [PATCH] [EDF30178] : Tetraedrize on MEDFileUMesh --- src/INTERP_KERNEL/BBTreeDiscrete.txx | 177 ++++++++++++++++ src/MEDCoupling/MEDCouplingMemArray.cxx | 21 -- src/MEDCoupling/MEDCouplingMemArray.hxx | 6 +- src/MEDCoupling/MEDCouplingMemArray.txx | 87 ++++++++ src/MEDCoupling_Swig/DataArrayInt.i | 11 + .../MEDCouplingBasicsTest7.py | 17 ++ src/MEDLoader/Swig/MEDLoaderFinalize.i | 1 + src/MEDLoader/Swig/MEDLoaderFinalize.py | 194 ++++++++++++++++++ src/MEDLoader/Swig/MEDLoaderTest4.py | 56 ++++- 9 files changed, 547 insertions(+), 23 deletions(-) create mode 100644 src/INTERP_KERNEL/BBTreeDiscrete.txx diff --git a/src/INTERP_KERNEL/BBTreeDiscrete.txx b/src/INTERP_KERNEL/BBTreeDiscrete.txx new file mode 100644 index 000000000..2312d4288 --- /dev/null +++ b/src/INTERP_KERNEL/BBTreeDiscrete.txx @@ -0,0 +1,177 @@ +// Copyright (C) 2024 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 +// + +#pragma once + +#include +#include + +#include +#include +#include + + +// DataType is the type of data to locate +// ConnType is the type of IDs returned +template +class BBTreeDiscrete +{ +private: + std::unique_ptr _left; + std::unique_ptr _right; + int _level; + DataType _max_left; + DataType _min_right; + const DataType *_bb; + typename std::vector _elems; + bool _terminal; + ConnType _nbelems; + + static const int MIN_NB_ELEMS=15; + static const int MAX_LEVEL=20; +public: + BBTreeDiscrete() = default; + /*! + Constructor of the bounding box tree + \param bbs pointer to the [x1 y1 x2 y2 ...] array containing the bounding boxes that are to be indexed. + \param elems array to the indices of the elements contained in the BBTreeDiscrete + \param level level in the BBTreeDiscrete recursive structure + \param nbelems nb of elements in the BBTreeDiscrete + */ + BBTreeDiscrete(const DataType* bbs, ConnType *elems, int level, ConnType nbelems): + _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems) + { + if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL) + { + _terminal=true; + + } + DataType median = std::numeric_limits::max(); + { + std::unique_ptr nodes( new DataType [nbelems] ); + _elems.resize(nbelems); + for (ConnType i=0; i(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems); + median = nodes[nbelems/2]; + } + + std::vector new_elems_left; + std::vector new_elems_right; + + new_elems_left.reserve(nbelems/2+1); + new_elems_right.reserve(nbelems/2+1); + DataType max_left = -std::numeric_limits::max(); + DataType min_right= std::numeric_limits::max(); + for (ConnType i=0; i= median) + { + new_elems_right.push_back(elem); + if (valuemax_left) max_left = value; + } + } + _max_left = max_left; + _min_right = min_right; + ConnType *tmp( nullptr ); + if(!new_elems_left.empty()) + tmp = new_elems_left.data(); + _left.reset(new BBTreeDiscrete(bbs, tmp, level+1, (ConnType)new_elems_left.size()) ); + tmp = nullptr; + if(!new_elems_right.empty()) + tmp = new_elems_right.data(); + _right.reset(new BBTreeDiscrete(bbs, tmp, level+1, (ConnType)new_elems_right.size()) ); + + } + + ~BBTreeDiscrete() = default; + + /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb + + \param bb pointer to query bounding box + \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box + */ + void getIntersectingElems(const DataType *bb, std::vector& elems) const + { + // terminal node : return list of elements intersecting bb + if (_terminal) + { + for (ConnType i=0; i<_nbelems; i++) + { + const DataType * const bb_ptr = _bb + _elems[i]*dim; + bool intersects = true; + for (int idim=0; idimgetIntersectingElems(bb, elems); + return; + } + if (value > _max_left) + { + _right->getIntersectingElems(bb,elems); + return; + } + _left->getIntersectingElems(bb,elems); + _right->getIntersectingElems(bb,elems); + } + + ConnType size() + { + if (_terminal) return _nbelems; + return _left->size()+_right->size(); + } +}; diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 84deb1ef3..9d43256dc 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -2533,27 +2533,6 @@ void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& an ang=ret->getAngle(); } -/*! - * Sorts value within every tuple of \a this array. - * \param [in] asc - if \a true, the values are sorted in ascending order, else, - * in descending order. - * \throw If \a this is not allocated. - */ -void DataArrayDouble::sortPerTuple(bool asc) -{ - checkAllocated(); - double *pt=getPointer(); - mcIdType nbOfTuple(getNumberOfTuples()); - std::size_t nbOfComp(getNumberOfComponents()); - if(asc) - for(mcIdType i=0;i()); - declareAsNew(); -} - /*! * Modify all elements of \a this array, so that * an element _x_ becomes \f$ numerator / x \f$. diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index bde9a3816..058502d2b 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -287,6 +287,7 @@ namespace MEDCoupling void renumberInPlace(const mcIdType *old2New); void renumberInPlaceR(const mcIdType *new2Old); void sort(bool asc=true); + void sortPerTuple(bool asc); typename Traits::ArrayType *renumber(const mcIdType *old2New) const; typename Traits::ArrayType *renumberR(const mcIdType *new2Old) const; typename Traits::ArrayType *renumberAndReduce(const mcIdType *old2New, mcIdType newNbOfTuple) const; @@ -505,7 +506,6 @@ namespace MEDCoupling DataArrayDouble *buildEuclidianDistanceDenseMatrix() const; DataArrayDouble *buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const; void asArcOfCircle(double center[2], double& radius, double& ang) const; - void sortPerTuple(bool asc); void applyInv(double numerator); void applyPow(double val); void applyRPow(double val); @@ -585,6 +585,7 @@ namespace MEDCoupling void writeVTK(std::ostream& ofs, mcIdType indent, const std::string& type, const std::string& nameInFile, DataArrayByte *byteArr) const; void transformWithIndArr(const T *indArrBg, const T *indArrEnd); void transformWithIndArr(const MapKeyVal& m); + void findCommonTuples(mcIdType limitTupleId, MCAuto &comm, MCAuto& commIndex) const; DataArrayIdType *findIdsEqual(T val) const; DataArrayIdType *transformWithIndArrR(const T *indArr2Bg, const T *indArrEnd) const; void splitByValueRange(const T *arrBg, const T *arrEnd, @@ -665,6 +666,9 @@ namespace MEDCoupling void modulusEqual(const DataArrayType *other); static DataArrayType *Pow(const DataArrayType *a1, const DataArrayType *a2); void powEqual(const DataArrayType *other); + public: + template + void findCommonTuplesAlg(const T *bbox, mcIdType nbNodes, mcIdType limitNodeId, DataArrayIdType *c, DataArrayIdType *cI) const; public: static DataArrayIdType *FindPermutationFromFirstToSecond(const DataArrayType *ids1, const DataArrayType *ids2); static DataArrayIdType *FindPermutationFromFirstToSecondDuplicate(const DataArrayType *ids1, const DataArrayType *ids2); diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx index b00a16418..46d5ac131 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.txx +++ b/src/MEDCoupling/MEDCouplingMemArray.txx @@ -29,6 +29,7 @@ #include "InterpKernelAutoPtr.hxx" #include "MCAuto.hxx" #include "MEDCouplingMap.txx" +#include "BBTreeDiscrete.txx" #include #include @@ -1093,6 +1094,28 @@ namespace MEDCoupling declareAsNew(); } + /*! + * Sorts value within every tuple of \a this array. + * \param [in] asc - if \a true, the values are sorted in ascending order, else, + * in descending order. + * \throw If \a this is not allocated. + */ + template + void DataArrayTemplate::sortPerTuple(bool asc) + { + this->checkAllocated(); + T *pt( this->getPointer() ); + mcIdType nbOfTuple(this->getNumberOfTuples()); + std::size_t nbOfComp(this->getNumberOfComponents()); + if(asc) + for(mcIdType i=0;i()); + this->declareAsNew(); + } + /*! * Sorts values of the array and put the result in a newly allocated returned array. * This method does not alterate \a this content. @@ -4179,6 +4202,70 @@ struct NotInRange this->declareAsNew(); } + template + template + void DataArrayDiscrete::findCommonTuplesAlg(const T *bbox, mcIdType nbNodes, mcIdType limitNodeId, DataArrayIdType *c, DataArrayIdType *cI) const + { + const T *coordsPtr(this->begin()); + BBTreeDiscrete myTree(bbox,nullptr,0,nbNodes); + std::vector isDone(nbNodes); + for(mcIdType i=0;i intersectingElems; + myTree.getIntersectingElems(coordsPtr+i*SPACEDIM,intersectingElems); + if(intersectingElems.size()>1) + { + std::vector commonNodes; + for(std::vector::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++) + if(*it!=i) + if(*it>=limitNodeId) + { + commonNodes.push_back(*it); + isDone[*it]=true; + } + if(!commonNodes.empty()) + { + cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size())+1); + c->pushBackSilent(i); + c->insertAtTheEnd(commonNodes.begin(),commonNodes.end()); + } + } + } + } + } + + template + void DataArrayDiscrete::findCommonTuples(mcIdType limitTupleId, MCAuto &comm, MCAuto& commIndex) const + { + this->checkAllocated(); + std::size_t nbOfCompo( this->getNumberOfComponents() ); + if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work + throw INTERP_KERNEL::Exception("DataArrayDiscrete::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4."); + + mcIdType nbOfTuples( this->getNumberOfTuples() ); + // + comm = DataArrayIdType::New(); commIndex = DataArrayIdType::New(); comm->alloc(0,1); commIndex->pushBackSilent(0); + switch(nbOfCompo) + { + case 4: + findCommonTuplesAlg<4>(this->begin(),nbOfTuples,limitTupleId,comm,commIndex); + break; + case 3: + findCommonTuplesAlg<3>(this->begin(),nbOfTuples,limitTupleId,comm,commIndex); + break; + case 2: + findCommonTuplesAlg<2>(this->begin(),nbOfTuples,limitTupleId,comm,commIndex); + break; + case 1: + findCommonTuplesAlg<1>(this->begin(),nbOfTuples,limitTupleId,comm,commIndex); + break; + default: + throw INTERP_KERNEL::Exception("DataArrayDiscrete::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !"); + } + } + /*! * Creates a new DataArrayInt containing IDs (indices) of tuples holding value equal to a * given one. The ids are sorted in the ascending order. diff --git a/src/MEDCoupling_Swig/DataArrayInt.i b/src/MEDCoupling_Swig/DataArrayInt.i index 69f6d2847..339c847d3 100644 --- a/src/MEDCoupling_Swig/DataArrayInt.i +++ b/src/MEDCoupling_Swig/DataArrayInt.i @@ -107,6 +107,7 @@ INT getMaxAbsValueInArray() const; INT getMinValueInArray() const; void abs(); + void sortPerTuple(bool asc); ARRAY *computeAbs() const; void applyLin(INT a, INT b, INT compoId); void applyLin(INT a, INT b); @@ -2091,6 +2092,16 @@ } } + PyObject *findCommonTuples(mcIdType limitNodeId=-1) const + { + MCAuto comm,commIndex; + self->findCommonTuples(limitNodeId,comm,commIndex); + PyObject *res = PyList_New(2); + PyList_SetItem(res,0,SWIG_NewPointerObj(SWIG_as_voidptr(comm.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + PyList_SetItem(res,1,SWIG_NewPointerObj(SWIG_as_voidptr(commIndex.retn()),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); + return res; + } + static PyObject *ExtractFromIndexedArrays(PyObject *li, const ARRAY *arrIn, const DataArrayIdType *arrIndxIn) throw(INTERP_KERNEL::Exception) { ARRAY *arrOut=0; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index d83227d2d..f827cbdf3 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1528,5 +1528,22 @@ class MEDCouplingBasicsTest7(unittest.TestCase): assert( b.isEqual( DataArrayInt([4,2,3,5,1,2]) ) ) assert( f.isEqual( DataArrayInt([0,3,6]) ) ) + def testDASortPerTuple1(self): + """ + EDF30178 : useful method DataArrayInt.sortPerTuple + """ + arr = DataArrayInt([(1,2,3),(5,4,6),(9,8,7)]) + arr.sortPerTuple( True ) + self.assertTrue( arr.isEqual( DataArrayInt([(1,2,3),(4,5,6),(7,8,9)]) ) ) + + def testDAIFindCommonTuples(self): + """ + EDF30178 : useful method DataArrayInt.findCommonTuples + """ + arr = DataArrayInt([(1,2,3),(4,5,6),(-1,2,3),(3,2,1),(1,2,3),(6,5,4),(4,5,6),(4,5,6)]) + c,ci = arr.findCommonTuples(2) + self.assertTrue( ci.isEqual( DataArrayInt([0,2,5]) ) ) + self.assertTrue( c.isEqual( DataArrayInt([0,4, 1,6,7]) ) ) + if __name__ == '__main__': unittest.main() diff --git a/src/MEDLoader/Swig/MEDLoaderFinalize.i b/src/MEDLoader/Swig/MEDLoaderFinalize.i index 75669a3c2..3238a963c 100644 --- a/src/MEDLoader/Swig/MEDLoaderFinalize.i +++ b/src/MEDLoader/Swig/MEDLoaderFinalize.i @@ -21,6 +21,7 @@ %pythoncode %{ import MEDLoaderFinalize MEDFileUMesh.reduceToCells = MEDLoaderFinalize.MEDFileUMeshReduceToCells +MEDFileUMesh.tetrahedrize = MEDLoaderFinalize.MEDFileUMeshTetrahedrize MEDFileUMesh.fuseNodesAndCells = MEDLoaderFinalize.MEDFileUMeshFuseNodesAndCells del MEDLoaderFinalize MEDFileMeshesIterator.__next__ = MEDFileMeshesIterator.next diff --git a/src/MEDLoader/Swig/MEDLoaderFinalize.py b/src/MEDLoader/Swig/MEDLoaderFinalize.py index ced823319..fa6de1564 100644 --- a/src/MEDLoader/Swig/MEDLoaderFinalize.py +++ b/src/MEDLoader/Swig/MEDLoaderFinalize.py @@ -99,6 +99,200 @@ def MEDFileUMeshFuseNodesAndCells(self, compType = 2 , eps = 1e-6, logLev = logg mmOut.setFamilyFieldArr(1, famsMergedNode) return mmOut +def MEDFileUMeshTetrahedrize(self, splitType, logLev = logging.INFO): + """ + [EDF30178] : Method splitting hexa,prisms and underlying quads into resp and underlying triangles + """ + import MEDLoader as ml + def getLogger( level = logging.INFO ): + FORMAT = '%(levelname)s : %(asctime)s : [%(filename)s:%(funcName)s:%(lineno)s] : %(message)s' + logging.basicConfig( format = FORMAT, level = level ) + return logging.getLogger() + logger = getLogger( logLev ) + + def HexaSpliter( splitType ): + """ + :param splitType : see MEDCouplingUMesh.simplexize + """ + m3 = ml.MEDCouplingUMesh("",3) ; m3.allocateCells() ; m3.insertNextCell(ml.NORM_HEXA8,list(range(8))) + m3.simplexize( splitType ) + m3 = ml.MEDCoupling1SGTUMesh( m3 ) + conn = m3.getNodalConnectivity() ; conn.rearrange(4) + return conn.getValuesAsTuple() + + def Penta6Spliter( splitType ): + return [(3,5,4,1),(1,3,5,0),(0,5,1,2)] + + def SplitByType( geoType, splitType ): + m = { ml.NORM_HEXA8 : HexaSpliter, ml.NORM_PENTA6:Penta6Spliter } + return m[ geoType ](splitType) + + def SplitMeshByType( splitType, m0st, famSt = None ): + """ + :param m0st: MEDCoupling1SGTUMesh instance to be split + :param famSt: DataArrayInt storing input family field attached to m0st + """ + conn = m0st.getNodalConnectivity()[:] ; conn.rearrange( m0st.getNumberOfNodesPerCell() ) + geoType = m0st.getCellModelEnum() + subTetra = SplitByType(geoType,splitType) + famOut = None + if famSt: + famOut = famSt.duplicateEachTupleNTimes( len(subTetra) ) + m0stTetras = ml.MEDCoupling1SGTUMesh(m0st.getName(),ml.NORM_TETRA4) + m0stTetras.setCoords( self.getCoords() ) + connTetras = ml.DataArrayInt.Meld([conn[:,elt] for elt in subTetra]) + connTetras.rearrange(1) + m0stTetras.setNodalConnectivity( connTetras ) + return m0stTetras.buildUnstructured(),famOut + + def LocateTwoTrisForEachQuad(quads,tris): + """ + This function locate for each quad in quads the 2 triangles among triangles into tris. + + :param quads: 4 components DataArrayInt storing nodal conn of quad4 + :param tris: 3 components DataArrayInt storing nodal conn containing division of quad4 to locate + """ + from itertools import combinations + quads.sortPerTuple(True) + tris.sortPerTuple(True) + curCompoId = ml.DataArrayInt(len(quads)) ; curCompoId[:] = 0 + res = ml.DataArrayInt(len(quads) * 2 ) ; res[:] = -1 + for elt in combinations(range(4),3): + arr = ml.DataArrayInt.Aggregate([quads[:,elt],tris]) + offset = len(quads) + c,ci = arr.findCommonTuples(offset) + if not ci.deltaShiftIndex().isUniform(2): + raise RuntimeError("Duplication of tris detected should never happen !") + c.rearrange(2) + if not c[:,0].findIdsGreaterOrEqualTo(offset).empty(): + raise RuntimeError("Duplication of tris detected should never happen !") + if not curCompoId[c[:,0]].findIdsGreaterOrEqualTo(2).empty(): + raise RuntimeError("Internal Error : Quad4 is mapped into more than 2 sub cell triangles ! Something is wrong ! Presence of 3D overlapping cells ?") + res[2*c[:,0]+curCompoId[c[:,0]]] = c[:,1] - offset + curCompoId[c[:,0]] += 1 + if not curCompoId.isUniform(2): + raise RuntimeError("It smells very bad ! Impossible to find 2 triangles for some of quadrangles !") + res.rearrange(2) + return res + + def deal3D( mmOut, splitType): + """ + : return : 3D cells in self having a QUAD4 as subcell candidate of spliting + """ + m0 = self[0] + m0s = [ml.MEDCoupling1SGTUMesh(elt) for elt in m0.splitByType()] + fams0 = self.getFamilyFieldAtLevel(0) + outSubMesh = [] + outFams = [] + startCellId = 0 + for m0st in m0s: + endCellId = startCellId + m0st.getNumberOfCells() + famSt = fams0[startCellId:endCellId] + geoType = m0st.getCellModelEnum() + if geoType == ml.NORM_TETRA4: + outSubMesh.append( m0st.buildUnstructured() ) + outFams.append( famSt ) + continue + m0StSplit, famStOut = SplitMeshByType(splitType,m0st,famSt) + outFams.append( famStOut ) + outSubMesh.append( m0StSplit ) + startCellId = endCellId + m0tetra = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords( outSubMesh ) + fam0tetra = ml.DataArrayInt.Aggregate( outFams ) + m0tetra.setDescription( self.getDescription() ) + m0tetra.setName( self.getName() ) + mmOut[0] = m0tetra + mmOut.setFamilyFieldArr(0,fam0tetra) + return ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords( [elt.buildUnstructured() for elt in m0s if elt.getCellModelEnum() != ml.NORM_TETRA4] ) + + def deal2D( mmOut, meshContainingQuadsAsSubCells, splitType): + m1 = self[-1] + m1s = [ml.MEDCoupling1SGTUMesh(elt) for elt in m1.splitByType()] + managed2DTypes = [ml.NORM_TRI3,ml.NORM_QUAD4] + quads4 = [ elt for elt in m1s if elt.getCellModelEnum()==ml.NORM_QUAD4 ] + if not all( [elt.getCellModelEnum() in [ml.NORM_TRI3,ml.NORM_QUAD4] for elt in m1s] ): + typesStr = [ ml.MEDCouplingUMesh.GetReprOfGeometricType( elt.getCellModelEnum() ) for elt in m1s ] + managedTypesStr = [ ml.MEDCouplingUMesh.GetReprOfGeometricType( elt ) for elt in managed2DTypes ] + raise RuntimeError( f"Some geotype in -1 level ( {typesStr} ) are not in managed types ( {managedTypesStr} )" ) + if len( quads4 ) == 1: + quads4 = quads4[0] + pass + logger.debug("Starting to deduce triangulation of quads in -1 level") + logger.debug("Starting to compute sub cells of 3D cells containing QUAD4 as subcell") + two2DCellContainingQuads,_,_,rd,rdi = meshContainingQuadsAsSubCells.buildDescendingConnectivity() + tmp = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords([quads4.buildUnstructured(),two2DCellContainingQuads]) + offset = quads4.getNumberOfCells() + logger.debug("Try to reduce list of 3D cells containing QUAD4 as subcell") + cce,ccei = tmp.findCommonCells(2,offset) + if not ccei.deltaShiftIndex().isUniform(2): + raise RuntimeError("Case of fusable quad4 not managed") + cce.rearrange(2) + if not cce[:,0].findIdsGreaterOrEqualTo( offset ).empty(): + raise RuntimeError("Case of fusable quad4 not managed") + cells3DToKeep,_ = ml.DataArrayInt.ExtractFromIndexedArrays(cce[:,1]-offset,rd,rdi) + cells3DToKeep.sort() + cells3DToKeep = cells3DToKeep.buildUnique() + threedCellsLyingOnQuads = meshContainingQuadsAsSubCells[cells3DToKeep] + threedCellsLyingOnQuads.sortCellsInMEDFileFrmt() + logger.debug("Start to compute the most compact list of tetras") + allSubTetras = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords( [ SplitMeshByType( splitType , ml.MEDCoupling1SGTUMesh( elt ) )[0] for elt in threedCellsLyingOnQuads.splitByType() ] ) + allSubTris = ml.MEDCoupling1SGTUMesh( allSubTetras.buildDescendingConnectivity()[0] ) + cSubTris = allSubTris.getNodalConnectivity()[:] + cSubTris.rearrange(3) + cQuads4 = quads4.getNodalConnectivity()[:] ; cQuads4.rearrange(4) + logger.debug("Start to find the right split of input quads to respect conformity with previous 3D splitting") + res = LocateTwoTrisForEachQuad(cQuads4,cSubTris) + + m1Out = ml.MEDCoupling1SGTUMesh(self.getName(),ml.NORM_TRI3) + m1Out.copyTinyInfoFrom(self[0]) + m1Out.setCoords( self.getCoords() ) + res.rearrange(1) + cSubTris[res] + connOut = cSubTris[res] + connOut.rearrange(1) + m1Out.setNodalConnectivity( connOut ) + m1Out = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords( [ elt.buildUnstructured() for elt in m1s if elt.getCellModelEnum()==ml.NORM_TRI3 ] + [ m1Out.buildUnstructured() ] ) + m1Out.copyTinyInfoFrom(self[0]) + mmOut[-1] = m1Out + famM1 = self.getFamilyFieldAtLevel(-1) + if famM1: + outFams = [] + logger.debug("Start dealing families of 2D cells") + startCellId = 0 + for m1st in m1s: + endCellId = startCellId + m1st.getNumberOfCells() + famSt = famM1[startCellId:endCellId] + startCellId = endCellId + geoType = m1st.getCellModelEnum() + if geoType == ml.NORM_TRI3: + outFams.append( famSt ) + elif geoType == ml.NORM_QUAD4: + outFams.append( famSt.duplicateEachTupleNTimes( 2 ) ) + else: + raise RuntimeError("Not managed geo type !") + mmOut.setFamilyFieldArr( -1, ml.DataArrayInt.Aggregate(outFams) ) + # + if self.getMeshDimension() != 3: + raise RuntimeError( f"Expecting mesh with dimension 3 ! Dimension is {self.getMeshDimension()}" ) + mmOut = ml.MEDFileUMesh() + levs = self.getNonEmptyLevels() + logger.info("Treating 3D level") + meshContainingQuadsAsSubCells = deal3D( mmOut, splitType) + if -1 in levs: + deal2D( mmOut, meshContainingQuadsAsSubCells, splitType) + # dealing remaining levs not impacting by tetrahedrization + for remainingLev in [elt for elt in self.getNonEmptyLevels() if elt not in [0,-1]]: + logger.debug(f"Dealing with level {remainingLev}") + mLev = self[ remainingLev ] + mmOut[ remainingLev ] = mLev + famField = self.getFamilyFieldAtLevel( remainingLev ) + if famField: + mmOut.setFamilyFieldArr(remainingLev,famField) + # + mmOut.copyFamGrpMapsFrom( self ) + return mmOut + + def MEDFileUMeshReduceToCells(self, level, keepCells, removeOrphanNodes=True): """ Method returning a new MEDFileUMesh, restriction of self to level and keepCell cells at this level. diff --git a/src/MEDLoader/Swig/MEDLoaderTest4.py b/src/MEDLoader/Swig/MEDLoaderTest4.py index 27a7bf265..e465fff54 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest4.py +++ b/src/MEDLoader/Swig/MEDLoaderTest4.py @@ -5823,7 +5823,61 @@ class MEDLoaderTest4(unittest.TestCase): self.assertTrue( mmOut.getGroupArr(1,"RedNode").isEqualWithoutConsideringStr( DataArrayInt([6, 7, 10, 11]) ) ) self.assertTrue( mmOut.getGroupArr(1,"GreenNode").isEqualWithoutConsideringStr( DataArrayInt([7, 11, 19, 22]) ) ) self.assertEqual( len( mmOut.getFamilyFieldAtLevel(1).findIdsGreaterOrEqualTo(0) ) , 28 ) - self.assertEqual( len( mmOut.getFamilyFieldAtLevel(-1).findIdsLowerOrEqualTo(0) ) , 45 ) + self.assertEqual( len( mmOut.getFamilyFieldAtLevel(-1).findIdsLowerOrEqualTo(0) ) , 45 ) + + def test50(self): + """ + [EDF30178] : test of MEDFileUMesh.tetrahedrize + """ + import logging + arr = DataArrayDouble([(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)]) + m = MEDCouplingUMesh("mesh",2) + m.setCoords( arr ) + m.allocateCells() + m.insertNextCell(NORM_TRI3,[1,2,4]) + m.insertNextCell(NORM_TRI3,[2,4,5]) + m.insertNextCell(NORM_QUAD4,[0,1,4,3]) + m.insertNextCell(NORM_QUAD4,[3,4,7,6]) + m.insertNextCell(NORM_QUAD4,[4,5,8,7]) + m.changeSpaceDimension(3,0.) + arr1D = DataArrayDouble([0,1,2]) + m1D = MEDCouplingCMesh() ; m1D.setCoords(arr1D) ; m1D = m1D.buildUnstructured() + m1D.changeSpaceDimension(3,0.) + m1D.setCoords( m1D.getCoords()[:,[1,2,0]] ) + m3D = m.buildExtrudedMesh(m1D,0) + m3D.setName( m.getName() ) + m3D.sortCellsInMEDFileFrmt() + m.setCoords( m3D.getCoords() ) + mm = MEDFileUMesh() + mm[0] = m3D + mm[-1] = m + grpTri3 = DataArrayInt([0,1]) ; grpTri3.setName("tri3") + grpQuad4 = DataArrayInt([2,3,4]) ; grpQuad4.setName("quad4") + mm.addGroup(-1,grpTri3) + mm.addGroup(-1,grpQuad4) + grpPrism = DataArrayInt([0,1,2,3]) ; grpPrism.setName("penta") + grpHexa = DataArrayInt([4,5,6,7,8,9]) ; grpHexa.setName("hexa") + mm.addGroup(0,grpPrism) + mm.addGroup(0,grpHexa) + splitType = PLANAR_FACE_6 + mmOut = mm.tetrahedrize(splitType,logging.WARNING) # <- hot point is here + self.assertAlmostEqual( mmOut[0].getMeasureField(True).getArray().accumulate()[0], 8.0, 10 ) + self.assertAlmostEqual( mmOut[-1].getMeasureField(True).getArray().accumulate()[0], 4.0, 10 ) + types0 = mmOut[0].splitByType() + self.assertEqual( len(types0) , 1 ) + types0 = MEDCoupling1SGTUMesh(types0[0]) + self.assertEqual( types0.getCellModelEnum() , NORM_TETRA4 ) + types1 = mmOut[-1].splitByType() + self.assertEqual( len(types1) , 1 ) + types1 = MEDCoupling1SGTUMesh(types1[0]) + self.assertEqual( types1.getCellModelEnum() , NORM_TRI3 ) + self.assertAlmostEqual( mmOut.getGroup(-1,"tri3").getMeasureField(True).getArray().accumulate()[0], 1.0, 10 ) + self.assertAlmostEqual( mmOut.getGroup(-1,"quad4").getMeasureField(True).getArray().accumulate()[0], 3.0, 10 ) + self.assertAlmostEqual( mmOut.getGroup(0,"penta").getMeasureField(True).getArray().accumulate()[0], 2.0, 10 ) + self.assertAlmostEqual( mmOut.getGroup(0,"hexa").getMeasureField(True).getArray().accumulate()[0], 6.0, 10 ) + self.assertEqual( mmOut.getCoords().getHiddenCppPointer() , mm.getCoords().getHiddenCppPointer() ) + self.assertTrue( mm.getFamilyFieldAtLevel(0).getDifferentValues().isEqual( mmOut.getFamilyFieldAtLevel(0).getDifferentValues() ) ) + self.assertTrue( mm.getFamilyFieldAtLevel(-1).getDifferentValues().isEqual( mmOut.getFamilyFieldAtLevel(-1).getDifferentValues() ) ) pass -- 2.39.2