From 4643866bcc351da857d08e9abce42aa58b6d6309 Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 11 Jan 2021 14:25:46 +0100 Subject: [PATCH] buildInnerBoundaryAlongM1Group: bug fix for singular point in 3D --- src/MEDCoupling/MEDCouplingUMesh.cxx | 34 ++++++++++++++++----- src/MEDLoader/Swig/MEDLoaderTest3.py | 45 ++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 8 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 53e2b76bc..c20a00b27 100755 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2322,7 +2322,6 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const * \param [out] cellIdsNeededToBeRenum cell ids in \b this in which the renumber of nodes should be performed. * \param [out] cellIdsNotModified cell ids mcIdType \b this that lies on \b otherDimM1OnSameCoords mesh whose connectivity do \b not need to be modified as it is the case for \b cellIdsNeededToBeRenum. * - * \warning This method modifies param \b otherDimM1OnSameCoords (for speed reasons). */ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayIdType *& nodeIdsToDuplicate, DataArrayIdType *& cellIdsNeededToBeRenum, DataArrayIdType *& cellIdsNotModified) const @@ -2339,14 +2338,15 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On // Checking star-shaped M1 group: DAInt dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New(); - MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); + MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); // 2D: a mesh of points, 3D: a mesh of segs DAInt dsi = rdit0->deltaShiftIndex(); - DAInt idsTmp0 = dsi->findIdsNotInRange(-1, 3); + DAInt idsTmp0 = dsi->findIdsNotInRange(-1, 3); // for 2D: if a point is connected to more than 2 segs. For 3D: if a seg is connected to more than two faces. if(idsTmp0->getNumberOfTuples()) throw INTERP_KERNEL::Exception("MEDFileUMesh::buildInnerBoundaryAlongM1Group: group is too complex: some points (or edges) have more than two connected segments (or faces)!"); dt0=0; dit0=0; rdt0=0; rdit0=0; idsTmp0=0; - // Get extreme nodes from the group (they won't be duplicated), ie nodes belonging to boundary cells of M1 + // Get extreme nodes from the group (they won't be duplicated except if they also lie on bound of M0 -- see below), + // ie nodes belonging to the boundary "cells" (might be points) of M1 DAInt xtremIdsM2 = dsi->findIdsEqual(1); dsi = 0; MCUMesh meshM2Part = static_cast(meshM2->buildPartOfMySelf(xtremIdsM2->begin(), xtremIdsM2->end(),true)); DAInt xtrem = meshM2Part->computeFetchedNodeIds(); @@ -2354,15 +2354,17 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New(); MCUMesh m0desc = buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0; dsi = rdit0->deltaShiftIndex(); - DAInt boundSegs = dsi->findIdsEqual(1); // boundary segs/faces of the M0 mesh + DAInt boundSegs = dsi->findIdsEqual(1); dsi = 0; // boundary segs/faces of the M0 mesh MCUMesh m0descSkin = static_cast(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true)); DAInt fNodes = m0descSkin->computeFetchedNodeIds(); - // In 3D, some points on the boundary of M0 still need duplication: + // In 3D, some points on the boundary of M0 will NOT be duplicated (where as in 2D, points on the boundary of M0 are always duplicated) + // Think of a partial (plane) crack in a cube: the points at the tip of the crack and not located inside the volume of the cube are not duplicated + // although they are technically on the skin of the cube. DAInt notDup = 0; if (getMeshDimension() == 3) { DAInt dnu1=DataArrayIdType::New(), dnu2=DataArrayIdType::New(), dnu3=DataArrayIdType::New(), dnu4=DataArrayIdType::New(); - MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4); + MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4); // all segments of the skin of the 3D (M0) mesh dnu1=0;dnu2=0;dnu3=0;dnu4=0; DataArrayIdType * corresp=0; meshM2->areCellsIncludedIn(m0descSkinDesc,2,corresp); @@ -2370,10 +2372,26 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On corresp->decrRef(); if (validIds->getNumberOfTuples()) { + // Build the set of segments which are: in the desc mesh of the skin of the 3D mesh (M0) **and** in the desc mesh of the M1 group: MCUMesh m1IntersecSkin = static_cast(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true)); + // Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above) DAInt notDuplSkin = m1IntersecSkin->findBoundaryNodes(); DAInt fNodes1 = fNodes->buildSubstraction(notDuplSkin); - notDup = xtrem->buildSubstraction(fNodes1); + + // Also, in this (segment) mesh, nodes connected to more than 3 segs should not be dup either (singular points - see testBuildInnerBoundary6()) + dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New(); + MCUMesh meshM2Desc = meshM2->buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0; // a mesh made of node cells + dsi = rdit0->deltaShiftIndex(); + DAInt singPoints = dsi->findIdsNotInRange(-1,4); // points connected to (strictly) more than 3 segments + const mcIdType *cc = meshM2Desc->getNodalConnectivity()->begin(), *ccI = meshM2Desc->getNodalConnectivityIndex()->begin(); + mcIdType * singPointsP = singPoints->rwBegin(); + for (mcIdType j=0; j < singPoints->getNumberOfTuples(); j++) // replace ids in singPoints by real coordinate index (was index of cells in notDuplSkin) + { + mcIdType nodeCellIdx = singPointsP[j]; + singPointsP[j] = cc[ccI[nodeCellIdx]+1]; // +1 to skip type + } + DAInt fNodes2 = fNodes1->buildSubstraction(singPoints); + notDup = xtrem->buildSubstraction(fNodes2); } else notDup = xtrem->buildSubstraction(fNodes); diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 5b89f8706..4773608e4 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -1580,6 +1580,51 @@ class MEDLoaderTest3(unittest.TestCase): m_bis0.checkDeepEquivalOnSameNodesWith(mfu.getMeshAtLevel(-1), 2, 9.9999999) pass + @WriteInTmpDir + def testBuildInnerBoundary6(self): + """ 3D test where the crack has a funny shape with a singular point (i.e. two faces of the M1 group are only connected by one point, not a full segment) + The singular point was wrongly duplicated. + """ + coo = DataArrayDouble([(-1.38778e-17,0.226,0),(-1.38778e-17,-1.38778e-17,0),(0.226,0.226,0),(0.226,-1.38778e-17,0),(0.452,0.226,0),(0.452,-1.38778e-17,0), + (-1.38778e-17,0.452,0),(0.226,0.452,0),(0.452,0.452,0),(-1.38778e-17,0.226,0.25),(0.226,0.226,0.25),(0.226,-1.38778e-17,0.25),(-1.38778e-17,-1.38778e-17,0.25), + (-1.38778e-17,0.226,0.779375),(0.226,0.226,0.779375),(0.226,-1.38778e-17,0.779375),(-1.38778e-17,-1.38778e-17,0.779375),(-1.38778e-17,0.226,1.30875), + (0.226,0.226,1.30875),(0.226,-1.38778e-17,1.30875),(-1.38778e-17,-1.38778e-17,1.30875),(0.452,0.226,0.25),(0.452,-1.38778e-17,0.25),(0.452,0.226,0.779375), + (0.452,-1.38778e-17,0.779375),(0.452,0.226,1.30875),(0.452,-1.38778e-17,1.30875),(-1.38778e-17,0.452,0.25),(0.226,0.452,0.25),(-1.38778e-17,0.452,0.779375), + (0.226,0.452,0.779375),(-1.38778e-17,0.452,1.30875),(0.226,0.452,1.30875),(0.452,0.452,0.25),(0.452,0.452,0.779375),(0.452,0.452,1.30875),(0.146,0.226,0.779375), + (0.146,-1.38778e-17,0.779375),(0.146,0.226,1.30875),(0.146,-1.38778e-17,1.30875),(0.146,0.452,0.779375),(0.146,0.452,1.30875)]) + c0 = [18, 0, 2, 3, 1, 9, 10, 11, 12, 18, 9, 10, 11, 12, 13, 36, 37, 16, 18, 13, 36, 37, 16, 17, 38, 39, 20, 18, 2, 4, 5, 3, 10, 21, 22, 11, 18, 10, 21, 22, 11, 14, 23, 24, 15, + 18, 14, 23, 24, 15, 18, 25, 26, 19, 18, 6, 7, 2, 0, 27, 28, 10, 9, 18, 27, + 28, 10, 9, 29, 40, 36, 13, 18, 29, 40, 36, 13, 31, 41, 38, 17, 18, 7, 8, 4, 2, 28, 33, 21, 10, 18, 28, 33, 21, 10, 30, 34, 23, 14, 18, 30, 34, 23, 14, 32, 35, 25, 18] + cI0 = [0, 9, 18, 27, 36, 45, 54, 63, 72, 81, 90, 99, 108] + m3 = MEDCouplingUMesh("3D", 3) + m3.setCoords(coo) + m3.setConnectivity(DataArrayInt(c0), DataArrayInt(cI0)) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([7,12,22,27]); grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m3.getNumberOfNodes() + nodesDup, cells1, cells2 = mfu.buildInnerBoundaryAlongM1Group("group") + m3_bis = mfu.getMeshAtLevel(0) + m3_bis.checkConsistency() + m2_bis = mfu.getMeshAtLevel(-1) + m2_bis.checkConsistency() + self.assertEqual(nNod+8, mfu.getNumberOfNodes()) + self.assertEqual(nNod+8, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod+8, m2_bis.getNumberOfNodes()) + self.assertEqual([13, 14, 17, 18, 23, 25, 36, 38], nodesDup.getValues()) + self.assertEqual(m3_bis.getCoords()[nodesDup].getValues(), m3_bis.getCoords()[nNod:].getValues()) + self.assertEqual(set([1, 2, 4, 5]), set(cells1.getValues())) + self.assertEqual(set([7, 8, 10, 11]), set(cells2.getValues())) + self.assertEqual([7, 12, 22, 27],mfu.getGroupArr(-1,"group").getValues()) + self.assertEqual([56, 57, 58, 59],mfu.getGroupArr(-1,"group_dup").getValues()) # here only one cell has been duplicated + m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity() + m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) + pass + @WriteInTmpDir def testBasicConstructors(self): GeneratePyfile18(self) -- 2.39.2