From 42ddcb536033a49bacf71eea7102c43829f3680d Mon Sep 17 00:00:00 2001 From: ageay Date: Wed, 13 Nov 2013 16:35:45 +0000 Subject: [PATCH] New method to rearrange nodal conn of HEXA8 for a tetrahedrization that generates conform mesh --- src/MEDCoupling/MEDCoupling1GTUMesh.cxx | 151 +++++++++++++++++++++++ src/MEDCoupling/MEDCoupling1GTUMesh.hxx | 4 + src/MEDCoupling/MEDCouplingMemArray.cxx | 13 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 7 +- src/MEDCoupling_Swig/MEDCouplingCommon.i | 4 + 5 files changed, 172 insertions(+), 7 deletions(-) diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx index 851d63734..68c81c8e5 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx @@ -26,6 +26,8 @@ using namespace ParaMEDMEM; +const int MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS[6]={0,1,2,4,3,5}; + MEDCoupling1GTUMesh::MEDCoupling1GTUMesh() { } @@ -1632,6 +1634,155 @@ MEDCoupling1GTUMesh *MEDCoupling1SGTUMesh::computeDualMesh() const } } +/*! + * This method explode each NORM_HEXA8 cells in \a this into 6 NORM_QUAD4 cells and put the result into the MEDCoupling1SGTUMesh returned instance. + * + * \return MEDCoupling1SGTUMesh * - a newly allocated instances (to be managed by the caller) storing the result of the explosion. + * \throw If \a this is not a mesh containing only NORM_HEXA8 cells. + * \throw If \a this is not properly allocated. + */ +MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4() const +{ + const INTERP_KERNEL::CellModel& cm(getCellModel()); + if(cm.getEnum()!=INTERP_KERNEL::NORM_HEXA8) + throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4 : this method can be applied only on HEXA8 mesh !"); + int nbHexa8(getNumberOfCells()); + const int *inConnPtr(getNodalConnectivity()->begin()); + MEDCouplingAutoRefCountObjectPtr ret(MEDCoupling1SGTUMesh::New(getName().c_str(),INTERP_KERNEL::NORM_QUAD4)); + MEDCouplingAutoRefCountObjectPtr c(DataArrayInt::New()); c->alloc(nbHexa8*6*4,1); + int *cPtr(c->getPointer()); + for(int i=0;isetCoords(getCoords()); + ret->setNodalConnectivity(c); + return ret.retn(); +} + +/// @cond INTERNAL + +bool UpdateHexa8Cell(int validAxis, int neighId, const int *validConnQuad4NeighSide, int *allFacesNodalConn, int *myNeighbours) +{ + if(myNeighbours[validAxis]==neighId && allFacesNodalConn[4*validAxis+0]==validConnQuad4NeighSide[0]) + return true; + int oldAxis((int)std::distance(myNeighbours,std::find(myNeighbours,myNeighbours+6,neighId))); + std::size_t pos(std::distance(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,std::find(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS+6,oldAxis))); + std::size_t pos0(pos/2),pos1(pos%2); + int oldAxisOpp(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2]); + int oldConn[8],myConn[8]={-1,-1,-1,-1,-1,-1,-1,-1},tmpConn[4],tmpConn2[8]={0,1,2,3,4,5,6,7},edgeConn[2],allFacesTmp[24],neighTmp[6]; + oldConn[0]=allFacesNodalConn[0]; oldConn[1]=allFacesNodalConn[1]; oldConn[2]=allFacesNodalConn[2]; oldConn[3]=allFacesNodalConn[3]; + oldConn[4]=allFacesNodalConn[4]; oldConn[5]=allFacesNodalConn[7]; oldConn[6]=allFacesNodalConn[6]; oldConn[7]=allFacesNodalConn[5]; + const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_HEXA8)); + cm.fillSonCellNodalConnectivity(validAxis,tmpConn2,tmpConn); + for(int i=0;i<4;i++) + myConn[tmpConn[i]]=validConnQuad4NeighSide[(4-i)%4]; + for(int i=0;i<4;i++) + { + int nodeId(myConn[i]!=-1?myConn[i]:myConn[4+i]);//the node id for which the opposite one will be found + bool found(false); + INTERP_KERNEL::NormalizedCellType typeOfSon; + for(int j=0;j<12 && !found;j++) + { + cm.fillSonEdgesNodalConnectivity3D(j,oldConn,-1,edgeConn,typeOfSon); + if(edgeConn[0]==nodeId || edgeConn[1]==nodeId) + { + if(std::find(allFacesNodalConn+4*oldAxisOpp,allFacesNodalConn+4*oldAxisOpp+4,edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0])!=allFacesNodalConn+4*oldAxisOpp+4) + { + if(myConn[i]==-1) + myConn[myConn[i]]=edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0]; + else + myConn[myConn[i]+4]=edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0]; + found=true; + } + } + } + if(!found) + throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error !"); + } + for(int i=0;i<6;i++) + { + cm.fillSonCellNodalConnectivity(i,myConn,allFacesTmp+4*i); + std::set s(allFacesTmp+4*i,allFacesTmp+4*i+4); + bool found(false); + for(int j=0;j<6 && !found;j++) + { + std::set s1(allFacesNodalConn+4*j,allFacesNodalConn+4*j+4); + if(s==s1) + { + neighTmp[i]=myNeighbours[j]; + found=true; + } + } + if(!found) + throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error #2 !"); + } + std::copy(allFacesTmp,allFacesTmp+24,allFacesNodalConn); + std::copy(neighTmp,neighTmp+6,myNeighbours); + return false; +} + +/// @endcond + +/*! + * This method expects the \a this contains NORM_HEXA8 cells only. This method will sort each cells in \a this so that their numbering was + * homogeneous. If it succeeds the result of MEDCouplingUMesh::tetrahedrize will return a conform mesh. + * + * \return DataArrayInt * - a newly allocated array (to be managed by the caller) containing renumbered cell ids. + * + * \throw If \a this is not a mesh containing only NORM_HEXA8 cells. + * \throw If \a this is not properly allocated. + * \sa MEDCouplingUMesh::tetrahedrize, MEDCouplingUMesh::simplexize. + */ +DataArrayInt *MEDCoupling1SGTUMesh::sortHexa8EachOther() +{ + MEDCouplingAutoRefCountObjectPtr quads(explodeEachHexa8To6Quad4());//checks that only hexa8 + int nbHexa8(getNumberOfCells()),*cQuads(quads->getNodalConnectivity()->getPointer()); + MEDCouplingAutoRefCountObjectPtr neighOfQuads(DataArrayInt::New()); neighOfQuads->alloc(nbHexa8*6,1); neighOfQuads->fillWithValue(-1); + int *ptNeigh(neighOfQuads->getPointer()); + {//neighOfQuads tells for each face of each Quad8 which cell (if!=-1) is connected to this face. + MEDCouplingAutoRefCountObjectPtr quadsTmp(quads->buildUnstructured()); + MEDCouplingAutoRefCountObjectPtr ccSafe,cciSafe; + DataArrayInt *cc(0),*cci(0); + quadsTmp->findCommonCells(3,0,cc,cci); + ccSafe=cc; cciSafe=cci; + const int *ccPtr(ccSafe->begin()),nbOfPair(cci->getNumberOfTuples()-1); + for(int i=0;i ret(DataArrayInt::New()); ret->alloc(0,1); + std::vector fetched(nbHexa8,false); + std::vector::iterator it(std::find(fetched.begin(),fetched.end(),true)); + while(it!=fetched.end())//it will turns as time as number of connected zones + { + int cellId((int)std::distance(fetched.begin(),it));//it is the seed of the connected zone. + std::set s; s.insert(cellId);//s contains already organized. + while(!s.empty()) + { + std::set sNext; + for(std::set::const_iterator it0=s.begin();it0!=s.end();it0++) + { + fetched[*it0]=true; + int *myNeighb(ptNeigh+6*(*it0)); + for(int i=0;i<6;i++) + { + std::size_t pos(std::distance(HEXA8_FACE_PAIRS,std::find(HEXA8_FACE_PAIRS,HEXA8_FACE_PAIRS+6,i))); + std::size_t pos0(pos/2),pos1(pos%2); + if(myNeighb[i]!=-1) + { + if(!UpdateHexa8Cell(HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2],*it0,cQuads+6*4*(*it0)+4*i,cQuads+6*4*myNeighb[i],ptNeigh+6*myNeighb[i])) + ret->pushBackSilent(myNeighb[i]); + sNext.insert(myNeighb[i]); + } + } + } + s=sNext; + } + } + return ret.retn(); +} + MEDCoupling1DGTUMesh *MEDCoupling1SGTUMesh::computeDualMesh3D() const { static const int DUAL_TETRA_0[36]={ diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx index 087ed5298..93ba77249 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx @@ -144,6 +144,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT static MEDCoupling1SGTUMesh *Merge1SGTUMeshesOnSameCoords(std::vector& a); MEDCOUPLING_EXPORT MEDCoupling1SGTUMesh *buildSetInstanceFromThis(int spaceDim) const; MEDCOUPLING_EXPORT MEDCoupling1GTUMesh *computeDualMesh() const; + MEDCOUPLING_EXPORT DataArrayInt *sortHexa8EachOther(); + MEDCOUPLING_EXPORT MEDCoupling1SGTUMesh *explodeEachHexa8To6Quad4() const; public://serialization MEDCOUPLING_EXPORT void getTinySerializationInformation(std::vector& tinyInfoD, std::vector& tinyInfo, std::vector& littleStrings) const; MEDCOUPLING_EXPORT void resizeForUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector& littleStrings) const; @@ -165,6 +167,8 @@ namespace ParaMEDMEM MEDCoupling1DGTUMesh *computeDualMesh2D() const; private: MEDCouplingAutoRefCountObjectPtr _conn; + public: + static const int HEXA8_FACE_PAIRS[6]; }; class MEDCoupling1DGTUMesh : public MEDCoupling1GTUMesh diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index ed815f4af..4775f960a 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -1959,7 +1959,7 @@ bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, * [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1 * gives the number of groups of coincident tuples. * \throw If \a this is not allocated. - * \throw If the number of components is not in [1,2,3]. + * \throw If the number of components is not in [1,2,3,4]. * * \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example". * @@ -1970,14 +1970,17 @@ void DataArrayDouble::findCommonTuples(double prec, int limitTupleId, DataArrayI { checkAllocated(); int nbOfCompo=getNumberOfComponents(); - if ((nbOfCompo<1) || (nbOfCompo>3)) //test before work - throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2 or 3."); + if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work + throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4."); int nbOfTuples=getNumberOfTuples(); // MEDCouplingAutoRefCountObjectPtr c(DataArrayInt::New()),cI(DataArrayInt::New()); c->alloc(0,1); cI->pushBackSilent(0); switch(nbOfCompo) { + case 4: + findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI); + break; case 3: findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI); break; @@ -1988,7 +1991,7 @@ void DataArrayDouble::findCommonTuples(double prec, int limitTupleId, DataArrayI findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI); break; default: - throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2 and 3 ! not implemented for other number of components !"); + throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !"); } comm=c.retn(); commIndex=cI.retn(); @@ -2185,7 +2188,7 @@ DataArrayInt *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble * \return DataArrayDouble * - the new instance of DataArrayDouble that the caller * is to delete using decrRef() as it is no more needed. * \throw If \a this is not allocated. - * \throw If the number of components is not in [1,2,3]. + * \throw If the number of components is not in [1,2,3,4]. * * \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example". */ diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 2ba48b592..6afb8bf8f 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -5245,7 +5245,7 @@ void MEDCouplingUMesh::tessellate2DCurve(double eps) * and \a this->getMeshDimension() != 3. * \throw If \a policy is not one of the four discussed above. * \throw If the nodal connectivity of cells is not defined. - * \sa MEDCouplingUMesh::tetrahedrize + * \sa MEDCouplingUMesh::tetrahedrize, MEDCoupling1SGTUMesh::sortHexa8EachOther */ DataArrayInt *MEDCouplingUMesh::simplexize(int policy) { @@ -9391,10 +9391,13 @@ DataArrayInt *MEDCouplingUMesh::ComputeRangesFromTypeDistribution(const std::vec * an id of old cell producing it. The caller is to delete this array using * decrRef() as it is no more needed. * \return MEDCoupling1SGTUMesh * - the mesh containing only INTERP_KERNEL::NORM_TETRA4 cells. + * + * \warning This method operates on each cells in this independantly ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output + * the policy PLANAR_FACE_6 should be used on a mesh sorted with MEDCoupling1SGTUMesh::sortHexa8EachOther. * * \throw If \a this is not a 3D mesh (spaceDim==3 and meshDim==3). * \throw If \a this is not fully constituted with linear 3D cells. - * \sa MEDCouplingUMesh::simplexize + * \sa MEDCouplingUMesh::simplexize, MEDCoupling1SGTUMesh::sortHexa8EachOther */ MEDCoupling1SGTUMesh *MEDCouplingUMesh::tetrahedrize(int policy, DataArrayInt *& n2oCells, int& nbOfAdditionalPoints) const { diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index dfc70a7e7..8924b3e72 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -273,6 +273,8 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCoupling1SGTUMesh::New; %newobject ParaMEDMEM::MEDCoupling1SGTUMesh::buildSetInstanceFromThis; %newobject ParaMEDMEM::MEDCoupling1SGTUMesh::computeDualMesh; +%newobject ParaMEDMEM::MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4; +%newobject ParaMEDMEM::MEDCoupling1SGTUMesh::sortHexa8EachOther; %newobject ParaMEDMEM::MEDCoupling1SGTUMesh::Merge1SGTUMeshes; %newobject ParaMEDMEM::MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords; %newobject ParaMEDMEM::MEDCoupling1DGTUMesh::New; @@ -2574,6 +2576,8 @@ namespace ParaMEDMEM static MEDCoupling1SGTUMesh *Merge1SGTUMeshes(const MEDCoupling1SGTUMesh *mesh1, const MEDCoupling1SGTUMesh *mesh2) throw(INTERP_KERNEL::Exception); MEDCoupling1SGTUMesh *buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception); MEDCoupling1GTUMesh *computeDualMesh() const throw(INTERP_KERNEL::Exception); + MEDCoupling1SGTUMesh *explodeEachHexa8To6Quad4() const throw(INTERP_KERNEL::Exception); + DataArrayInt *sortHexa8EachOther() throw(INTERP_KERNEL::Exception); %extend { MEDCoupling1SGTUMesh(const char *name, INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception) -- 2.30.2