Commit 2fa4c1b3 authored by Edward Andò's avatar Edward Andò
Browse files

geers parallelised

parent 03a90685
......@@ -173,7 +173,7 @@ def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault
return Ffield
def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask=True, bruteForceDistance=False):
def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask=True, bruteForceDistance=False, nProcesses=nProcessesDefault, verbose=True):
"""
This function computes the transformation gradient field F from a given displacement field.
Please note: the transformation gradient tensor: F = I + du/dx.
......@@ -204,6 +204,14 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
or explicitly compute distance to each point (bruteForceDistance=True).
Default = False
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
Ffield: nz x ny x nx x 3x3 array of n cells
......@@ -213,9 +221,7 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
----
Taken from the implementation in "TomoWarp2: A local digital volume correlation code", Tudisco et al., 2017
"""
import numpy
import scipy.spatial
pbar = progressbar.ProgressBar()
##Define dimensions
dims = displacementField.shape[0:3]
......@@ -255,10 +261,15 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
displacementFieldFlatGood = displacementFieldFlat
#build KD-tree for neighbour identification
if not bruteForceDistance:
treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
# Output array for good points
FfieldFlatGood = numpy.zeros_like(FfieldFlat[goodPointsMask])
for goodPoint in pbar(range(fieldCoordsFlatGood.shape[0])):
# Function for parallel mode
global geersOnePoint
def geersOnePoint(goodPoint):
#This is for the linear model, equation 15 in Geers
centralNodePosition = fieldCoordsFlatGood[goodPoint]
centralNodeDisplacement = displacementFieldFlatGood[goodPoint]
......@@ -332,7 +343,22 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
#print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
pass
FfieldFlatGood[goodPoint] = F
return goodPoint, F
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldCoordsFlatGood.shape[0]).start()
finishedPoints = 0
# Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(geersOnePoint, range(fieldCoordsFlatGood.shape[0])):
if verbose:
finishedPoints += 1
pbar.update(finishedPoints)
FfieldFlatGood[returns[0]] = returns[1]
if verbose: pbar.finish()
FfieldFlat[goodPointsMask] = FfieldFlatGood
return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
......
......@@ -370,11 +370,11 @@ class testAll(unittest.TestCase):
displacements = (pointsDef - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
displacements2D = (pointsDef2D - pointsRef2D).reshape(dims2D[0], dims2D[1], dims2D[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, verbose=True)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D, verbose=False)
ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"]
ev2D = spam.deformation.decomposeFfield(Ffield2D, ["vol"])["vol"]
ev = spam.deformation.decomposeFfield(Ffield, ["vol"], verbose=False)["vol"]
ev2D = spam.deformation.decomposeFfield(Ffield2D, ["vol"], verbose=False)["vol"]
self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3)
self.assertAlmostEqual(ev2D.mean(), 2. * (dilation - 1), places=3)
......@@ -386,7 +386,7 @@ class testAll(unittest.TestCase):
# Deform the points by spreading them apart -- this is dilating so positive volumetric strain
displacements = (pointsRef * [zStretch, 1.0, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, verbose=False)
# only compute on internal nodes, we know there's a problem at the edges
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), zStretch, places=3)
......@@ -401,8 +401,8 @@ class testAll(unittest.TestCase):
# Deform the points by spreading them apart -- this is dilating so positive volumetric strain
displacements = (pointsRef * [1.0, yStretch, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
displacements2D = (pointsRef2D * [1.0, yStretch, 1.0] - pointsRef2D).reshape(dims2D[0], dims2D[1], dims2D[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, verbose=False)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D, verbose=False)
# only compute on internal nodes, we know there's a problem at the edges
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), 1.0, places=3)
......@@ -433,12 +433,12 @@ class testAll(unittest.TestCase):
displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
displacements2D = (pointsRotated2D - pointsRef2D).reshape(dims2D[0], dims2D[1], dims2D[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, verbose=False)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D, verbose=False)
# Compute strains
decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"])
decomposedF2D = spam.deformation.decomposeFfield(Ffield2D, ["r", "vol", "dev"])
decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"], verbose=False)
decomposedF2D = spam.deformation.decomposeFfield(Ffield2D, ["r", "vol", "dev"], verbose=False)
self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 0].mean(), rotAngle, places=3)
self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 1].mean(), 0.0, places=3)
......@@ -555,6 +555,7 @@ class testAll(unittest.TestCase):
self.assertEqual(numpy.isnan(Ffield[:,:,:,2,2]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield2D[:,:,:,2,2]).sum(), 1)
def test_bagiStrain(self):
#######################################################
# Check triangulation and strains
......
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