--- /dev/null
+""" 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
+
