]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Fix for duplicateNodesOnM1Group() so that non connex cracking
authorabn <adrien.bruneton@cea.fr>
Mon, 14 Sep 2015 12:12:18 +0000 (14:12 +0200)
committerabn <adrien.bruneton@cea.fr>
Mon, 14 Sep 2015 12:12:18 +0000 (14:12 +0200)
are correctly handled.

src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDLoader/MEDFileMesh.cxx
src/MEDLoader/Swig/MEDLoaderTest.py
src/MEDLoader/Swig/MEDLoaderTest2.py
src/MEDLoader/Swig/MEDLoaderTest3.py

index 4efb00a5e11a6443b3e246f9099efe8c6b26e3a3..887ea5ecada141aff69c0a4086deb1c5d7270958 100644 (file)
@@ -1396,9 +1396,9 @@ void MEDCouplingUMesh::simplifyPolyhedra(double eps)
 }
 
 /*!
- * This method returns all node ids used in \b this. The data array returned has to be dealt by the caller.
- * The returned node ids are sortes ascendingly. This method is closed to MEDCouplingUMesh::getNodeIdsInUse except
- * the format of returned DataArrayInt instance.
+ * This method returns all node ids used in the connectivity of \b this. The data array returned has to be dealt by the caller.
+ * The returned node ids are sorted ascendingly. This method is close to MEDCouplingUMesh::getNodeIdsInUse except
+ * the format of the returned DataArrayInt instance.
  * 
  * \return a newly allocated DataArrayInt sorted ascendingly of fetched node ids.
  * \sa MEDCouplingUMesh::getNodeIdsInUse, areAllNodesFetched
@@ -2325,13 +2325,13 @@ DataArrayInt *MEDCouplingUMesh::findCellIdsOnBoundary() const
 }
 
 /*!
- * This method find in \b this cells ids that lie on mesh \b otherDimM1OnSameCoords.
+ * This method finds in \b this the cell ids that lie on mesh \b otherDimM1OnSameCoords.
  * \b this and \b otherDimM1OnSameCoords have to lie on the same coordinate array pointer. The coherency of that coords array with connectivity
  * of \b this and \b otherDimM1OnSameCoords is not important here because this method works only on connectivity.
  * this->getMeshDimension() - 1 must be equal to otherDimM1OnSameCoords.getMeshDimension()
  *
- * s0 is the cells ids set in \b this lying on at least one node in fetched nodes in \b otherDimM1OnSameCoords.
- * This method method returns cells ids set s = s1 + s2 where :
+ * s0 is the cell ids set in \b this lying on at least one node in the fetched nodes in \b otherDimM1OnSameCoords.
+ * This method returns the cells ids set s = s1 + s2 where:
  * 
  *  - s1 are cells ids in \b this whose dim-1 constituent equals a cell in \b otherDimM1OnSameCoords.
  *  - s2 are cells ids in \b s0 - \b s1 whose at least two neighbors are in s1.
@@ -2339,8 +2339,8 @@ DataArrayInt *MEDCouplingUMesh::findCellIdsOnBoundary() const
  * \throw if \b otherDimM1OnSameCoords is not part of constituent of \b this, or if coordinate pointer of \b this and \b otherDimM1OnSameCoords
  *        are not same, or if this->getMeshDimension()-1!=otherDimM1OnSameCoords.getMeshDimension()
  *
- * \param [out] cellIdsRk0 a newly allocated array containing cells ids in \b this containg s0 in above algorithm.
- * \param [out] cellIdsRk1 a newly allocated array containing cells ids of s1+s2 \b into \b cellIdsRk0 subset. To get absolute ids of s1+s2 simply invoke
+ * \param [out] cellIdsRk0 a newly allocated array containing the cell ids of s0 (which are cell ids of \b this) in the above algorithm.
+ * \param [out] cellIdsRk1 a newly allocated array containing the cell ids of s1+s2 \b into the \b cellIdsRk0 subset. To get the absolute ids of s1+s2, simply invoke
  *              cellIdsRk1->transformWithIndArr(cellIdsRk0->begin(),cellIdsRk0->end());
  */
 void MEDCouplingUMesh::findCellIdsLyingOn(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayInt *&cellIdsRk0, DataArrayInt *&cellIdsRk1) const
@@ -2466,19 +2466,60 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
   cellIdsRk1->transformWithIndArr(cellIdsRk0Auto->begin(),cellIdsRk0Auto->end());
   //
   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m0Part2=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellIdsRk1->begin(),cellIdsRk1->end(),true));
+  int nCells2 = m0Part2->getNumberOfCells();
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc00=DataArrayInt::New(),descI00=DataArrayInt::New(),revDesc00=DataArrayInt::New(),revDescI00=DataArrayInt::New();
   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m01=m0Part2->buildDescendingConnectivity(desc00,descI00,revDesc00,revDescI00);
+  // Neighbor information of the mesh without considering the crack (serves to count how many connex pieces it is made of)
+  DataArrayInt *tmp00=0,*tmp11=0;
+  MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00, tmp00, tmp11);
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> neighInit00(tmp00);
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> neighIInit00(tmp11);
+  // Neighbor information of the mesh WITH the crack (some neighbors are removed):
   DataArrayInt *idsTmp=0;
   bool b=m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids(idsTmp);
   if(!b)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : the given mdim-1 mesh in other is not a constituent of this !");
+  // In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
+  // of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
   MEDCouplingUMesh::RemoveIdsFromIndexedArrays(ids->begin(),ids->end(),desc00,descI00);
   DataArrayInt *tmp0=0,*tmp1=0;
+  // Compute the neighbor of each cell in m0Part2, taking into account the broken link above. Two
+  // cells on either side of the crack (defined by the mesh of low dimension) are not neighbor anymore.
   ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00,tmp0,tmp1);
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> neigh00(tmp0);
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> neighI00(tmp1);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellsToModifyConn0_torenum=MEDCouplingUMesh::ComputeSpreadZoneGradually(neigh00,neighI00);
+
+
+  // For each initial connex part of the sub-mesh (or said differently for each independent crack):
+  int seed = 0, nIter = 0;
+  int nIterMax = nCells2+1; // Safety net for the loop
+  bool * hitCells = new bool[nCells2];
+  for (int iii=0; iii<nCells2; iii++) hitCells[iii] = false;
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellsToModifyConn0_torenum = DataArrayInt::New();
+  cellsToModifyConn0_torenum->alloc(0,1);
+  while (nIter < nIterMax)
+    {
+      if (std::find(hitCells, hitCells+nCells2, false) == hitCells+nCells2)
+        break;
+      // Connex zone without the crack (to compute the next seed really)
+      int dnu;
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connexCheck = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neighInit00,neighIInit00, -1, dnu);
+      int cnt = 0;
+      for (int * ptr = connexCheck->getPointer(); cnt < connexCheck->getNumberOfTuples(); ptr++, cnt++)
+        hitCells[*ptr] = true;
+      // Connex zone WITH the crack (to identify cells lying on either part of the crack)
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> spreadZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neigh00,neighI00, -1, dnu);
+      cellsToModifyConn0_torenum = DataArrayInt::Aggregate(cellsToModifyConn0_torenum, spreadZone, 0);
+      // Compute next seed, i.e. a cell in another connex part:
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> comple = cellsToModifyConn0_torenum->buildComplement(nCells2);
+      seed = comple->getIJ(0,0);
+      nIter++;
+    }
+  delete[] hitCells;
+  if (nIter >= nIterMax)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate(): internal error - too many iterations.");
+
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1);
   cellsToModifyConn0_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end());
   cellsToModifyConn1_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end());
@@ -11021,7 +11062,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx(const int *idsOfSelectBg, c
 /*!
  * This method works on a pair input (\b arrIn, \b arrIndxIn) where \b arr indexes is in \b arrIndxIn.
  * This method expects that these two input arrays come from the output of MEDCouplingUMesh::computeNeighborsOfCells method.
- * This method start from id 0 that will be contained in output DataArrayInt. It searches then all neighbors of id0 regarding arrIn[arrIndxIn[0]:arrIndxIn[0+1]].
+ * This method start from id 0 that will be contained in output DataArrayInt. It searches then all neighbors of id0 looking at arrIn[arrIndxIn[0]:arrIndxIn[0+1]].
  * Then it is repeated recursively until either all ids are fetched or no more ids are reachable step by step.
  * A negative value in \b arrIn means that it is ignored.
  * This method is useful to see if a mesh is contiguous regarding its connectivity. If it is not the case the size of returned array is different from arrIndxIn->getNumberOfTuples()-1.
index b2a610fbbee11fef798cf3fcaff0931ae75d57da..b0956d7ed1fa81e6d26354a345f5f6e2f3bdd6f2 100644 (file)
@@ -3581,7 +3581,8 @@ void MEDFileUMesh::optimizeFamilies()
     _groups.erase(*it);
 }
 
-void MEDFileUMesh::duplicateNodesOnM1Group(const std::string& grpNameM1, DataArrayInt *&nodesDuplicated, DataArrayInt *&cellsModified, DataArrayInt *&cellsNotModified)
+void MEDFileUMesh::duplicateNodesOnM1Group(const std::string& grpNameM1, DataArrayInt *&nodesDuplicated,
+                                           DataArrayInt *&cellsModified, DataArrayInt *&cellsNotModified)
 {
   std::vector<int> levs=getNonEmptyLevels();
   if(std::find(levs.begin(),levs.end(),0)==levs.end() || std::find(levs.begin(),levs.end(),-1)==levs.end())
index b1d386ce0569aa774217404c3452af170a0d4fa2..38ee4afa9671e1d4251aa830e5644ab2fb5b6a13 100644 (file)
@@ -752,4 +752,5 @@ class MEDLoaderTest(unittest.TestCase):
         pass
     pass
 
-unittest.main()
+if __name__ == "__main__":
+  unittest.main()
index 7da9d18fa963d5747be95a3deb827fd81db6839b..06cc645d08f78f66ede23570cf04f0eebf117ffc 100644 (file)
@@ -337,4 +337,5 @@ class MEDLoaderTest(unittest.TestCase):
         pass
     pass
 
-unittest.main()
+if __name__ == "__main__":
+  unittest.main()
index 71fc2494d69d4bfb72595ecc569110ded0cf63d4..2847ba66fc94b760976417e4773ef3942e0ce411 100644 (file)
@@ -1332,6 +1332,48 @@ class MEDLoaderTest(unittest.TestCase):
         mm.write(fname,2)       
         pass
 
+    def testDuplicateNodesOnM1Group3(self):
+        """ Test duplicateNodesOnM1Group() with *non-connex* cracks """
+        fname = "Pyfile73.med"
+        m = MEDCouplingCMesh.New()
+        m.setCoordsAt(0, DataArrayDouble([0.0,1.1,2.3,3.6,5.0]))
+        m.setCoordsAt(1, DataArrayDouble([0.,1.,2.]))
+        m = m.buildUnstructured(); m.setName("simple")
+        m2 = m.buildDescendingConnectivity()[0]
+        m2.setName(m.getName())
+            
+        # A crack in two non connected parts of the mesh:
+        grpSeg = DataArrayInt([2,11]) ; grpSeg.setName("Grp") 
+
+        mm = MEDFileUMesh.New()
+        mm.setMeshAtLevel(0,m)
+        mm.setMeshAtLevel(-1,m2)
+        mm.setGroupsAtLevel(-1,[grpSeg])
+        nodes, cellsMod, cellsNotMod = mm.duplicateNodesOnM1Group("Grp")
+        self.assertEqual([5,9],nodes.getValues());
+        self.assertEqual([0,3],cellsMod.getValues());
+        self.assertEqual([4,7],cellsNotMod.getValues());
+        self.assertEqual(17,mm.getNumberOfNodes())
+        self.assertEqual([2,11],mm.getGroupArr(-1,"Grp").getValues())
+        self.assertEqual([22,23],mm.getGroupArr(-1,"Grp_dup").getValues())
+        ref0=[4, 1, 0, 15, 6, 4, 4, 3, 8, 16]
+        ref1=[4, 6, 5, 10, 11, 4, 9, 8, 13, 14]
+        self.assertEqual(ref0,mm.getMeshAtLevel(0)[[0,3]].getNodalConnectivity().getValues())
+        self.assertEqual(ref1,mm.getMeshAtLevel(0)[[4,7]].getNodalConnectivity().getValues())
+        self.assertRaises(InterpKernelException,mm.getGroup(-1,"Grp_dup").checkGeoEquivalWith,mm.getGroup(-1,"Grp"),2,1e-12);# Grp_dup and Grp are not equal considering connectivity only
+        mm.getGroup(-1,"Grp_dup").checkGeoEquivalWith(mm.getGroup(-1,"Grp"),12,1e-12)# Grp_dup and Grp are equal considering connectivity and coordinates
+
+        refValues=DataArrayDouble([1.1, 1.2, 1.3, 1.4, 1.1, 1.2, 1.3, 1.4])
+        valsToTest=mm.getMeshAtLevel(0).getMeasureField(True).getArray() ; delta=(valsToTest-refValues) ; delta.abs()
+        self.assertTrue(delta.getMaxValue()[0]<1e-12)
+        #
+        mm.getCoords()[-len(nodes):]+=[0.,-0.3]
+        self.assertRaises(InterpKernelException,mm.getGroup(-1,"Grp_dup").checkGeoEquivalWith,mm.getGroup(-1,"Grp"),12,1e-12);
+        refValues2=refValues[:] ; refValues2[0] = 0.935; refValues2[3] = 1.19
+        valsToTest=mm.getMeshAtLevel(0).getMeasureField(True).getArray() ; delta=(valsToTest-refValues2) ; delta.abs()
+        self.assertTrue(delta.getMaxValue()[0]<1e-12)
+        mm.write(fname,2)   
+
     def testBasicConstructors(self):
         fname="Pyfile18.med"
         m=MEDFileMesh.New(fname)
@@ -4555,4 +4597,5 @@ class MEDLoaderTest(unittest.TestCase):
 
     pass
 
-unittest.main()
+if __name__ == "__main__":
+  unittest.main()