Commit ee5dbdb9 authored by Edward Andò's avatar Edward Andò
Browse files

interpolateDisplacementsOnGrid() function in unstructured mesh + test

parent 828cd14f
Pipeline #61648 passed with stages
in 25 minutes and 49 seconds
......@@ -34,6 +34,10 @@ WARNING
"""
import numpy
import spam.mesh
import multiprocessing
# Global number of processes
nProcessesDefault = multiprocessing.cpu_count()
def createCuboid(lengths, lcar, origin=[0., 0., 0.], verbose=False, gmshFile=None, vtkFile=None):
......@@ -689,8 +693,8 @@ def tetVolumes(points, connectivity):
return volumes
def create2Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gmshFile=None, vtkFile=None):
def create2Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gmshFile=None, vtkFile=None):
import pygmsh
import meshio
......@@ -876,8 +880,9 @@ def create2Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gms
return
def create8Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gmshFile=None, vtkFile=None):
# NOT USED
# NOT USED
import pygmsh
import meshio
......@@ -1061,3 +1066,113 @@ def create8Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gms
return
def interpolateDisplacementsOnGrid(ns, points, displacements, imShape=None, nNeighbours=6, nProcesses=nProcessesDefault, verbose=False):
"""
This function takes displacements defined on a cloud of points
(typically particle centres from a labelled image and displacements
measured from ddic) and interpolates the displacements onto a 3D
grid defined by "ns".
Typical use case is to use the regular strain calcuation tools on
a discrete field.
This is similar to the approach used in (which uses triangular
shape functions to interpolate displacements):
Desrues, J., & Viggiani, G. (2004). Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry. International Journal for Numerical and Analytical Methods in Geomechanics, 28(4), 279-321.
Parameters
----------
ns : int
Node spacing for computation of regular grid
points : 3xN numpy array of floats
Points at which displacements are defined
displacements : 3xN numpy array of floats
Displacements defined at points
imShape : 3-component list of ints, optional
Size of domain to grid,
if not set, the grid will be the max of the point cloud + 2*ns
nNeighbours : int, optional
Number of neighbours to use for interpolation, should be at least 4.
Default = 6
nProcesses : int, optional
Number of parallel processes to use.
Default = number of CPUs detected by multiprocessing.cpu_count()
verbose : bool, optional
Show progress bar?
Default = False
Returns
-------
Dictionary containing:
- "fieldDims": fieldDims
- "fieldCoords": fieldCoords
- "displacementField": displacementField
"""
import spam.DIC
from scipy.spatial import KDTree
if imShape is None:
imShape = (numpy.max(points, axis=0) + 2*ns).astype(int)
print(f"imShape is {imShape}")
### 1. Generate grid with ns
nodePositions, nodesDim = spam.DIC.grid.makeGrid(imShape, ns)
numberOfNodes = nodesDim[0]*nodesDim[1]*nodesDim[2]
### 0. Create displacement field output array
gridDisplacements = numpy.zeros((numberOfNodes, 3), dtype='<f4')
### 2. Prepare KD-tree with displacements and points
treeCoord = KDTree(points)
### 3. for each grid point:
### see if inside or outside of data cloud
### if inside, interpolate displacement, fill in output array
if verbose:
import progressbar
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
pbar.start()
finishedNodes = 0
global interpolateOneNode
def interpolateOneNode(nodeNumber):
#print(f"looking up point {nodeNumber} at coord {nodePositions[nodeNumber]}")
distance, ind = treeCoord.query(nodePositions[nodeNumber], k=nNeighbours)
particlePositions = points[ind]
# 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:
#print("point accepted!")
displacement = numpy.zeros(3)
for n in range(nNeighbours):
displacement += displacements[ind[n]]*distance[n]
displacement /= distance.sum()
return nodeNumber, displacement
else:
#print("point rejected")
return nodeNumber, numpy.array([numpy.nan, numpy.nan, numpy.nan])
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(interpolateOneNode, range(numberOfNodes)):
if verbose:
finishedNodes += 1
pbar.update(finishedNodes)
gridDisplacements[returns[0]] = returns[1]
if verbose:
pbar.finish()
return {'fieldDims': nodesDim,
'fieldCoords': nodePositions,
'displacementField': gridDisplacements}
......@@ -222,6 +222,16 @@ class testAll(unittest.TestCase):
fields = {"origin": [0, 0, 0], "lengths": [10, 10, 10], "nCells": [5, 5, 5], "values": values}
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])
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))
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