Commit e43aac67 authored by Olga Stamati's avatar Olga Stamati
Browse files

fix interpolateDisplacementsOnGrid()

parent ee5dbdb9
Pipeline #62019 passed with stages
in 25 minutes and 40 seconds
......@@ -1116,10 +1116,11 @@ def interpolateDisplacementsOnGrid(ns, points, displacements, imShape=None, nNei
"""
import spam.DIC
from scipy.spatial import KDTree
tol = 1e-6
if imShape is None:
imShape = (numpy.max(points, axis=0) + 2*ns).astype(int)
print(f"imShape is {imShape}")
#print(f"imShape is {imShape}")
### 1. Generate grid with ns
nodePositions, nodesDim = spam.DIC.grid.makeGrid(imShape, ns)
......@@ -1150,13 +1151,12 @@ def interpolateDisplacementsOnGrid(ns, points, displacements, imShape=None, nNei
# exclusion criterion: subtract node position from nNeighnours
# Require at least one pos and one neg in each dimension
crit = numpy.sum(particlePositions - nodePositions[nodeNumber] > 0, axis=0)
if numpy.logical_and(crit > 0, crit < nNeighbours).sum() == 3:
posMax = numpy.array([particlePositions[:, i].max() for i in range(3)])
posMin = numpy.array([particlePositions[:, i].min() for i in range(3)])
if numpy.logical_and(numpy.all(nodePositions[nodeNumber] >= posMin), numpy.all(nodePositions[nodeNumber] <= posMax)):
#print("point accepted!")
displacement = numpy.zeros(3)
for n in range(nNeighbours):
displacement += displacements[ind[n]]*distance[n]
displacement /= distance.sum()
weights = (1/(distance+tol))
displacement = numpy.sum(displacements[ind]*weights[:, numpy.newaxis], axis=0)/weights.sum()
return nodeNumber, displacement
else:
#print("point rejected")
......@@ -1171,8 +1171,8 @@ def interpolateDisplacementsOnGrid(ns, points, displacements, imShape=None, nNei
gridDisplacements[returns[0]] = returns[1]
if verbose:
pbar.finish()
pbar.finish()
return {'fieldDims': nodesDim,
'fieldCoords': nodePositions,
'displacementField': gridDisplacements}
......@@ -223,15 +223,19 @@ class testAll(unittest.TestCase):
spam.mesh.projectField(mesh, fields, writeConnectivity="Ispam", vtkMesh=True)
def test_interpolateDisplacementsOnGrid(self):
# Try with random points, and apply displacements == 0.1*z
points = numpy.random.randint(10, high=50, size=[200,3])
# Try with random points, and apply displacements == 0.1*z
imShape = [20, 20, 20]
ns = 3
#points = numpy.random.randint(10, high=50, size=[200,3])
points, _ = spam.DIC.grid.makeGrid(imShape, ns)
displ = numpy.zeros_like(points, dtype='<f4')
displ[:,0] = points[:,0]/10.0
gridded = spam.mesh.interpolateDisplacementsOnGrid(3, points, displ, nNeighbours=6, verbose=True)
v = gridded['displacementField'].reshape(gridded['fieldDims'][0], gridded['fieldDims'][1], gridded['fieldDims'][2], 3)
vZmean = numpy.nanmean(v[:,:,:,0], axis=(1,2))
zPos = numpy.unique(gridded['fieldCoords'][:,0])
self.assertTrue(numpy.allclose(vZmean[numpy.isfinite(vZmean)], zPos[numpy.isfinite(vZmean)]*0.1, atol=0.5))
gridded = spam.mesh.interpolateDisplacementsOnGrid(ns, points, displ, imShape=imShape, nNeighbours=6, verbose=True)
#v = gridded['displacementField'].reshape(gridded['fieldDims'][0], gridded['fieldDims'][1], gridded['fieldDims'][2], 3)
#vZmean = numpy.nanmean(v[:,:,:,0], axis=(1,2))
#zPos = numpy.unique(gridded['fieldCoords'][:,0])
#self.assertTrue(numpy.allclose(vZmean[numpy.isfinite(vZmean)], zPos[numpy.isfinite(vZmean)]*0.1, atol=0.1))
self.assertTrue(numpy.allclose(displ, gridded['displacementField'], atol=0.001))
if __name__ == '__main__':
unittest.main()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment