Skip to content

GitLab

Commit 65710e44 authored by Olga Stamati's avatar Olga Stamati
Browse files

implement Geers transformation gradient in 2D (closes #191)

parent 894bf309
Pipeline #53366 passed with stages
in 24 minutes and 17 seconds
Hide whitespace changes
Inline Side-by-side
scripts/spam-regularStrain
View file @ 65710e44
......@@ -64,9 +64,6 @@ if dims[0] == 1:
twoD = True
#spaceZ = 0
print("spam-regularStrain: detected 2D field")
if not args.Q8:
print("\tIMPORTANT: overriding the strain calculation mode to Q8")
args.Q8 = True
else:
twoD = False
#spaceZ = fieldCoords[dims[2]*dims[1],0] - fieldCoords[0,0]
......@@ -131,7 +128,6 @@ else:
# Save strain fields
print("\nspam-regularStrain: Saving strain fields...")
if args.TSV:
# Positions for Geers are the measurement points
if not args.Q8:
if twoD:
......
tools/deformation/deformationField.py
View file @ 65710e44
......@@ -177,19 +177,27 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
import scipy.spatial
pbar = progressbar.ProgressBar()
#Define dimensions
##Define dimensions
dims = displacementField.shape[0:3]
nNodes = dims[0]*dims[1]*dims[2]
displacementFieldFlat = displacementField.reshape(nNodes, 3)
# Check if a 2D field is passed
twoD = False
if dims[0] == 1: twoD = True
#Deformation gradient tensor F = du/dx +I
Ffield = numpy.zeros((dims[0], dims[1], dims[2], 3, 3))
FfieldFlat = numpy.zeros((nNodes, 3, 3))
#Define field coordinates for KD-tree distance calculation, warning, the top point is at zero
fieldCoordsFlat = numpy.mgrid[nodeSpacing[0]:dims[0]*nodeSpacing[0]+nodeSpacing[0]:nodeSpacing[0],
nodeSpacing[1]:dims[1]*nodeSpacing[1]+nodeSpacing[1]:nodeSpacing[1],
nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
if twoD:
fieldCoordsFlat = numpy.mgrid[0:1:1,
nodeSpacing[1]:dims[1]*nodeSpacing[1]+nodeSpacing[1]:nodeSpacing[1],
nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
else:
fieldCoordsFlat = numpy.mgrid[nodeSpacing[0]:dims[0]*nodeSpacing[0]+nodeSpacing[0]:nodeSpacing[0],
nodeSpacing[1]:dims[1]*nodeSpacing[1]+nodeSpacing[1]:nodeSpacing[1],
nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
#Get non-nan displacements
if mask:
......@@ -214,8 +222,8 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
#This is for the linear model, equation 15 in Geers
centralNodePosition = fieldCoordsFlatGood[goodPoint]
centralNodeDisplacement = displacementFieldFlatGood[goodPoint]
sX0X0 = numpy.zeros((3,3))
sX0Xt = numpy.zeros((3,3))
sX0X0 = numpy.zeros((3, 3))
sX0Xt = numpy.zeros((3, 3))
m0 = numpy.zeros((3))
mt = numpy.zeros((3))
......@@ -266,13 +274,25 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
A = sX0X0 - numpy.dot(m0, m0)
C = sX0Xt - numpy.dot(m0, mt)
F = numpy.zeros((3, 3))
if twoD:
A = A[1:, 1:]
C = C[1:, 1:]
try:
F[1:, 1:] = numpy.dot(numpy.linalg.inv(A), C)
F[0, 0] = 1.0
except numpy.linalg.linalg.LinAlgError:
#print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
pass
else:
try:
F = numpy.dot(numpy.linalg.inv(A), C)
except numpy.linalg.linalg.LinAlgError:
#print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
pass
try:
FfieldFlatGood[goodPoint] = numpy.dot(numpy.linalg.inv(A), C)
except numpy.linalg.linalg.LinAlgError:
pass
#print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
FfieldFlatGood[goodPoint] = F
FfieldFlat[goodPointsMask] = FfieldFlatGood
return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
......
tools/tests/test_deformation.py
View file @ 65710e44
......@@ -342,65 +342,74 @@ class testAll(unittest.TestCase):
def test_geers(self):
#nodeSpacing = [150,150,150]
nodeSpacing = [200,200,200]
imSize = [2000,1000,1000]
nodeSpacing = [200, 200, 200]
imSize = [2000, 1000, 1000]
# 2D
nodeSpacing2D = [1, 100, 100]
imSize2D = [1, 1000, 1000]
# Random points
pointsRef, dims = spam.DIC.grid.makeGrid(imSize, nodeSpacing=nodeSpacing)
pointsRef, dims = spam.DIC.makeGrid(imSize, nodeSpacing=nodeSpacing)
pointsRef2D, dims2D = spam.DIC.makeGrid(imSize2D, nodeSpacing=nodeSpacing2D)
#######################################################
########################################################
# Dilation of points for grain strain calc
#######################################################
# Test with 1% dilation
dilation = 1.01
# make sure average volumetric strain is the same as the projected one
pointsDef = pointsRef * dilation
displacements = (pointsDef - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
pointsDef = pointsRef * dilation
pointsDef2D = pointsRef2D * dilation
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)
ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"]
# print("\tU trace mean new:",numpy.array(U)[:,[0,1,2],[0,1,2]].mean())
# print("\tev",numpy.array(ev))
# print("\tev",numpy.array(ev).mean())
ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"]
ev2D = spam.deformation.decomposeFfield(Ffield2D, ["vol"])["vol"]
self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3)
#self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3)
self.assertAlmostEqual(ev2D.mean(), 2. * (dilation - 1), places=3)
########################################################
## Z-stretch of points for strain calc
########################################################
#########################################################
### Z-stretch of points for strain calc
#########################################################
zStretch = 1.02
print("\n\nSpreading points in Z by by 1.2")
# 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)
displacements = (pointsRef * [zStretch, 1.0, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
# 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)
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), 1.0, places=3)
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3)
########################################################
## Y-stretch of points for strain calc
########################################################
yStretch = 1.02
print("\n\nSpreading points in Y by by 1.2")
# 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)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
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)
# 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)
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), yStretch, places=3)
self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3)
self.assertAlmostEqual(Ffield2D[:,1:-1,1:-1,1,1].mean(), yStretch, places=3)
self.assertAlmostEqual(Ffield2D[:,1:-1,1:-1,2,2].mean(), 1.0, places=3)
#######################################################
# Rotation of points for grain strain calc
#######################################################
pointsRotated = pointsRef.copy()
########################################################
## Rotation of points for grain strain calc
########################################################
pointsRotated = pointsRef.copy()
pointsRotated2D = pointsRef2D.copy()
rotAngle = 3.0
for n, point in enumerate(pointsRef):
Phi = spam.deformation.computePhi({'r': [rotAngle, 0.0, 0.0]},
......@@ -409,49 +418,78 @@ class testAll(unittest.TestCase):
pointsRotated[n] += Phi[0:3, -1]
# print( point, pointsRotated[n], '\n\n' )
displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
for n, point in enumerate(pointsRef2D):
Phi = spam.deformation.computePhi({'r': [rotAngle, 0.0, 0.0]},
PhiCentre = numpy.array(((numpy.array(imSize2D)-1)/2)-1),
PhiPoint = point)
pointsRotated2D[n] += Phi[0:3, -1]
# print( point, pointsRotated[n], '\n\n' )
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)
# Compute strains
decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"])
decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"])
decomposedF2D = spam.deformation.decomposeFfield(Ffield2D, ["r", "vol", "dev"])
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)
self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 2].mean(), 0.0, places=3)
self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 0, places=3)
self.assertAlmostEqual(decomposedF['dev'][1:-1,1:-1,1:-1].mean(), 0, places=3)
self.assertAlmostEqual(decomposedF2D['r'][:,1:-1,1:-1, 0].mean(), rotAngle, places=3)
self.assertAlmostEqual(decomposedF2D['r'][:,1:-1,1:-1, 1].mean(), 0.0, places=3)
self.assertAlmostEqual(decomposedF2D['r'][:,1:-1,1:-1, 2].mean(), 0.0, places=3)
self.assertAlmostEqual(decomposedF2D['vol'][:,1:-1,1:-1].mean(), 0, places=3)
self.assertAlmostEqual(decomposedF2D['dev'][:,1:-1,1:-1].mean(), 0, places=3)
#######################################################
# Shear of points for grain strain calc:
# N.B.: this is a shear which results in a SYMMETRIC F
#######################################################
pointsRotated = pointsRef.copy()
########################################################
## Shear of points for grain strain calc:
## N.B.: this is a shear which results in a SYMMETRIC F
########################################################
pointsSheared = pointsRef.copy()
shearVal = 0.1
for n, point in enumerate(pointsRef):
Phi = spam.deformation.computePhi({'s': [shearVal, 0.0, 0.0]},
PhiCentre = numpy.array(((numpy.array(imSize)-1)/2)-1),
PhiPoint = point)
pointsRotated[n] += Phi[0:3, -1]
# print( point, pointsRotated[n], '\n\n' )
pointsSheared[n] += Phi[0:3, -1]
displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
pointsSheared2D = pointsRef2D.copy()
for n, point in enumerate(pointsRef2D):
Phi = spam.deformation.computePhi({'s': [0.0, 0.0, shearVal]},
PhiCentre = numpy.array(((numpy.array(imSize2D)-1)/2)-1),
PhiPoint = point)
pointsSheared2D[n] += Phi[0:3, -1]
displacements = (pointsSheared - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
displacements2D = (pointsSheared2D - pointsRef2D).reshape(dims2D[0], dims2D[1], dims2D[2], 3)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
# Compute strains
decomposedF = spam.deformation.decomposeFfield(Ffield, ["U"])
decomposedF = spam.deformation.decomposeFfield(Ffield, ["U"])
decomposedF2D = spam.deformation.decomposeFfield(Ffield2D, ["U"])
UmeanWithoutBoundaries = numpy.mean(decomposedF['U'][1:-1,1:-1,1:-1], axis=(0, 1, 2))
UmeanWithoutBoundaries2D = numpy.mean(decomposedF2D['U'][0,1:-1,1:-1], axis=(0, 1))
UmeanWithoutBoundaries = numpy.mean(decomposedF['U'][1:-1,1:-1,1:-1], axis=(0, 1, 2))
for i in range(3):
for j in range(i,3):
self.assertAlmostEqual(UmeanWithoutBoundaries[i,j], UmeanWithoutBoundaries[j,i], places=3)
self.assertAlmostEqual(UmeanWithoutBoundaries[0,1], shearVal, places=3)
for i in range(3):
for j in range(i,3):
self.assertAlmostEqual(UmeanWithoutBoundaries2D[i,j], UmeanWithoutBoundaries2D[j,i], places=3)
self.assertAlmostEqual(UmeanWithoutBoundaries2D[2,1], shearVal, places=3)
#######################################################
# Homogeneous shrinking together
#######################################################
########################################################
## Homogeneous shrinking together
########################################################
shrink = 0.99
# Deform the points by pushing them closer together -- this is compressing so negative volumetric strain
displacements = (pointsRef * shrink - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
......@@ -464,11 +502,9 @@ class testAll(unittest.TestCase):
decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"])
self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3)
#######################################################
# Increasing neighbourhood radius should decrease noise
#######################################################
########################################################
## Increasing neighbourhood radius should decrease noise
########################################################
#Add noise to displacements
displacementAmplitude = displacements.max()-displacements.min()
displacements += numpy.random.random(displacements.shape)*0.1*displacementAmplitude
......@@ -490,25 +526,28 @@ class testAll(unittest.TestCase):
self.assertEqual(std1 > std1p5, True)
self.assertEqual(std1p5 > std2, True)
#######################################################
# Application of mask
#######################################################
zStretch = 1.02
print("\n\nSpreading points in Z by by 1.2")
########################################################
## Application of mask
########################################################
xStretch = 1.02
# 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)
displacements = (pointsRef * [1.0, 1.0, xStretch] - pointsRef).reshape(dims[0], dims[1], dims[2], 3)
displacements2D = (pointsRef2D * [1.0, 1.0, xStretch] - pointsRef2D).reshape(dims2D[0], dims2D[1], dims2D[2], 3)
# manually add a NaN in the middle of the displacement field
displacements[dims[0]//2, dims[1]//2, dims[2]//2] = [numpy.nan, numpy.nan, numpy.nan]
displacements2D[0, dims2D[1]//2, dims2D[2]//2] = [numpy.nan, numpy.nan, numpy.nan]
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
# only compute on internal nodes, we know there's a problem at the edges
self.assertEqual(numpy.isnan(Ffield[:,:,:,0,0]).sum(), 1)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, mask=True)
self.assertEqual(numpy.isnan(Ffield[:,:,:,0,0]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield[:,:,:,2,2]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield2D[:,:,:,2,2]).sum(), 1)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, mask=True)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D, mask=True)
self.assertEqual(numpy.isnan(Ffield[:,:,:,2,2]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield2D[:,:,:,2,2]).sum(), 1)
def test_bagiStrain(self):
#######################################################
......@@ -524,7 +563,6 @@ class testAll(unittest.TestCase):
# [ 0, 100, 100],
# [100, 100, 100]], dtype='<f4')
#######################################################
# Check connectivityFromVoronoi
#######################################################
......
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