From c5572012b21cff7a170adf4f00b347132a526718 Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 14 Sep 2015 14:12:18 +0200 Subject: [PATCH] Fix for duplicateNodesOnM1Group() so that non connex cracking are correctly handled. --- src/MEDCoupling/MEDCouplingUMesh.cxx | 61 +++++++++++++++++++++++----- src/MEDLoader/MEDFileMesh.cxx | 3 +- src/MEDLoader/Swig/MEDLoaderTest.py | 3 +- src/MEDLoader/Swig/MEDLoaderTest2.py | 3 +- src/MEDLoader/Swig/MEDLoaderTest3.py | 45 +++++++++++++++++++- 5 files changed, 101 insertions(+), 14 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 4efb00a5e..887ea5eca 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -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 m0Part2=static_cast(buildPartOfMySelf(cellIdsRk1->begin(),cellIdsRk1->end(),true)); + int nCells2 = m0Part2->getNumberOfCells(); MEDCouplingAutoRefCountObjectPtr desc00=DataArrayInt::New(),descI00=DataArrayInt::New(),revDesc00=DataArrayInt::New(),revDescI00=DataArrayInt::New(); MEDCouplingAutoRefCountObjectPtr 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 neighInit00(tmp00); + MEDCouplingAutoRefCountObjectPtr 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 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 neigh00(tmp0); MEDCouplingAutoRefCountObjectPtr neighI00(tmp1); - MEDCouplingAutoRefCountObjectPtr 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 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 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 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 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 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. diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index b2a610fbb..b0956d7ed 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -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 levs=getNonEmptyLevels(); if(std::find(levs.begin(),levs.end(),0)==levs.end() || std::find(levs.begin(),levs.end(),-1)==levs.end()) diff --git a/src/MEDLoader/Swig/MEDLoaderTest.py b/src/MEDLoader/Swig/MEDLoaderTest.py index b1d386ce0..38ee4afa9 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest.py +++ b/src/MEDLoader/Swig/MEDLoaderTest.py @@ -752,4 +752,5 @@ class MEDLoaderTest(unittest.TestCase): pass pass -unittest.main() +if __name__ == "__main__": + unittest.main() diff --git a/src/MEDLoader/Swig/MEDLoaderTest2.py b/src/MEDLoader/Swig/MEDLoaderTest2.py index 7da9d18fa..06cc645d0 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest2.py +++ b/src/MEDLoader/Swig/MEDLoaderTest2.py @@ -337,4 +337,5 @@ class MEDLoaderTest(unittest.TestCase): pass pass -unittest.main() +if __name__ == "__main__": + unittest.main() diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 71fc2494d..2847ba66f 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -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() -- 2.39.2