[com5] Fix com5 bug on windows plus inconsistencies in code

avaframe/com1DFA/DFAToolsCython.pyx

@@ -573,8 +573,8 @@ cpdef (double, double, int, int, int, double, double, double, double) samosProje
     # are we in the same cell?
     if Lxi==Lxj and Lyi==Lyj:
       return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3]
-    Lxj = Lxi
-    Lyj = Lyi
+    Lxi = Lxj
+    Lyi = Lyj
   return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3]
 
 
@@ -675,6 +675,9 @@ cpdef (double, double, double, int, int, int, double, double, double, double) di
     # second step: conserve the distance dist (correction ov the position)
     # measure distance between this new point and the previous time step position
     distn = norm(xNew-xPrev, yNew-yPrev, zNew-zPrev)
+    # guard against stopped particle (distn == 0)
+    if distn == 0:
+      break
     # adjust position on the Xprev, Xnew line
     xNew = xPrev + (xNew-xPrev) * dist / distn
     yNew = yPrev + (yNew-yPrev) * dist / distn

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -481,7 +481,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
     if entrMassCell < 0:
       entrMassCell = 0
     entrMassRaster[indCellY, indCellX] = entrMassCell
-  fields['entrMassRaster'] = np.asarray(entrMassRaster)
+  fields['entrMassRaster'] = np.asarray(entrMassRaster).copy()
 
   return particles, force, fields
 
@@ -1115,7 +1115,6 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0):
 
   if stop:
     particles['iterate'] = False
-    particles['stoppedParticles'] = particles
     massStopped = massStopped + massFlowing
 
     if typeStop == 1:
@@ -1398,27 +1397,27 @@ def updateFieldsC(cfg, particles, dem, fields):
       FTStopBilinear[j, i] = - MassStopBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i] 
       FTEntBilinear[j, i] = - MassEntBilinear[j, i] / (areaRaster[j, i] * rhoEnt)
 
-  fields['FM'] = np.asarray(MassBilinear)
-  fields['FV'] = np.asarray(VBilinear)
-  fields['Vx'] = np.asarray(VXBilinear)
-  fields['Vy'] = np.asarray(VYBilinear)
-  fields['Vz'] = np.asarray(VZBilinear)
-  fields['FT'] = np.asarray(FTBilinear)
-  fields['pfv'] = np.asarray(PFV)
-  fields['pft'] = np.asarray(PFT)
-  fields['dmDet'] = np.asarray(DMDet)
-  fields['FTStop'] = np.asarray(FTStopBilinear)
-  fields['FTDet'] = np.asarray(FTDetBilinear)
-  fields['FTEnt'] = np.asarray(FTEntBilinear)
+  fields['FM'] = np.asarray(MassBilinear).copy()
+  fields['FV'] = np.asarray(VBilinear).copy()
+  fields['Vx'] = np.asarray(VXBilinear).copy()
+  fields['Vy'] = np.asarray(VYBilinear).copy()
+  fields['Vz'] = np.asarray(VZBilinear).copy()
+  fields['FT'] = np.asarray(FTBilinear).copy()
+  fields['pfv'] = np.asarray(PFV).copy()
+  fields['pft'] = np.asarray(PFT).copy()
+  fields['dmDet'] = np.asarray(DMDet).copy()
+  fields['FTStop'] = np.asarray(FTStopBilinear).copy()
+  fields['FTDet'] = np.asarray(FTDetBilinear).copy()
+  fields['FTEnt'] = np.asarray(FTEntBilinear).copy()
 
   if computeP:
-    fields['ppr'] = np.asarray(PP)
-    fields['P'] = np.asarray(PBilinear)
+    fields['ppr'] = np.asarray(PP).copy()
+    fields['P'] = np.asarray(PBilinear).copy()
   if computeTA:
-    fields['TA'] = np.asarray(travelAngleField)
-    fields['pta'] = np.asarray(PTA)
+    fields['TA'] = np.asarray(travelAngleField).copy()
+    fields['pta'] = np.asarray(PTA).copy()
   if computeKE:
-    fields['pke'] = np.asarray(PKE)
+    fields['pke'] = np.asarray(PKE).copy()
 
 
   for k in range(nPart):
@@ -1674,7 +1673,9 @@ def computeCohesionForceC(cfg, particles, force):
   force['forceSPHX'] = np.asarray(forceSPHX)
   force['forceSPHY'] = np.asarray(forceSPHY)
   force['forceSPHZ'] = np.asarray(forceSPHZ)
-  particles['bondDist'] = np.asarray(bondDist)
+  # bondDist is modified in-place via the memoryview; do not reassign
+  # particles['bondDist'] here, otherwise numpy will create a VIEW with
+  # base=memoryview which can hang on deallocation on Windows.
   return force, particles
 
 
@@ -1775,9 +1776,9 @@ def initializeBondsC(particles, triangles):
     bondDist[bondStart2[p2+1]-1] = distanceIni
     bondStart2[p2+1] = bondStart2[p2+1] - 1
 
-  particles['bondStart'] = np.asarray(bondStart)
-  particles['bondPart'] = np.asarray(bondPart)
-  particles['bondDist'] = np.asarray(bondDist)
+  particles['bondStart'] = np.asarray(bondStart).copy()
+  particles['bondPart'] = np.asarray(bondPart).copy()
+  particles['bondDist'] = np.asarray(bondDist).copy()
   return particles
 
 def countRemovedBonds(particles, mask, nRemove):
@@ -1870,9 +1871,35 @@ def removedBonds(particles, mask, nRemove, nBondRemove):
          countBond = countBond + 1
 
 
-  particles['bondStart'] = np.asarray(bondStartNew)
-  particles['bondPart'] = np.asarray(bondPartNew)
-  particles['bondDist'] = np.asarray(bondDistNew)
+  # Remap bondStart and bondPart to use compacted (new) particle indices
+  # After removePart compacts particle arrays via np.extract, old indices become stale
+  cdef int[:] oldToNew = np.full(nPart, -1, dtype=np.int32)
+  cdef int newIdx = 0
+  cdef int kk
+  for kk in range(nPart):
+    if keepParticle[kk] == 1:
+      oldToNew[kk] = newIdx
+      newIdx = newIdx + 1
+  cdef int nPartNew = newIdx
+
+  # Rebuild bondStart indexed by new particle indices
+  cdef int[:] bondStartFinal = np.zeros(nPartNew + 1, dtype=np.int32)
+  cdef int bondCumul = 0
+  newIdx = 0
+  for kk in range(nPart):
+    if keepParticle[kk] == 1:
+      bondCumul = bondCumul + (bondStartNew[kk + 1] - bondStartNew[kk])
+      bondStartFinal[newIdx + 1] = bondCumul
+      newIdx = newIdx + 1
+
+  # Remap bondPartNew target indices from old to new
+  cdef int ibb
+  for ibb in range(countBondNew):
+    bondPartNew[ibb] = oldToNew[bondPartNew[ibb]]
+
+  particles['bondStart'] = np.asarray(bondStartFinal).copy()
+  particles['bondPart'] = np.asarray(bondPartNew).copy()
+  particles['bondDist'] = np.asarray(bondDistNew).copy()
   return particles
 
 

avaframe/com1DFA/particleTools.py

@@ -229,7 +229,8 @@ def removePart(particles, mask, nRemove, reasonString="", snowSlide=0):
     for key in particles:
         if key == "nPart":
             particles["nPart"] = particles["nPart"] - nRemove
-        elif key == "stoppedParticles":
+        elif key == "stoppedParticles" or key == "bondStart" or key == "bondPart" or key == "bondDist":
+            # bond arrays are already rebuilt by removedBonds with new-indexed data
             continue
         # for all keys in particles that are arrays of size nPart do:
         elif type(particles[key]).__module__ == np.__name__:

avaframe/tests/test_DFAToolsCython.py

@@ -310,6 +310,105 @@ def test_reprojectionC(capfd):
     assert abs(dist-dist0) <= threshold*(dist0 + csz)
 
 
+def test_samosProjectionConvergence(capfd):
+    """Test samosProjectionIteratrive convergence: result is stable across iteration counts."""
+    ncols = 100
+    nrows = 50
+    csz = 5
+    header = {}
+    header['ncols'] = ncols
+    header['nrows'] = nrows
+    header['cellsize'] = csz
+    dem = {}
+    dem['header'] = header
+    x0 = 300
+    z0 = 500
+    X = np.linspace(0, csz * (ncols - 1), ncols)
+    Y = np.linspace(0, csz * (nrows - 1), nrows)
+    XX, YY = np.meshgrid(X, Y)
+    ZZ = z0 / (x0 * x0) * (XX - x0) * (XX - x0)
+    dem['rasterData'] = ZZ
+    num = 1
+    interpOption = 2
+    dem = gT.getNormalMesh(dem, num=num)
+    Nx = dem['Nx']
+    Ny = dem['Ny']
+    Nz = dem['Nz']
+
+    # Point above the parabola surface, offset to trigger cell change during reprojection
+    xIn = 350
+    yIn = 125
+    Lx0, Ly0, iCell, w0, w1, w2, w3 = DFAfunC.getCellAndWeights(
+        xIn, yIn, ncols, nrows, csz, interpOption)
+    zSurf = DFAfunC.getScalar(Lx0, Ly0, w0, w1, w2, w3, ZZ)
+    zIn = zSurf + 20
+
+    # Reference: many iterations
+    xRef, yRef, iCellRef, LxR0, LyR0, wR0, wR1, wR2, wR3 = DFAfunC.samosProjectionIteratrive(
+        xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, 50)
+    zRef = DFAfunC.getScalar(LxR0, LyR0, wR0, wR1, wR2, wR3, ZZ)
+
+    # With enough iterations, result matches the reference (non-regression)
+    for nIter in [29, 30, 50]:
+        xN, yN, iCellN, LxN0, LyN0, wN0, wN1, wN2, wN3 = DFAfunC.samosProjectionIteratrive(
+            xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, nIter)
+        zN = DFAfunC.getScalar(LxN0, LyN0, wN0, wN1, wN2, wN3, ZZ)
+        assert xN == pytest.approx(xRef, rel=1e-4)
+        assert yN == pytest.approx(yRef, rel=1e-4)
+        assert zN == pytest.approx(zRef, rel=1e-4)
+
+
+def test_distConservStoppedParticle(capfd):
+    """distConservProjectionIteratrive handles distn=0 inside loop without NaN.
+
+    prev is on a flat surface, current is above the surface at the same (x,y).
+    Initial distn is non-zero (current is above surface), so the loop enters.
+    After orthogonal reprojection, the particle lands on the surface and distn
+    becomes zero, triggering the division by zero in the distance conservation step."""
+    ncols = 20
+    nrows = 20
+    csz = 5
+    header = {}
+    header['ncols'] = ncols
+    header['nrows'] = nrows
+    header['cellsize'] = csz
+    dem = {}
+    dem['header'] = header
+    X = np.linspace(0, csz * (ncols - 1), ncols)
+    Y = np.linspace(0, csz * (nrows - 1), nrows)
+    XX, YY = np.meshgrid(X, Y)
+    # Flat surface so normal is straight up, no lateral shift
+    ZZ = np.zeros_like(XX)
+    dem['rasterData'] = ZZ
+    num = 1
+    interpOption = 2
+    dem = gT.getNormalMesh(dem, num=num)
+    Nx = dem['Nx']
+    Ny = dem['Ny']
+    Nz = dem['Nz']
+
+    xPrev = 50.0
+    yPrev = 50.0
+    zPrev = 0.0  # on flat surface
+    xCur = xPrev  # same x, y — no lateral movement
+    yCur = yPrev
+    zCur = 20.0  # above surface
+    threshold = 0.001
+
+    xNew, yNew, zNew, iCell, Lx0, Ly0, w0, w1, w2, w3 = DFAfunC.distConservProjectionIteratrive(
+        xPrev, yPrev, zPrev, ZZ, Nx, Ny, Nz,
+        xCur, yCur, zCur,
+        csz, ncols, nrows, interpOption, 10, threshold)
+
+    # Must not be NaN or inf (would be if division by zero is unguarded)
+    assert not np.isnan(xNew)
+    assert not np.isnan(yNew)
+    assert not np.isnan(zNew)
+    assert np.isfinite(xNew)
+    assert np.isfinite(yNew)
+    assert np.isfinite(zNew)
+
+
 def test_SamosATfric(capfd):
     """ Test the account4FrictionForce function"""
     tau0 = 0