From 8be762fe70bd74f2c86e0576020db8a1a8fcd947 Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 16 Sep 2015 15:55:08 +0200 Subject: [PATCH] Further "improvement" on duplicateNodesOnM1Group(). Cells identification was still buggy. --- src/MEDCoupling/MEDCouplingUMesh.cxx | 79 +++++++++++++--------------- src/MEDLoader/Swig/MEDLoaderTest3.py | 24 ++++----- 2 files changed, 48 insertions(+), 55 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 96d4b34ea..e736e57e5 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2331,16 +2331,14 @@ DataArrayInt *MEDCouplingUMesh::findCellIdsOnBoundary() const * this->getMeshDimension() - 1 must be equal to otherDimM1OnSameCoords.getMeshDimension() * * 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. + * This method also returns the cells ids set s1 which contains the cell ids in \b this for which one of the dim-1 constituent + * equals a cell in \b otherDimM1OnSameCoords. * * \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 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 + * \param [out] cellIdsRk1 a newly allocated array containing the cell ids of s1 \b indexed into the \b cellIdsRk0 subset. To get the absolute ids of s1, simply invoke * cellIdsRk1->transformWithIndArr(cellIdsRk0->begin(),cellIdsRk0->end()); */ void MEDCouplingUMesh::findCellIdsLyingOn(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayInt *&cellIdsRk0, DataArrayInt *&cellIdsRk1) const @@ -2366,23 +2364,10 @@ void MEDCouplingUMesh::findCellIdsLyingOn(const MEDCouplingUMesh& otherDimM1OnSa for(const int *idOther=idsOtherInConsti->begin();idOther!=idsOtherInConsti->end();idOther++) s1.insert(revDescThisPartPtr+revDescIThisPartPtr[*idOther],revDescThisPartPtr+revDescIThisPartPtr[*idOther+1]); MEDCouplingAutoRefCountObjectPtr s1arr_renum1=DataArrayInt::New(); s1arr_renum1->alloc((int)s1.size(),1); std::copy(s1.begin(),s1.end(),s1arr_renum1->getPointer()); - MEDCouplingAutoRefCountObjectPtr s1Comparr_renum1=s1arr_renum1->buildComplement(s0arr->getNumberOfTuples()); - DataArrayInt *neighThisPart=0,*neighIThisPart=0; - ComputeNeighborsOfCellsAdv(descThisPart,descIThisPart,revDescThisPart,revDescIThisPart,neighThisPart,neighIThisPart); - MEDCouplingAutoRefCountObjectPtr neighThisPartAuto(neighThisPart),neighIThisPartAuto(neighIThisPart); - ExtractFromIndexedArrays(s1Comparr_renum1->begin(),s1Comparr_renum1->end(),neighThisPart,neighIThisPart,neighThisPart,neighIThisPart);// reuse of neighThisPart and neighIThisPart - neighThisPartAuto=neighThisPart; neighIThisPartAuto=neighIThisPart; - RemoveIdsFromIndexedArrays(s1Comparr_renum1->begin(),s1Comparr_renum1->end(),neighThisPart,neighIThisPart); - neighThisPartAuto=0; - MEDCouplingAutoRefCountObjectPtr s2_tmp=neighIThisPart->deltaShiftIndex(); - const int li[2]={0,1}; - MEDCouplingAutoRefCountObjectPtr s2_renum2=s2_tmp->getIdsNotEqualList(li,li+2); - s2_renum2->transformWithIndArr(s1Comparr_renum1->begin(),s1Comparr_renum1->end());//s2_renum2==s2_renum1 - MEDCouplingAutoRefCountObjectPtr s_renum1=DataArrayInt::Aggregate(s2_renum2,s1arr_renum1,0); - s_renum1->sort(); - // + s1arr_renum1->sort(); cellIdsRk0=s0arr.retn(); - cellIdsRk1=s_renum1.retn(); + //cellIdsRk1=s_renum1.retn(); + cellIdsRk1=s1arr_renum1.retn(); } /*! @@ -2448,6 +2433,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayInt *& nodeIdsToDuplicate, DataArrayInt *& cellIdsNeededToBeRenum, DataArrayInt *& cellIdsNotModified) const { + typedef MEDCouplingAutoRefCountObjectPtr DAInt; + checkFullyDefined(); otherDimM1OnSameCoords.checkFullyDefined(); if(getCoords()!=otherDimM1OnSameCoords.getCoords()) @@ -2456,28 +2443,28 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : the mesh given in other parameter must have this->getMeshDimension()-1 !"); DataArrayInt *cellIdsRk0=0,*cellIdsRk1=0; findCellIdsLyingOn(otherDimM1OnSameCoords,cellIdsRk0,cellIdsRk1); - MEDCouplingAutoRefCountObjectPtr cellIdsRk0Auto(cellIdsRk0),cellIdsRk1Auto(cellIdsRk1); - MEDCouplingAutoRefCountObjectPtr s0=cellIdsRk1->buildComplement(cellIdsRk0->getNumberOfTuples()); + DAInt cellIdsRk0Auto(cellIdsRk0),cellIdsRk1Auto(cellIdsRk1); + DAInt s0=cellIdsRk1->buildComplement(cellIdsRk0->getNumberOfTuples()); s0->transformWithIndArr(cellIdsRk0Auto->begin(),cellIdsRk0Auto->end()); MEDCouplingAutoRefCountObjectPtr m0Part=static_cast(buildPartOfMySelf(s0->begin(),s0->end(),true)); - MEDCouplingAutoRefCountObjectPtr s1=m0Part->computeFetchedNodeIds(); - MEDCouplingAutoRefCountObjectPtr s2=otherDimM1OnSameCoords.computeFetchedNodeIds(); - MEDCouplingAutoRefCountObjectPtr s3=s2->buildSubstraction(s1); + DAInt s1=m0Part->computeFetchedNodeIds(); + DAInt s2=otherDimM1OnSameCoords.computeFetchedNodeIds(); + DAInt s3=s2->buildSubstraction(s1); 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(); + DAInt 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); + DAInt neighInit00(tmp00); + DAInt 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); + DAInt 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 @@ -2487,38 +2474,44 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On // 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); + DAInt neigh00(tmp0); + DAInt neighI00(tmp1); // 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 - std::vector hitCells(nCells2); - hitCells.assign(nCells2, false); - MEDCouplingAutoRefCountObjectPtr cellsToModifyConn0_torenum = DataArrayInt::New(); + DAInt hitCells = DataArrayInt::New(); hitCells->alloc(nCells2); + hitCells->fillWithValue(-1); + DAInt cellsToModifyConn0_torenum = DataArrayInt::New(); cellsToModifyConn0_torenum->alloc(0,1); while (nIter < nIterMax) { - if (std::find(hitCells.begin(), hitCells.end(), false) == hitCells.end()) + DAInt t = hitCells->getIdsEqual(-1); + if (!t->getNumberOfTuples()) 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); + DAInt 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; + hitCells->setIJ(*ptr,0,1); // 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); + DAInt 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); + // Compute next seed, i.e. a cell in another connex part, which was not covered by the previous iterations + DAInt comple = cellsToModifyConn0_torenum->buildComplement(nCells2); + DAInt nonHitCells = hitCells->getIdsEqual(-1); + DAInt intersec = nonHitCells->buildIntersection(comple); + if (intersec->getNumberOfTuples()) + { seed = intersec->getIJ(0,0); } + else + { break; } nIter++; } if (nIter >= nIterMax) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate(): internal error - too many iterations."); - MEDCouplingAutoRefCountObjectPtr cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1); + DAInt cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1); cellsToModifyConn0_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end()); cellsToModifyConn1_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end()); // diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 2847ba66f..baa860a56 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -1343,34 +1343,34 @@ class MEDLoaderTest(unittest.TestCase): m2.setName(m.getName()) # A crack in two non connected parts of the mesh: - grpSeg = DataArrayInt([2,11]) ; grpSeg.setName("Grp") + grpSeg = DataArrayInt([3,19]) ; 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([1,13],nodes.getValues()); + self.assertEqual([0,6],cellsMod.getValues()); + self.assertEqual([1,7],cellsNotMod.getValues()); self.assertEqual(17,mm.getNumberOfNodes()) - self.assertEqual([2,11],mm.getGroupArr(-1,"Grp").getValues()) + self.assertEqual([3,19],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()) + ref0=[4, 15, 0, 5, 6, 4, 8, 7, 12, 16] + ref1=[4, 2, 1, 6, 7, 4, 9, 8, 13, 14] + self.assertEqual(ref0,mm.getMeshAtLevel(0)[[0,6]].getNodalConnectivity().getValues()) + self.assertEqual(ref1,mm.getMeshAtLevel(0)[[1,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) + self.assertTrue(delta.getMaxValue()[0]<1e-10) # 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() + refValues2=refValues[:] ; refValues2[0] = 1.265; refValues2[6] = 1.105 + valsToTest=mm.getMeshAtLevel(0).getMeasureField(True).getArray() ; delta=(valsToTest-refValues2) ; delta.abs() self.assertTrue(delta.getMaxValue()[0]<1e-12) mm.write(fname,2) -- 2.39.2