Skip to content

GitLab

Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit 5a73e117 authored by Edward Andò's avatar Edward Andò
Browse files

Closes #190, spam-regularStrains spam-discreteStrains and the... 

Closes #190, spam-regularStrains spam-discreteStrains and the spam.deformation.Ffield... and spam.deformation.decompose...Field families now multiprocessed, bye bye coverage for now
parent 2fa4c1b3
Pipeline #63346 passed with stages
in 25 minutes and 47 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 244 additions and 183 deletions
scripts/spam-discreteStrain
View file @ 5a73e117
......@@ -30,10 +30,13 @@ either projected back to grains (whereby the value at each grain is a weighted l
...or projected onto a regular grid.
"""
from __future__ import print_function
import numpy
import os
os.environ['OPENBLAS_NUM_THREADS'] = '1'
import argparse
import numpy
import multiprocessing
import spam.helpers
import spam.DIC
import spam.deformation
......@@ -66,6 +69,9 @@ else: start = 1
# if args.SMALL_STRAINS:
# largeStrains = False
# Figure out processes if not passed
if args.PROCESSES is None: args.PROCESSES = multiprocessing.cpu_count()
print("\nCurrent Settings:")
argsDict = vars(args)
for key in sorted(argsDict):
......@@ -179,7 +185,7 @@ if not triangulationAvailable:
#connectivity = connectivity[goodTets]
print("spam-discreteStrain: Computing F=I+du/dx for all tetrahedra")
Ffield = spam.deformation.FfieldBagi(points, connectivity, displacements)
Ffield = spam.deformation.FfieldBagi(points, connectivity, displacements, verbose=True, nProcesses=args.PROCESSES)
#strainMatrix, F, R, volStrain, devStrain = spam.deformation.FfieldBagi(points, connectivity, displacements, onlyStrain=False)
# Compute bagi strains with initial and deformed centres of mass.
......@@ -189,7 +195,7 @@ if args.PROJECT_TO_GRAINS:
# We need to project F, since it is in ZYX, U is in the eigendirections and cannot be summed.
Fgrains = spam.mesh.projectTetFieldToGrains(points+displacements, connectivity, Ffield)
decomposedFfield = spam.deformation.decomposeFfield(Fgrains, args.COMPONENTS, twoD=False)
decomposedFfield = spam.deformation.decomposeFfield(Fgrains, args.COMPONENTS, twoD=False, verbose=True, nProcesses=args.PROCESSES)
spam.helpers.writeStrainTSV(args.OUT_DIR+"/"+args.PREFIX+"-grainProjection.tsv",
points, decomposedFfield, firstColumn="Label", startRow=start)
......@@ -197,7 +203,7 @@ if args.PROJECT_TO_GRAINS:
print("\nspam-discreteStrain: Decomposing F into ", args.COMPONENTS, "for all tetrahedra")
decomposedFfield = spam.deformation.decomposeFfield(Ffield, args.COMPONENTS, twoD=False)
decomposedFfield = spam.deformation.decomposeFfield(Ffield, args.COMPONENTS, twoD=False, verbose=True, nProcesses=args.PROCESSES)
#if args.VTK:
......
scripts/spam-regularStrain
View file @ 5a73e117
......@@ -18,16 +18,20 @@ You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
"""
from __future__ import print_function
import os
os.environ['OPENBLAS_NUM_THREADS'] = '1'
import argparse
import tifffile
import numpy
import multiprocessing
import argparse
import spam.helpers
import spam.DIC
import spam.deformation
import spam.mesh
# Define argument parser object
parser = argparse.ArgumentParser(description="spam-regularStrain "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
"This script computes different components of strain, given a regularly-spaced displacement"+\
......@@ -37,6 +41,9 @@ parser = argparse.ArgumentParser(description="spam-regularStrain "+spam.helpers.
# Parse arguments with external helper function
args = spam.helpers.optionsParser.regularStrainParser(parser)
# Figure out processes if not passed
if args.PROCESSES is None: args.PROCESSES = multiprocessing.cpu_count()
print("+----------------------------+")
print("| Regular Strain Calculation |")
print("+----------------------------+")
......@@ -54,19 +61,15 @@ dims = f["fieldDims"]
fieldCoords = f["fieldCoords"]
# Calculate node spacing for each direction
#spaceX = fieldCoords[1,2] - fieldCoords[0,2]
#spaceY = fieldCoords[dims[2],1] - fieldCoords[0,1]
# 2020-08-31 OS: safer calculation of node spacing
nodeSpacing = numpy.array([numpy.unique(fieldCoords[:, i])[1] - numpy.unique(fieldCoords[:, i])[0] if len(numpy.unique(fieldCoords[:, i])) > 1 else numpy.unique(fieldCoords[:, i])[0] for i in range(3)])
# Catch 2D case
if dims[0] == 1:
twoD = True
#spaceZ = 0
print("spam-regularStrain: detected 2D field")
else:
twoD = False
#spaceZ = fieldCoords[dims[2]*dims[1],0] - fieldCoords[0,0]
# Check if a mask of points (based on return status of the correlation) is asked
......@@ -103,15 +106,25 @@ else:
print("\nspam-regularStrain: Computing F=I+du/dx")
if not args.Q8:
Ffield = spam.deformation.FfieldRegularGeers(disp, nodeSpacing=nodeSpacing,
Ffield = spam.deformation.FfieldRegularGeers(disp,
nodeSpacing=nodeSpacing,
mask=args.MASK,
neighbourRadius=args.STRAIN_NEIGHBOUR_RADIUS)
neighbourRadius=args.STRAIN_NEIGHBOUR_RADIUS,
nProcesses=args.PROCESSES,
verbose=True)
else:
Ffield = spam.deformation.FfieldRegularQ8( disp, nodeSpacing=nodeSpacing)
Ffield = spam.deformation.FfieldRegularQ8( disp,
nodeSpacing=nodeSpacing,
nProcesses=args.PROCESSES,
verbose=True)
# Now compute what's been asked for...
print("\nspam-regularStrain: Decomposing F into ", args.COMPONENTS)
decomposedFfield = spam.deformation.decomposeFfield(Ffield, args.COMPONENTS, twoD=twoD)
decomposedFfield = spam.deformation.decomposeFfield(Ffield,
args.COMPONENTS,
twoD=twoD,
nProcesses=args.PROCESSES,
verbose=True)
# Define base fileName
......@@ -172,7 +185,7 @@ if args.VTK:
else: aspectRatio = [ 1, nodeSpacing[1], nodeSpacing[2]]
# For geers strains are at the measurement points
# Ad per the displacements coming out of spam-ldic this will plot nicely if 2xHWS = NS
# As per the displacements coming out of spam-ldic this will plot nicely if 2xHWS = NS
if not args.Q8:
origin=fieldCoords[0]-numpy.array(aspectRatio)/2.0
# Q8's centre is between measurement points, but corners fall on displacement points, obviously
......
tests/test_scripts.py
View file @ 5a73e117
# -*- coding: utf-8 -*-
from __future__ import print_function
import unittest
import numpy
......@@ -17,8 +16,12 @@ import math
import spam.kalisphera
import coverage
# This for hiding script output
FNULL = open(os.devnull, 'w')
# This for showing it
#FNULL = None
# We also check 2D slice so don't displace in Z or rotate around an axis other than Z
refTranslation = [0.0, 2.0, 0.0]
refRotation = [5.0, 0.0, 0.0]
......@@ -159,12 +162,14 @@ class testAll(unittest.TestCase):
## OK onto the regularStrain
#######################################################
# make sure it runs the help without error
exitCode = subprocess.call(["spam-regularStrain", "--help"], stdout=FNULL, stderr=FNULL)
exitCode = subprocess.call(["spam-regularStrain", "--help"],
stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
# Now let's calculate rotation vector with Geers elements (radius=1) and TSV output
exitCode = subprocess.call(["spam-regularStrain", "snow-ref-snow-def.tsv", "-comp", "r", "-tsv", "-mask", "-r", "1"],
stdout=FNULL, stderr=FNULL)
exitCode = subprocess.call(["spam-regularStrain", "snow-ref-snow-def.tsv", "-comp", "r", "-mask", "-r", "1"],
stdout=FNULL, stderr=FNULL
)
self.assertEqual(exitCode, 0)
# Now read output
......@@ -180,7 +185,7 @@ class testAll(unittest.TestCase):
# Now let's calculate small and large strains with Q8 elements and TSV output
exitCode = subprocess.call(["spam-regularStrain", "--Q8", "snow-ref-snow-def.tsv",
"-comp", "U", "e", "vol", "volss", "dev", "devss",
"-tsv", "-vtk", "-mask"],
"-vtk", "-mask"],
stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
strain = spam.helpers.readStrainTSV("snow-ref-snow-def-strain-Q8.tsv")
......@@ -496,7 +501,9 @@ class testAll(unittest.TestCase):
testFolder + "Step1.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step1-discreteDVC.tsv",
"-pfd"])
"-pfd"],
stdout=FNULL,
stderr=FNULL)
self.assertEqual(exitCode, 0)
#Check number of converged grains
......@@ -563,7 +570,9 @@ class testAll(unittest.TestCase):
testFolder + "Step0.tif",
testFolder + "Lab0.tif",
testFolder + "Step1.tif",
"-od", testFolder + ""])
"-od", testFolder + ""],
stdout=FNULL,
stderr=FNULL)
self.assertEqual(exitCode, 0)
# Check TSV for local results
......@@ -601,7 +610,9 @@ class testAll(unittest.TestCase):
testFolder + "Step0.tif",
testFolder + "Lab0.tif",
testFolder + "Step2.tif",
"-od", testFolder + ""])
"-od", testFolder + ""],
stdout=FNULL,
stderr=FNULL)
self.assertEqual(exitCode, 0)
# Check TSV for local results
......@@ -642,7 +653,9 @@ class testAll(unittest.TestCase):
testFolder + "Step2.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv",
"-pfd"])
"-pfd"],
stdout=FNULL,
stderr=FNULL)
self.assertEqual(exitCode, 0)
# Check TSV for local results
......@@ -680,7 +693,9 @@ class testAll(unittest.TestCase):
testFolder + "Step2.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv",
"-pfd", "-skp"])
"-pfd", "-skp"],
stdout=FNULL,
stderr=FNULL)
# Read new TSV
TSV2 = spam.helpers.readCorrelationTSV(testFolder + "Step0-Step2-discreteDVC.tsv")
# Compare "deltaPhiNorm" for all but the last one (they should be the same since they were skipped)
......@@ -706,7 +721,9 @@ class testAll(unittest.TestCase):
testFolder + "extreme-im1lab.tif",
testFolder + "extreme-im2.tif",
"-regbb", "4", "-regbe", "2",
"-lcm", "10", "-ld", "2"])
"-lcm", "10", "-ld", "2"],
stdout=FNULL,
stderr=FNULL)
self.assertEqual(exitCode, 0)
TSVextreme = spam.helpers.readCorrelationTSV(testFolder + "extreme-im1-extreme-im2-discreteDVC.tsv")
#self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
......@@ -739,7 +756,8 @@ class testAll(unittest.TestCase):
testFolder + "test.vtk",
'-pre', 'test-strain', '-comp', 'dev', 'vol', 'U'],
stdout=FNULL,
stderr=FNULL)
stderr=FNULL
)
self.assertEqual(exitCode, 0)
VTK = spam.helpers.readUnstructuredVTK(testFolder + "test-strain.vtk")
self.assertAlmostEqual(VTK[3]['dev'].mean(), 0, places=3)
......
tools/deformation/deformationField.py
View file @ 5a73e117
......@@ -62,20 +62,25 @@ def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault
dims = list(displacementField.shape[0:3])
# Q8 has 1 element fewer than the number of displacement points
nCells = [n - 1 for n in dims]
cellDims = [n - 1 for n in dims]
# Check if a 2D field is passed
if dims[0] == 1:
# Add a ficticious layer of nodes and cells in Z direction
dims[0] += 1
nCells[0] += 1
cellDims[0] += 1
nodeSpacing[0] += 1
# Add a ficticious layer of equal displacements so that the strain in z is null
displacementField = numpy.concatenate((displacementField, displacementField))
numberOfCells = cellDims[0]*cellDims[1]*cellDims[2]
dims = tuple(dims)
cellDims = tuple(cellDims)
# Transformation gradient tensor F = du/dx +I
Ffield = numpy.zeros((nCells[0], nCells[1], nCells[2], 3, 3))
Ffield = numpy.zeros((cellDims[0], cellDims[1], cellDims[2], 3, 3))
FfieldFlat = Ffield.reshape((numberOfCells, 3, 3))
# Define the coordinates of the Parent Element
# we're using isoparametric Q8 elements
......@@ -104,71 +109,52 @@ def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault
jacY = 2.0 / float(nodeSpacing[1])
jacX = 2.0 / float(nodeSpacing[2])
if verbose: bar = progressbar.ProgressBar(maxval=nCells[0]*nCells[1]*nCells[2]).start()
nodesDone = 0
## Loop over the cells
#global computeConvexVolume
#def computeConvexVolume(label):
#labelI = spam.label.getLabel(lab, label, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
#subvol = labelI['subvol']
#points = numpy.transpose(numpy.where(subvol))
#try:
#hull = scipy.spatial.ConvexHull(points)
#deln = scipy.spatial.Delaunay(points[hull.vertices])
#idx = numpy.stack(numpy.indices(subvol.shape), axis = -1)
#out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)
#hullIm = numpy.zeros(subvol.shape)
#hullIm[out_idx] = 1
#hullVol = spam.label.volumes(hullIm)
#return label, hullVol[-1]
#except:
#return label, 0
## Run multiprocessing
#with multiprocessing.Pool(processes=nProcesses) as pool:
#for returns in pool.imap_unordered(computeConvexVolume, range(1, nLabels+1)):
#if verbose:
#finishedNodes += 1
#widgets[0] = progressbar.FormatLabel(" vol={:0>7.5f} ".format(returns[1]))
#pbar.update(finishedNodes)
#convexVolume[returns[0]] = returns[1]
for kCell in range(nCells[0]):
for jCell in range(nCells[1]):
for iCell in range(nCells[2]):
# Check for nans in one of the 8 nodes of the cell
dCell = displacementField[kCell:kCell + 2, jCell:jCell + 2, iCell:iCell + 2]
if not numpy.all(numpy.isfinite(dCell)):
Ffield[kCell, jCell, iCell, :, :] = numpy.zeros((3, 3)) * numpy.nan
# If no nans start the strain calculation
else:
# Initialise the gradient of the displacement tensor
dudx = numpy.zeros((3, 3))
# Loop over each node of the cell
for node in range(8):
# Get the displacement value
d = displacementField[int(kCell + lid[node, 0]), int(jCell + lid[node, 1]), int(iCell + lid[node, 2]), :]
# Compute the influence of each node to the displacement gradient tensor
dudx[0, 0] += jacZ * SFderivative[node, 0] * d[0]
dudx[1, 1] += jacY * SFderivative[node, 1] * d[1]
dudx[2, 2] += jacX * SFderivative[node, 2] * d[2]
dudx[1, 0] += jacY * SFderivative[node, 1] * d[0]
dudx[0, 1] += jacZ * SFderivative[node, 0] * d[1]
dudx[2, 1] += jacX * SFderivative[node, 2] * d[1]
dudx[1, 2] += jacY * SFderivative[node, 1] * d[2]
dudx[2, 0] += jacX * SFderivative[node, 2] * d[0]
dudx[0, 2] += jacZ * SFderivative[node, 0] * d[2]
# Adding a transpose to dudx, it's ugly but allows us to pass #142
Ffield[kCell, jCell, iCell, :, :] = numpy.eye(3)+dudx.T
bar.update(nodesDone)
nodesDone += 1
bar.finish()
if verbose:
pbar = progressbar.ProgressBar(maxval=numberOfCells).start()
finishedCells = 0
# Loop over the cells
global computeOneQ8
def computeOneQ8(cell):
zCell, yCell, xCell = numpy.unravel_index(cell, cellDims)
# Check for nans in one of the 8 nodes of the cell
if not numpy.all(numpy.isfinite(displacementField[zCell:zCell + 2, yCell:yCell + 2, xCell:xCell + 2])):
F = numpy.zeros((3, 3)) * numpy.nan
# If no nans start the strain calculation
else:
# Initialise the gradient of the displacement tensor
dudx = numpy.zeros((3, 3))
# Loop over each node of the cell
for node in range(8):
# Get the displacement value
d = displacementField[int(zCell + lid[node, 0]), int(yCell + lid[node, 1]), int(xCell + lid[node, 2]), :]
# Compute the influence of each node to the displacement gradient tensor
dudx[0, 0] += jacZ * SFderivative[node, 0] * d[0]
dudx[1, 1] += jacY * SFderivative[node, 1] * d[1]
dudx[2, 2] += jacX * SFderivative[node, 2] * d[2]
dudx[1, 0] += jacY * SFderivative[node, 1] * d[0]
dudx[0, 1] += jacZ * SFderivative[node, 0] * d[1]
dudx[2, 1] += jacX * SFderivative[node, 2] * d[1]
dudx[1, 2] += jacY * SFderivative[node, 1] * d[2]
dudx[2, 0] += jacX * SFderivative[node, 2] * d[0]
dudx[0, 2] += jacZ * SFderivative[node, 0] * d[2]
# Adding a transpose to dudx, it's ugly but allows us to pass #142
F = numpy.eye(3)+dudx.T
return cell, F
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(computeOneQ8, range(numberOfCells)):
if verbose:
finishedCells += 1
pbar.update(finishedCells)
FfieldFlat[returns[0]] = returns[1]
if verbose: pbar.finish()
return Ffield
......@@ -346,9 +332,9 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
return goodPoint, F
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldCoordsFlatGood.shape[0]).start()
finishedPoints = 0
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:
......@@ -364,7 +350,7 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
def FfieldBagi(points, connectivity, displacements):
def FfieldBagi(points, connectivity, displacements, nProcesses=nProcessesDefault, verbose=True):
"""
Calculates transformation gradient function using Bagi's 1996 paper, especially equation 3 on page 174.
Required inputs are connectivity matrix for tetrahedra (for example from spam.mesh.triangulate) and
......@@ -381,14 +367,20 @@ def FfieldBagi(points, connectivity, displacements):
displacements : m x 3 numpy array
M Particles' displacement
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
Ffield: nx3x3 array of n cells
The field of the transformation gradient tensors
"""
from progressbar import ProgressBar
import spam.mesh
pbar = ProgressBar()
Ffield = numpy.zeros([connectivity.shape[0], 3, 3], dtype='<f4')
......@@ -413,22 +405,23 @@ def FfieldBagi(points, connectivity, displacements):
# print("spam.mesh.bagiStrain(): Constructing strain from Delaunay and Displacements")
# Loop through tetrahdra to get avec1, uPos1
# for tet in range(connectivity.shape[0]):
for tet in pbar(range(connectivity.shape[0])):
global computeOneTet
def computeOneTet(tet):
# Get the list of IDs, centroids, center of tet
tet_ids = connectivity[tet, :]
tetIDs = connectivity[tet, :]
# 2019-10-07 EA: Skip references to missing particles
if max(tet_ids) >= points.shape[0]:
print("spam.mesh.unstructured.bagiStrain(): this tet has node > points.shape[0], skipping")
pass
else:
tetCoords = points[tet_ids, :]
tetDisp = displacements[tet_ids, :]
#if max(tetIDs) >= points.shape[0]:
#print("spam.mesh.unstructured.bagiStrain(): this tet has node > points.shape[0], skipping")
#pass
#else:
if True:
tetCoords = points[tetIDs, :]
tetDisp = displacements[tetIDs, :]
tetCen = numpy.average(tetCoords, axis=0)
if numpy.isfinite(tetCoords).sum() + numpy.isfinite(tetDisp).sum() != 3*4*2:
print("spam.mesh.unstructured.bagiStrain(): nans in position or displacement, skipping")
if verbose: print("spam.mesh.unstructured.bagiStrain(): nans in position or displacement, skipping")
# Compute strains
Ffield[tet] = numpy.zeros((3,3))*numpy.nan
F = numpy.zeros((3,3))*numpy.nan
else:
# Loop through each face of tet to get avec, upos (Bagi, 1996, pg. 172)
# aVec1 = numpy.zeros([4, 3], dtype='<f4')
......@@ -460,7 +453,24 @@ def FfieldBagi(points, connectivity, displacements):
dudx /= float(tetVolumes[tet])
Ffield[tet] = numpy.eye(3) + dudx
F = numpy.eye(3) + dudx
return tet, F
# Iterate through flat field of Fs
if verbose:
pbar = progressbar.ProgressBar(maxval=connectivity.shape[0]).start()
finishedTets = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(computeOneTet, range(connectivity.shape[0])):
if verbose:
finishedTets += 1
pbar.update(finishedTets)
Ffield[returns[0]] = returns[1]
if verbose: pbar.finish()
return Ffield
......@@ -530,9 +540,9 @@ def decomposeFfield(Ffield, components, twoD=False, nProcesses=nProcessesDefault
'devss': numpy.nan}
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
if verbose:
pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
......@@ -626,9 +636,9 @@ def decomposePhiField(PhiField, components, twoD=False, nProcesses=nProcessesDef
'devss': numpy.nan}
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
if verbose:
pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
......
tools/deformation/deformationFunction.py
View file @ 5a73e117
......@@ -144,7 +144,7 @@ def computePhi(transformation, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.
###########################################################
# Polar Decomposition of a given F into human readable components
###########################################################
def decomposeF(F, twoD=False):
def decomposeF(F, twoD=False, verbose=False):
"""
Get components out of a transformation gradient tensor F
......@@ -157,6 +157,10 @@ def decomposeF(F, twoD=False):
is the F defined in 2D? This applies corrections to the strain outputs
Default = False
verbose : boolean, optional
Print errors?
Default = True
Returns
-------
transformation : dictionary of arrays
......@@ -186,19 +190,19 @@ def decomposeF(F, twoD=False):
# Check for NaNs if any quit
if numpy.isnan(F).sum() > 0:
print("deformationFunction.decomposeF(): Nan value in F. Exiting")
if verbose: print("deformationFunction.decomposeF(): Nan value in F. Exiting")
return transformation
# Check for inf if any quit
if numpy.isinf(F).sum() > 0:
print("deformationFunction.decomposeF(): Inf value in F. Exiting")
if verbose: print("deformationFunction.decomposeF(): Inf value in F. Exiting")
return transformation
# Check for singular F if yes quit
try:
numpy.linalg.inv(F)
except numpy.linalg.linalg.LinAlgError:
print("deformationFunction.decomposeF(): F is singular. Exiting")
if verbose: print("deformationFunction.decomposeF(): F is singular. Exiting")
return transformation
###########################################################
......@@ -335,7 +339,7 @@ def decomposeF(F, twoD=False):
return transformation
def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=False):
def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=False, verbose=False):
"""
Get components out of a linear deformation function "Phi"
......@@ -354,6 +358,10 @@ def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=
is the Phi defined in 2D? This applies corrections to the strain outputs
Default = False
verbose : boolean, optional
Print errors?
Default = True
Returns
-------
transformation : dictionary of arrays
......@@ -414,7 +422,7 @@ def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=
# apply Phi to the given point and calculate its displacement
tra -= dist - numpy.dot(F, dist)
decomposedF = decomposeF(F)
decomposedF = decomposeF(F, verbose=verbose)