From 4bf4c782826685ea09801c22c4ff8a39124bbb2b Mon Sep 17 00:00:00 2001 From: abn Date: Fri, 8 Oct 2021 15:52:07 +0200 Subject: [PATCH] Python version of some of the core algorithms of buildInnerBoundaryAlongM1Group() --- resources/dev/my_findNodesToDup.py | 197 +++++++++++++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 resources/dev/my_findNodesToDup.py diff --git a/resources/dev/my_findNodesToDup.py b/resources/dev/my_findNodesToDup.py new file mode 100644 index 000000000..6ba2b41ea --- /dev/null +++ b/resources/dev/my_findNodesToDup.py @@ -0,0 +1,197 @@ +""" Python version of findNodesToDuplicate() and findCellsToRenumber() which are at the core of the + buildInnerBoundaryAlongM1Group() algorithm. + This greatly helps algorithm tuning ... +""" + +from medcoupling import * + +def findNodesToDuplicate(this, otherDimM1OnSameCoords): + # Checking star-shaped M1 group: + meshM2, _,_,_,rdit0 = otherDimM1OnSameCoords.buildDescendingConnectivity() # 2D: a mesh of points, 3D: a mesh of segs + dsi = rdit0.deltaShiftIndex() + 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()): + raise ValueError("") + + # 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 + xtremIdsM2 = dsi.findIdsEqual(1) + meshM2Part = meshM2[xtremIdsM2] + xtrem = meshM2Part.computeFetchedNodeIds() + # Remove from the list points on the boundary of the M0 mesh (those need duplication!) + m0desc, dt0, dit0, rdt0, rdit0 = this.buildDescendingConnectivity() + dsi = rdit0.deltaShiftIndex() + boundSegs = dsi.findIdsEqual(1) # boundary segs/faces of the M0 mesh + m0descSkin = m0desc[boundSegs] + fNodes = m0descSkin.computeFetchedNodeIds() + # 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. + if this.getMeshDimension() == 3 : + m0descSkinDesc, _, _, _, _ = m0descSkin.buildDescendingConnectivity() # all segments of the skin of the 3D (M0) mesh + _, corresp = meshM2.areCellsIncludedIn(m0descSkinDesc,2) + # validIds is the list of segments which are on both the skin of *this*, and in the segments of the M1 group + # In the cube example above, this is a U shape polyline. + validIds = corresp.findIdsInRange(0, meshM2.getNumberOfCells()) + 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: + # (the U-shaped polyline described above) + m1IntersecSkin = m0descSkinDesc[validIds] + # Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above) + notDuplSkin = m1IntersecSkin.findBoundaryNodes() + fNodes1 = fNodes.buildSubstraction(notDuplSkin) + + # Specific logic to handle singular points : + # - take the sub-mesh of M1 (dim N-1) built keeping only elements touching the m1IntersecSkin + # - build its spread zones + # - a node common to two different spread zones actually connects two cells of M1 by just one point : this is a singular point + # which should not be duplicated. + # This is a costly algorithm so only go into it if a simple (non sufficient) criteria is met: a node connected to more than 3 segs in meshM2: + meshM2Desc, _, _, _, rdit0 = meshM2.buildDescendingConnectivity() # a mesh made of node cells + dsi = rdit0.deltaShiftIndex() + singPoints = dsi.findIdsNotInRange(-1,4) # points connected to (strictly) more than 3 segments + if singPoints.getNumberOfTuples(): + print ("Hitting singular point logic") + boundNodes = m1IntersecSkin.computeFetchedNodeIds() + boundCells = otherDimM1OnSameCoords.getCellIdsLyingOnNodes(boundNodes, False) # false= take cell in, even if not all nodes are in dupl + m1Bound = otherDimM1OnSameCoords[boundCells] + m1Bound.writeVTK("/tmp/m1Bound.vtu") + _, d,dI,revD,revDI = m1Bound.buildDescendingConnectivity() + neighSing, neighISing = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(d,dI,revD,revDI) + seed = 0 + hitCells = DataArrayInt.New(); hitCells.alloc(m1Bound.getNumberOfCells()); hitCells.fillWithValue(-1) + nIter, nIterMax = 0, m1Bound.getNumberOfCells()+1 + while nIter < nIterMax: + sprdZone, _ = MEDCouplingUMesh.ComputeSpreadZoneGraduallyFromSeed([seed], neighSing,neighISing, -1) + hitCells[sprdZone] = nIter + nIter += 1 + noHit = hitCells.findIdsEqual(-1) + if noHit.getNumberOfTuples(): + seed = noHit[0,0] + else: + break + if nIter >= nIterMax: + raise ValueError("Internal error should not happen") + print("dbg", hitCells.getValues()) + if hitCells.buildUniqueNotSorted().getNumberOfTuples() > 1: + # All points common to two (or more) different spread zones are singular: + cc, ccI = m1Bound.getReverseNodalConnectivity() + spreadCC = hitCells[cc] + singPoints = DataArrayInt.New(); singPoints.alloc(0,1) + for j in range(ccI.getNumberOfTuples()-1): + ext = spreadCC[ccI[j]:ccI[j+1]] + if ext.getNumberOfTuples() and not ext.isUniform(ext[0,0]): + singPoints.pushBackSilent(j) + fNodes1 = fNodes1.buildSubstraction(singPoints) + print ("sing,", singPoints.getValues()) + notDup = xtrem.buildSubstraction(fNodes1) + else: # if validIds ... + notDup = xtrem.buildSubstraction(fNodes) + else: # if 3D ... + notDup = xtrem.buildSubstraction(fNodes) + + # Now compute cells around group (i.e. cells where we will do the propagation to identify the two sub-sets delimited by the group) + m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds() + dupl = m1Nodes.buildSubstraction(notDup) + return dupl + + +def findCellsToRenumber(this, otherDimM1OnSameCoords, dupl): + """ Find cells to renumber + """ + cellsAroundGroup = this.getCellIdsLyingOnNodes(dupl, False) # false= take cell in, even if not all nodes are in dupl + + # + mAroundGrp=this[cellsAroundGroup] + nCells2 = mAroundGrp.getNumberOfCells() + mArGrpDesc,desc00,descI00,revDesc00,revDescI00=mAroundGrp.buildDescendingConnectivity() + + mAroundGrp.writeVTK("/tmp/mAr.vtu") + mArGrpDesc.writeVTK("/tmp/mAr_desc.vtu") + + # Identify cells of M1 group in sub-mesh mAroundGrp + _, idsOfM1 = mArGrpDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2) + nCellsArGrpDesc = mArGrpDesc.getNumberOfCells() +# print(idsOfM1.getValues()) + + # Build map giving for each cell ID in mAroundGrp the corresponding cell ID on the other side of the crack: + # Note that this does not cover all cells around the crack (a cell like a triangle might touche the crack only with its tip) + toOtherSide = {} + for i, v in enumerate(idsOfM1): + if v[0] >= nCellsArGrpDesc: # Keep valid match only + continue + idx0 = revDescI00[v[0], 0] + c1, c2 = revDesc00[idx0, 0], revDesc00[idx0+1,0] + toOtherSide[c1] = c2 + toOtherSide[c2] = c1 + + # Neighbor information of the mesh WITH the crack (some neighbors are removed): + # 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): + DataArrayInt.RemoveIdsFromIndexedArrays(idsOfM1,desc00,descI00) + # Compute the neighbor of each cell in mAroundGrp, 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. + neigh00, neighI00 = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00) + + # For each initial connex part of the M1 mesh (or said differently for each independent crack): + seed, nIter = 0, 0 + nIterMax = nCells2+1 # Safety net for the loop + hitCells = DataArrayInt.New(); hitCells.alloc(nCells2) + hitCells.fillWithValue(-1) # -1 : not hit, 0: one side of the crack, 1: other side of the crack + cellsToModifyConn0_torenum = DataArrayInt.New() + cellsToModifyConn0_torenum.alloc(0,1) + while nIter < nIterMax: +# print("dbg ", hitCells.getValues()) + t = hitCells.findIdsEqual(-1) + if not t.getNumberOfTuples(): + break + seed = t[0,0] + done = False + PING_FULL, PING_PART = 1,11 + PONG_FULL, PONG_PART = 2,22 + while not done and nIter < nIterMax: # Start of the ping-pong + nIter += 1 + # Identify connex zone around the seed + spreadZone, _ = MEDCouplingUMesh.ComputeSpreadZoneGraduallyFromSeed([seed], neigh00,neighI00, -1) + done = True + for i, s in enumerate(spreadZone.getValues()): + hitCells[s] = PING_FULL + if s in toOtherSide: + other = toOtherSide[s] + if hitCells[other] != PONG_FULL: + done = False + hitCells[other] = PONG_PART + # Compute next seed, i.e. a cell on the other side of the crack + seed = other + if done: + # we might have several disjoing PONG parts in front of a single PING connex part: + idsPong = hitCells.findIdsEqual(PONG_PART) + if idsPong.getNumberOfTuples(): + seed = idsPong[0,0] + done = False + continue # continue without switching side (or break if done remains false) + else: + # Go the other side + PING_FULL, PONG_FULL = PONG_FULL, PING_FULL + PING_PART, PONG_PART = PONG_PART, PING_PART + + nonHitCells = hitCells.findIdsEqual(-1) + if nonHitCells.getNumberOfTuples(): + seed = nonHitCells[0,0] + else: + break + if nIter >= nIterMax: + raise ValueError("Too many iterations - should not happen") + + cellsToModifyConn0_torenum = hitCells.findIdsEqual(PONG_FULL) + cellsToModifyConn1_torenum = hitCells.findIdsEqual(PING_FULL) + if cellsToModifyConn0_torenum.getNumberOfTuples() + cellsToModifyConn1_torenum.getNumberOfTuples() != cellsAroundGroup.getNumberOfTuples(): + raise ValueError("Some cells not hit - Internal error should not happen") + cellsToModifyConn0_torenum.transformWithIndArr(cellsAroundGroup) + cellsToModifyConn1_torenum.transformWithIndArr(cellsAroundGroup) + # + cellIdsNeededToBeRenum=cellsToModifyConn0_torenum + cellIdsNotModified=cellsToModifyConn1_torenum + + return cellIdsNeededToBeRenum, cellIdsNotModified + -- 2.39.2