--- /dev/null
+// 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 <vector>
+#include <algorithm>
+
+#include <memory>
+#include <limits>
+#include <cmath>
+
+
+// DataType is the type of data to locate
+// ConnType is the type of IDs returned
+template <int dim, class DataType, class ConnType>
+class BBTreeDiscrete
+{
+private:
+ std::unique_ptr<BBTreeDiscrete> _left;
+ std::unique_ptr<BBTreeDiscrete> _right;
+ int _level;
+ DataType _max_left;
+ DataType _min_right;
+ const DataType *_bb;
+ typename std::vector<ConnType> _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<DataType>::max();
+ {
+ std::unique_ptr<DataType[]> nodes( new DataType [nbelems] );
+ _elems.resize(nbelems);
+ for (ConnType i=0; i<nbelems; i++)
+ {
+ ConnType elem;
+ if (elems)
+ elem= elems[i];
+ else
+ elem=i;
+
+ _elems[i]=elem;
+ nodes[i]=bbs[elem*dim+(level%dim)];
+ }
+ if (_terminal) { return; }
+
+ std::nth_element<DataType*>(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems);
+ median = nodes[nbelems/2];
+ }
+
+ std::vector<ConnType> new_elems_left;
+ std::vector<ConnType> new_elems_right;
+
+ new_elems_left.reserve(nbelems/2+1);
+ new_elems_right.reserve(nbelems/2+1);
+ DataType max_left = -std::numeric_limits<DataType>::max();
+ DataType min_right= std::numeric_limits<DataType>::max();
+ for (ConnType i=0; i<nbelems;i++)
+ {
+ ConnType elem;
+ if( elems )
+ elem= elems[i];
+ else
+ elem=i;
+
+ DataType value = bbs[elem*dim+(level%dim)];
+
+ if (value >= median)
+ {
+ new_elems_right.push_back(elem);
+ if (value<min_right) min_right = value;
+ }
+ else
+ {
+ new_elems_left.push_back(elem);
+ if (value>max_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<ConnType>& 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; idim<dim; idim++)
+ {
+ if( bb_ptr[idim] != bb[idim] )
+ intersects=false;
+ }
+ if(intersects)
+ {
+ elems.push_back(_elems[i]);
+ }
+ }
+ return ;
+ }
+
+ //non terminal node
+ DataType value = bb[_level%dim];
+ if (value < _min_right)
+ {
+ _left->getIntersectingElems(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();
+ }
+};
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<nbOfTuple;i++,pt+=nbOfComp)
- std::sort(pt,pt+nbOfComp);
- else
- for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
- std::sort(pt,pt+nbOfComp,std::greater<double>());
- declareAsNew();
-}
-
/*!
* Modify all elements of \a this array, so that
* an element _x_ becomes \f$ numerator / x \f$.
void renumberInPlace(const mcIdType *old2New);
void renumberInPlaceR(const mcIdType *new2Old);
void sort(bool asc=true);
+ void sortPerTuple(bool asc);
typename Traits<T>::ArrayType *renumber(const mcIdType *old2New) const;
typename Traits<T>::ArrayType *renumberR(const mcIdType *new2Old) const;
typename Traits<T>::ArrayType *renumberAndReduce(const mcIdType *old2New, mcIdType newNbOfTuple) const;
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);
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<T, T>& m);
+ void findCommonTuples(mcIdType limitTupleId, MCAuto<DataArrayIdType> &comm, MCAuto<DataArrayIdType>& commIndex) const;
DataArrayIdType *findIdsEqual(T val) const;
DataArrayIdType *transformWithIndArrR(const T *indArr2Bg, const T *indArrEnd) const;
void splitByValueRange(const T *arrBg, const T *arrEnd,
void modulusEqual(const DataArrayType *other);
static DataArrayType *Pow(const DataArrayType *a1, const DataArrayType *a2);
void powEqual(const DataArrayType *other);
+ public:
+ template<int SPACEDIM>
+ 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);
#include "InterpKernelAutoPtr.hxx"
#include "MCAuto.hxx"
#include "MEDCouplingMap.txx"
+#include "BBTreeDiscrete.txx"
#include <set>
#include <sstream>
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<class T>
+ void DataArrayTemplate<T>::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<nbOfTuple;i++,pt+=nbOfComp)
+ std::sort(pt,pt+nbOfComp);
+ else
+ for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
+ std::sort(pt,pt+nbOfComp,std::greater<double>());
+ 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.
this->declareAsNew();
}
+ template<class T>
+ template<int SPACEDIM>
+ void DataArrayDiscrete<T>::findCommonTuplesAlg(const T *bbox, mcIdType nbNodes, mcIdType limitNodeId, DataArrayIdType *c, DataArrayIdType *cI) const
+ {
+ const T *coordsPtr(this->begin());
+ BBTreeDiscrete<SPACEDIM,T,mcIdType> myTree(bbox,nullptr,0,nbNodes);
+ std::vector<bool> isDone(nbNodes);
+ for(mcIdType i=0;i<nbNodes;i++)
+ {
+ if(!isDone[i])
+ {
+ std::vector<mcIdType> intersectingElems;
+ myTree.getIntersectingElems(coordsPtr+i*SPACEDIM,intersectingElems);
+ if(intersectingElems.size()>1)
+ {
+ std::vector<mcIdType> commonNodes;
+ for(std::vector<mcIdType>::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<class T>
+ void DataArrayDiscrete<T>::findCommonTuples(mcIdType limitTupleId, MCAuto<DataArrayIdType> &comm, MCAuto<DataArrayIdType>& 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.
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);
}
}
+ PyObject *findCommonTuples(mcIdType limitNodeId=-1) const
+ {
+ MCAuto<DataArrayIdType> comm,commIndex;
+ self->findCommonTuples(limitNodeId,comm,commIndex);
+ PyObject *res = PyList_New(2);
+ PyList_SetItem(res,0,SWIG_NewPointerObj(SWIG_as_voidptr(comm.retn()),SWIGTITraits<mcIdType>::TI, SWIG_POINTER_OWN | 0 ));
+ PyList_SetItem(res,1,SWIG_NewPointerObj(SWIG_as_voidptr(commIndex.retn()),SWIGTITraits<mcIdType>::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;
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()
%pythoncode %{
import MEDLoaderFinalize
MEDFileUMesh.reduceToCells = MEDLoaderFinalize.MEDFileUMeshReduceToCells
+MEDFileUMesh.tetrahedrize = MEDLoaderFinalize.MEDFileUMeshTetrahedrize
MEDFileUMesh.fuseNodesAndCells = MEDLoaderFinalize.MEDFileUMeshFuseNodesAndCells
del MEDLoaderFinalize
MEDFileMeshesIterator.__next__ = MEDFileMeshesIterator.next
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.
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
