Salome HOME
New method to rearrange nodal conn of HEXA8 for a tetrahedrization that generates...
authorageay <ageay>
Wed, 13 Nov 2013 16:35:45 +0000 (16:35 +0000)
committerageay <ageay>
Wed, 13 Nov 2013 16:35:45 +0000 (16:35 +0000)
src/MEDCoupling/MEDCoupling1GTUMesh.cxx
src/MEDCoupling/MEDCoupling1GTUMesh.hxx
src/MEDCoupling/MEDCouplingMemArray.cxx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 851d637343df10cc6f43464f38cb2f6bd1b79d5a..68c81c8e58d9aed71476a660081e4b3c12f8fc4c 100644 (file)
@@ -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<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName().c_str(),INTERP_KERNEL::NORM_QUAD4));
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c(DataArrayInt::New()); c->alloc(nbHexa8*6*4,1);
+  int *cPtr(c->getPointer());
+  for(int i=0;i<nbHexa8;i++,inConnPtr+=8)
+    {
+      for(int j=0;j<6;j++,cPtr+=4)
+        cm.fillSonCellNodalConnectivity(j,inConnPtr,cPtr);
+    }
+  ret->setCoords(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<int> s(allFacesTmp+4*i,allFacesTmp+4*i+4);
+      bool found(false);
+      for(int j=0;j<6 && !found;j++)
+        {
+          std::set<int> 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<MEDCoupling1SGTUMesh> quads(explodeEachHexa8To6Quad4());//checks that only hexa8
+  int nbHexa8(getNumberOfCells()),*cQuads(quads->getNodalConnectivity()->getPointer());
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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<MEDCouplingUMesh> quadsTmp(quads->buildUnstructured());
+    MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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<nbOfPair;i++)
+      { ptNeigh[ccPtr[2*i+0]]=ccPtr[2*i+1]/6; ptNeigh[ccPtr[2*i+1]]=ccPtr[2*i+0]/6; }
+  }
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
+  std::vector<bool> fetched(nbHexa8,false);
+  std::vector<bool>::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<int> s; s.insert(cellId);//s contains already organized.
+      while(!s.empty())
+        {
+          std::set<int> sNext;
+          for(std::set<int>::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]={
index 087ed529890aedbe8ac2f35d95fdd1e2b22112fb..93ba77249f3eadcd8139b89fb5c5dc2c32325720 100644 (file)
@@ -144,6 +144,8 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT static MEDCoupling1SGTUMesh *Merge1SGTUMeshesOnSameCoords(std::vector<const MEDCoupling1SGTUMesh *>& 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<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const;
     MEDCOUPLING_EXPORT void resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const;
@@ -165,6 +167,8 @@ namespace ParaMEDMEM
     MEDCoupling1DGTUMesh *computeDualMesh2D() const;
   private:
     MEDCouplingAutoRefCountObjectPtr<DataArrayInt> _conn;
+  public:
+    static const int HEXA8_FACE_PAIRS[6];
   };
 
   class MEDCoupling1DGTUMesh : public MEDCoupling1GTUMesh
index ed815f4aff4a62f6a7435cb81038b3dff1fe73e2..4775f960a15ec3ef8cb9eb642ec79e3d5e41bb76 100644 (file)
@@ -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<DataArrayInt> 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".
  */
index 2ba48b592b804217f45e6ad985a10a9790d98f96..6afb8bf8f162066d89ac76a48c41a5af5836cfc6 100644 (file)
@@ -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
 {
index dfc70a7e70ade944a77aaf5b6097ec6551358f56..8924b3e72ca7d252f7a8f109c72fea3c2c425ef8 100644 (file)
@@ -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)