Coverage for /usr/local/lib/python3.8/site-packages/spam/deformation/deformationField.py: 60%
262 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1# Library of SPAM functions for dealing with fields of Phi or fields of F
2# Copyright (C) 2020 SPAM Contributors
3#
4# This program is free software: you can redistribute it and/or modify it
5# under the terms of the GNU General Public License as published by the Free
6# Software Foundation, either version 3 of the License, or (at your option)
7# any later version.
8#
9# This program is distributed in the hope that it will be useful, but WITHOUT
10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12# more details.
13#
14# You should have received a copy of the GNU General Public License along with
15# this program. If not, see <http://www.gnu.org/licenses/>.
17import multiprocessing
19import numpy
20import progressbar
22# 2020-02-24 Olga Stamati and Edward Ando
23import spam.deformation
24import spam.label
26# Global number of processes
27nProcessesDefault = multiprocessing.cpu_count()
30def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault, verbose=False):
31 """
32 This function computes the transformation gradient field F from a given displacement field.
33 Please note: the transformation gradient tensor: F = I + du/dx.
35 This function computes du/dx in the centre of an 8-node cell (Q8 in Finite Elements terminology) using order one (linear) shape functions.
37 Parameters
38 ----------
39 displacementField : 4D array of floats
40 The vector field to compute the derivatives.
41 #Its shape is (nz, ny, nx, 3)
43 nodeSpacing : 3-component list of floats
44 Length between two nodes in every direction (*i.e.,* size of a cell)
46 nProcesses : integer, optional
47 Number of processes for multiprocessing
48 Default = number of CPUs in the system
50 verbose : boolean, optional
51 Print progress?
52 Default = True
54 Returns
55 -------
56 F : (nz-1) x (ny-1) x (nx-1) x 3x3 array of n cells
57 The field of the transformation gradient tensors
58 """
59 # import spam.DIC.deformationFunction
60 # import spam.mesh.strain
62 # Define dimensions
63 dims = list(displacementField.shape[0:3])
65 # Q8 has 1 element fewer than the number of displacement points
66 cellDims = [n - 1 for n in dims]
68 # Check if a 2D field is passed
69 if dims[0] == 1:
70 # Add a ficticious layer of nodes and cells in Z direction
71 dims[0] += 1
72 cellDims[0] += 1
73 nodeSpacing[0] += 1
75 # Add a ficticious layer of equal displacements so that the strain in z is null
76 displacementField = numpy.concatenate((displacementField, displacementField))
78 numberOfCells = cellDims[0] * cellDims[1] * cellDims[2]
79 dims = tuple(dims)
80 cellDims = tuple(cellDims)
82 # Transformation gradient tensor F = du/dx +I
83 Ffield = numpy.zeros((cellDims[0], cellDims[1], cellDims[2], 3, 3))
84 FfieldFlat = Ffield.reshape((numberOfCells, 3, 3))
86 # Define the coordinates of the Parent Element
87 # we're using isoparametric Q8 elements
88 lid = numpy.zeros((8, 3)).astype("<u1") # local index
89 lid[0] = [0, 0, 0]
90 lid[1] = [0, 0, 1]
91 lid[2] = [0, 1, 0]
92 lid[3] = [0, 1, 1]
93 lid[4] = [1, 0, 0]
94 lid[5] = [1, 0, 1]
95 lid[6] = [1, 1, 0]
96 lid[7] = [1, 1, 1]
98 # Calculate the derivatives of the shape functions
99 # Since the center is equidistant from all 8 nodes, each one gets equal weighting
100 SFderivative = numpy.zeros((8, 3))
101 for node in range(8):
102 # (local nodes coordinates) / weighting of each node
103 SFderivative[node, 0] = (2.0 * (float(lid[node, 0]) - 0.5)) / 8.0
104 SFderivative[node, 1] = (2.0 * (float(lid[node, 1]) - 0.5)) / 8.0
105 SFderivative[node, 2] = (2.0 * (float(lid[node, 2]) - 0.5)) / 8.0
107 # Compute the jacobian to go from local(Parent Element) to global base
108 jacZ = 2.0 / float(nodeSpacing[0])
109 jacY = 2.0 / float(nodeSpacing[1])
110 jacX = 2.0 / float(nodeSpacing[2])
112 if verbose:
113 pbar = progressbar.ProgressBar(maxval=numberOfCells).start()
114 finishedCells = 0
116 # Loop over the cells
117 global computeOneQ8
119 def computeOneQ8(cell):
120 zCell, yCell, xCell = numpy.unravel_index(cell, cellDims)
122 # Check for nans in one of the 8 nodes of the cell
123 if not numpy.all(numpy.isfinite(displacementField[zCell : zCell + 2, yCell : yCell + 2, xCell : xCell + 2])):
124 F = numpy.zeros((3, 3)) * numpy.nan
126 # If no nans start the strain calculation
127 else:
128 # Initialise the gradient of the displacement tensor
129 dudx = numpy.zeros((3, 3))
131 # Loop over each node of the cell
132 for node in range(8):
133 # Get the displacement value
134 d = displacementField[
135 int(zCell + lid[node, 0]),
136 int(yCell + lid[node, 1]),
137 int(xCell + lid[node, 2]),
138 :,
139 ]
141 # Compute the influence of each node to the displacement gradient tensor
142 dudx[0, 0] += jacZ * SFderivative[node, 0] * d[0]
143 dudx[1, 1] += jacY * SFderivative[node, 1] * d[1]
144 dudx[2, 2] += jacX * SFderivative[node, 2] * d[2]
145 dudx[1, 0] += jacY * SFderivative[node, 1] * d[0]
146 dudx[0, 1] += jacZ * SFderivative[node, 0] * d[1]
147 dudx[2, 1] += jacX * SFderivative[node, 2] * d[1]
148 dudx[1, 2] += jacY * SFderivative[node, 1] * d[2]
149 dudx[2, 0] += jacX * SFderivative[node, 2] * d[0]
150 dudx[0, 2] += jacZ * SFderivative[node, 0] * d[2]
151 # Adding a transpose to dudx, it's ugly but allows us to pass #142
152 F = numpy.eye(3) + dudx.T
153 return cell, F
155 # Run multiprocessing
156 with multiprocessing.Pool(processes=nProcesses) as pool:
157 for returns in pool.imap_unordered(computeOneQ8, range(numberOfCells)):
158 if verbose:
159 finishedCells += 1
160 pbar.update(finishedCells)
161 FfieldFlat[returns[0]] = returns[1]
162 pool.close()
163 pool.join()
165 if verbose:
166 pbar.finish()
168 return Ffield
171def FfieldRegularGeers(
172 displacementField,
173 nodeSpacing,
174 neighbourRadius=1.5,
175 nProcesses=nProcessesDefault,
176 verbose=False,
177):
178 """
179 This function computes the transformation gradient field F from a given displacement field.
180 Please note: the transformation gradient tensor: F = I + du/dx.
182 This function computes du/dx as a weighted function of neighbouring points.
183 Here is implemented the linear model proposed in:
184 "Computing strain fields from discrete displacement fields in 2D-solids", Geers et al., 1996
186 Parameters
187 ----------
188 displacementField : 4D array of floats
189 The vector field to compute the derivatives.
190 Its shape is (nz, ny, nx, 3).
192 nodeSpacing : 3-component list of floats
193 Length between two nodes in every direction (*i.e.,* size of a cell)
195 neighbourRadius : float, optional
196 Distance in nodeSpacings to include neighbours in the strain calcuation.
197 Default = 1.5*nodeSpacing which will give radius = 1.5*min(nodeSpacing)
199 mask : bool, optional
200 Avoid non-correlated NaN points in the displacement field?
201 Default = True
203 nProcesses : integer, optional
204 Number of processes for multiprocessing
205 Default = number of CPUs in the system
207 verbose : boolean, optional
208 Print progress?
209 Default = True
211 Returns
212 -------
213 Ffield: nz x ny x nx x 3x3 array of n cells
214 The field of the transformation gradient tensors
216 Note
217 ----
218 Taken from the implementation in "TomoWarp2: A local digital volume correlation code", Tudisco et al., 2017
219 """
220 import scipy.spatial
222 # Define dimensions
223 dims = displacementField.shape[0:3]
224 nNodes = dims[0] * dims[1] * dims[2]
225 displacementFieldFlat = displacementField.reshape(nNodes, 3)
227 # Check if a 2D field is passed
228 twoD = False
229 if dims[0] == 1:
230 twoD = True
232 # Deformation gradient tensor F = du/dx +I
233 # Ffield = numpy.zeros((dims[0], dims[1], dims[2], 3, 3))
234 FfieldFlat = numpy.zeros((nNodes, 3, 3))
236 if twoD:
237 fieldCoordsFlat = (
238 numpy.mgrid[
239 0:1:1,
240 nodeSpacing[1] : dims[1] * nodeSpacing[1] + nodeSpacing[1] : nodeSpacing[1],
241 nodeSpacing[2] : dims[2] * nodeSpacing[2] + nodeSpacing[2] : nodeSpacing[2],
242 ]
243 .reshape(3, nNodes)
244 .T
245 )
246 else:
247 fieldCoordsFlat = (
248 numpy.mgrid[
249 nodeSpacing[0] : dims[0] * nodeSpacing[0] + nodeSpacing[0] : nodeSpacing[0],
250 nodeSpacing[1] : dims[1] * nodeSpacing[1] + nodeSpacing[1] : nodeSpacing[1],
251 nodeSpacing[2] : dims[2] * nodeSpacing[2] + nodeSpacing[2] : nodeSpacing[2],
252 ]
253 .reshape(3, nNodes)
254 .T
255 )
257 # Get non-nan displacements
258 goodPointsMask = numpy.isfinite(displacementField[:, :, :, 0].reshape(nNodes))
259 badPointsMask = numpy.isnan(displacementField[:, :, :, 0].reshape(nNodes))
260 # Flattened variables
261 fieldCoordsFlatGood = fieldCoordsFlat[goodPointsMask]
262 displacementFieldFlatGood = displacementFieldFlat[goodPointsMask]
263 # set bad points to nan
264 FfieldFlat[badPointsMask] = numpy.eye(3) * numpy.nan
266 # build KD-tree for neighbour identification
267 treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
269 # Output array for good points
270 FfieldFlatGood = numpy.zeros_like(FfieldFlat[goodPointsMask])
272 # Function for parallel mode
273 global geersOnePoint
275 def geersOnePoint(goodPoint):
276 # This is for the linear model, equation 15 in Geers
277 centralNodePosition = fieldCoordsFlatGood[goodPoint]
278 centralNodeDisplacement = displacementFieldFlatGood[goodPoint]
279 sX0X0 = numpy.zeros((3, 3))
280 sX0Xt = numpy.zeros((3, 3))
281 m0 = numpy.zeros(3)
282 mt = numpy.zeros(3)
284 # Option 2: KDTree on distance
285 # KD-tree will always give the current point as zero-distance
286 ind = treeCoord.query_ball_point(centralNodePosition, neighbourRadius * max(nodeSpacing))
288 # We know that the current point will also be included, so remove it from the index list.
289 ind = numpy.array(ind)
290 ind = ind[ind != goodPoint]
291 nNeighbours = len(ind)
292 nodalRelativePositionsRef = numpy.zeros((nNeighbours, 3)) # Delta_X_0 in paper
293 nodalRelativePositionsDef = numpy.zeros((nNeighbours, 3)) # Delta_X_t in paper
295 for neighbour, i in enumerate(ind):
296 # Relative position in reference configuration
297 # absolute position of this neighbour node
298 # minus abs position of central node
299 nodalRelativePositionsRef[neighbour, :] = fieldCoordsFlatGood[i] - centralNodePosition
301 # Relative position in deformed configuration (i.e., plus displacements)
302 # absolute position of this neighbour node
303 # plus displacement of this neighbour node
304 # minus abs position of central node
305 # minus displacement of central node
306 nodalRelativePositionsDef[neighbour, :] = fieldCoordsFlatGood[i] + displacementFieldFlatGood[i] - centralNodePosition - centralNodeDisplacement
308 for u in range(3):
309 for v in range(3):
310 # sX0X0[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
311 # sX0Xt[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
312 # Proposed solution for #142 for direction of rotation
313 sX0X0[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
314 sX0Xt[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
316 m0 += nodalRelativePositionsRef[neighbour, :]
317 mt += nodalRelativePositionsDef[neighbour, :]
319 sX0X0 = nNeighbours * sX0X0
320 sX0Xt = nNeighbours * sX0Xt
322 A = sX0X0 - numpy.dot(m0, m0)
323 C = sX0Xt - numpy.dot(m0, mt)
324 F = numpy.zeros((3, 3))
326 if twoD:
327 A = A[1:, 1:]
328 C = C[1:, 1:]
329 try:
330 F[1:, 1:] = numpy.dot(numpy.linalg.inv(A), C)
331 F[0, 0] = 1.0
332 except numpy.linalg.linalg.LinAlgError:
333 # print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
334 pass
335 else:
336 try:
337 F = numpy.dot(numpy.linalg.inv(A), C)
338 except numpy.linalg.linalg.LinAlgError:
339 # print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
340 pass
342 return goodPoint, F
344 # Iterate through flat field of Fs
345 if verbose:
346 pbar = progressbar.ProgressBar(maxval=fieldCoordsFlatGood.shape[0]).start()
347 finishedPoints = 0
349 # Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat
350 with multiprocessing.Pool(processes=nProcesses) as pool:
351 for returns in pool.imap_unordered(geersOnePoint, range(fieldCoordsFlatGood.shape[0])):
352 if verbose:
353 finishedPoints += 1
354 pbar.update(finishedPoints)
355 FfieldFlatGood[returns[0]] = returns[1]
357 if verbose:
358 pbar.finish()
360 FfieldFlat[goodPointsMask] = FfieldFlatGood
361 return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
364def FfieldBagi(points, connectivity, displacements, nProcesses=nProcessesDefault, verbose=False):
365 """
366 Calculates transformation gradient function using Bagi's 1996 paper, especially equation 3 on page 174.
367 Required inputs are connectivity matrix for tetrahedra (for example from spam.mesh.triangulate) and
368 nodal positions in reference and deformed configurations.
370 Parameters
371 ----------
372 points : m x 3 numpy array
373 M Particles' points in reference configuration
375 connectivity : n x 4 numpy array
376 Delaunay triangulation connectivity generated by spam.mesh.triangulate for example
378 displacements : m x 3 numpy array
379 M Particles' displacement
381 nProcesses : integer, optional
382 Number of processes for multiprocessing
383 Default = number of CPUs in the system
385 verbose : boolean, optional
386 Print progress?
387 Default = True
389 Returns
390 -------
391 Ffield: nx3x3 array of n cells
392 The field of the transformation gradient tensors
393 """
394 import spam.mesh
396 Ffield = numpy.zeros([connectivity.shape[0], 3, 3], dtype="<f4")
398 connectivity = connectivity.astype(numpy.uint)
400 # define dimension
401 # D = 3.0
403 # Import modules
405 # Construct 4-list of 3-lists of combinations constituting a face of the tet
406 combs = [[0, 1, 2], [1, 2, 3], [2, 3, 0], [0, 1, 3]]
407 unode = [3, 0, 1, 2]
409 # Precompute tetrahedron Volumes
410 tetVolumes = spam.mesh.tetVolumes(points, connectivity)
412 # Initialize arrays for tet strains
413 # print("spam.mesh.bagiStrain(): Constructing strain from Delaunay and Displacements")
415 # Loop through tetrahdra to get avec1, uPos1
416 global computeOneTet
418 def computeOneTet(tet):
419 # Get the list of IDs, centroids, center of tet
420 tetIDs = connectivity[tet, :]
421 # 2019-10-07 EA: Skip references to missing particles
422 # if max(tetIDs) >= points.shape[0]:
423 # print("spam.mesh.unstructured.bagiStrain(): this tet has node > points.shape[0], skipping")
424 # pass
425 # else:
426 if True:
427 tetCoords = points[tetIDs, :]
428 tetDisp = displacements[tetIDs, :]
429 tetCen = numpy.average(tetCoords, axis=0)
430 if numpy.isfinite(tetCoords).sum() + numpy.isfinite(tetDisp).sum() != 3 * 4 * 2:
431 if verbose:
432 print("spam.mesh.unstructured.bagiStrain(): nans in position or displacement, skipping")
433 # Compute strains
434 F = numpy.zeros((3, 3)) * numpy.nan
435 else:
436 # Loop through each face of tet to get avec, upos (Bagi, 1996, pg. 172)
437 # aVec1 = numpy.zeros([4, 3], dtype='<f4')
438 # uPos1 = numpy.zeros([4, 3], dtype='<f4')
439 # uPos2 = numpy.zeros([4, 3], dtype='<f4')
440 dudx = numpy.zeros((3, 3), dtype="<f4")
442 for face in range(4):
443 faceNorm = numpy.cross(
444 tetCoords[combs[face][0]] - tetCoords[combs[face][1]],
445 tetCoords[combs[face][0]] - tetCoords[combs[face][2]],
446 )
448 # Get a norm vector to face point towards center of tet
449 faceCen = numpy.average(tetCoords[combs[face]], axis=0)
450 tmpnorm = faceNorm / (numpy.linalg.norm(faceNorm))
451 facetocen = tetCen - faceCen
452 if numpy.dot(facetocen, tmpnorm) < 0:
453 tmpnorm = -tmpnorm
455 # Divide by 6 (1/3 for 1/Dimension; 1/2 for area from cross product)
456 # See first eqn., Bagi, 1996, pg. 172.
457 # aVec1[face] = tmpnorm*numpy.linalg.norm(faceNorm)/6
459 # Undeformed positions
460 # uPos1[face] = tetCoords[unode[face]]
461 # Deformed positions
462 # uPos2[face] = tetComs2[unode[face]]
464 dudx += numpy.tensordot(
465 tetDisp[unode[face]],
466 tmpnorm * numpy.linalg.norm(faceNorm) / 6,
467 axes=0,
468 )
470 dudx /= float(tetVolumes[tet])
472 F = numpy.eye(3) + dudx
473 return tet, F
475 # Iterate through flat field of Fs
476 if verbose:
477 pbar = progressbar.ProgressBar(maxval=connectivity.shape[0]).start()
478 finishedTets = 0
480 # Run multiprocessing
481 with multiprocessing.Pool(processes=nProcesses) as pool:
482 for returns in pool.imap_unordered(computeOneTet, range(connectivity.shape[0])):
483 if verbose:
484 finishedTets += 1
485 pbar.update(finishedTets)
486 Ffield[returns[0]] = returns[1]
487 pool.close()
488 pool.join()
490 if verbose:
491 pbar.finish()
493 return Ffield
496def decomposeFfield(Ffield, components, twoD=False, nProcesses=nProcessesDefault, verbose=False):
497 """
498 This function takes in an F field (from either FfieldRegularQ8, FfieldRegularGeers, FfieldBagi) and
499 returns fields of desired transformation components.
501 Parameters
502 ----------
503 Ffield : multidimensional x 3 x 3 numpy array of floats
504 Spatial field of Fs
506 components : list of strings
507 These indicate the desired components consistent with spam.deformation.decomposeF or decomposePhi
509 twoD : bool, optional
510 Is the Ffield in 2D? This changes the strain calculation.
511 Default = False
513 nProcesses : integer, optional
514 Number of processes for multiprocessing
515 Default = number of CPUs in the system
517 verbose : boolean, optional
518 Print progress?
519 Default = True
521 Returns
522 -------
523 Dictionary containing appropriately reshaped fields of the transformation components requested.
525 Keys:
526 vol, dev, volss, devss are scalars
527 r and z are 3-component vectors
528 e and U are 3x3 tensors
529 """
530 # The last two are for the 3x3 F field
531 fieldDimensions = Ffield.shape[0:-2]
532 fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
533 FfieldFlat = Ffield.reshape(fieldRavelLength, 3, 3)
535 output = {}
536 for component in components:
537 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
538 output[component] = numpy.zeros(fieldRavelLength)
539 if component == "r" or component == "z":
540 output[component] = numpy.zeros((fieldRavelLength, 3))
541 if component == "U" or component == "e":
542 output[component] = numpy.zeros((fieldRavelLength, 3, 3))
544 # Function for parallel mode
545 global decomposeOneF
547 def decomposeOneF(n):
548 F = FfieldFlat[n]
549 if numpy.isfinite(F).sum() == 9:
550 decomposedF = spam.deformation.decomposeF(F, twoD=twoD)
551 return n, decomposedF
552 else:
553 return n, {
554 "r": numpy.array([numpy.nan] * 3),
555 "z": numpy.array([numpy.nan] * 3),
556 "U": numpy.eye(3) * numpy.nan,
557 "e": numpy.eye(3) * numpy.nan,
558 "vol": numpy.nan,
559 "dev": numpy.nan,
560 "volss": numpy.nan,
561 "devss": numpy.nan,
562 }
564 # Iterate through flat field of Fs
565 if verbose:
566 pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
567 finishedPoints = 0
569 # Run multiprocessing
570 with multiprocessing.Pool(processes=nProcesses) as pool:
571 for returns in pool.imap_unordered(decomposeOneF, range(fieldRavelLength)):
572 if verbose:
573 finishedPoints += 1
574 pbar.update(finishedPoints)
575 for component in components:
576 output[component][returns[0]] = returns[1][component]
577 pool.close()
578 pool.join()
580 if verbose:
581 pbar.finish()
583 # Reshape on the output
584 for component in components:
585 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
586 output[component] = numpy.array(output[component]).reshape(fieldDimensions)
588 if component == "r" or component == "z":
589 output[component] = numpy.array(output[component]).reshape(Ffield.shape[0:-1])
591 if component == "U" or component == "e":
592 output[component] = numpy.array(output[component]).reshape(Ffield.shape)
594 return output
597def decomposePhiField(PhiField, components, twoD=False, nProcesses=nProcessesDefault, verbose=False):
598 """
599 This function takes in a Phi field (from readCorrelationTSV?) and
600 returns fields of desired transformation components.
602 Parameters
603 ----------
604 PhiField : multidimensional x 4 x 4 numpy array of floats
605 Spatial field of Phis
607 components : list of strings
608 These indicate the desired components consistent with spam.deformation.decomposePhi
610 twoD : bool, optional
611 Is the PhiField in 2D? This changes the strain calculation.
612 Default = False
614 nProcesses : integer, optional
615 Number of processes for multiprocessing
616 Default = number of CPUs in the system
618 verbose : boolean, optional
619 Print progress?
620 Default = True
622 Returns
623 -------
624 Dictionary containing appropriately reshaped fields of the transformation components requested.
626 Keys:
627 vol, dev, volss, devss are scalars
628 t, r, and z are 3-component vectors
629 e and U are 3x3 tensors
630 """
631 # The last two are for the 4x4 Phi field
632 fieldDimensions = PhiField.shape[0:-2]
633 fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
634 PhiFieldFlat = PhiField.reshape(fieldRavelLength, 4, 4)
636 output = {}
637 for component in components:
638 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
639 output[component] = numpy.zeros(fieldRavelLength)
640 if component == "t" or component == "r" or component == "z":
641 output[component] = numpy.zeros((fieldRavelLength, 3))
642 if component == "U" or component == "e":
643 output[component] = numpy.zeros((fieldRavelLength, 3, 3))
645 # Function for parallel mode
646 global decomposeOnePhi
648 def decomposeOnePhi(n):
649 Phi = PhiFieldFlat[n]
650 if numpy.isfinite(Phi).sum() == 16:
651 decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
652 return n, decomposedPhi
653 else:
654 return n, {
655 "t": numpy.array([numpy.nan] * 3),
656 "r": numpy.array([numpy.nan] * 3),
657 "z": numpy.array([numpy.nan] * 3),
658 "U": numpy.eye(3) * numpy.nan,
659 "e": numpy.eye(3) * numpy.nan,
660 "vol": numpy.nan,
661 "dev": numpy.nan,
662 "volss": numpy.nan,
663 "devss": numpy.nan,
664 }
666 # Iterate through flat field of Fs
667 if verbose:
668 pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
669 finishedPoints = 0
671 # Run multiprocessing
672 with multiprocessing.Pool(processes=nProcesses) as pool:
673 for returns in pool.imap_unordered(decomposeOnePhi, range(fieldRavelLength)):
674 if verbose:
675 finishedPoints += 1
676 pbar.update(finishedPoints)
677 for component in components:
678 output[component][returns[0]] = returns[1][component]
679 pool.close()
680 pool.join()
682 if verbose:
683 pbar.finish()
685 # Reshape on the output
686 for component in components:
687 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
688 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2])
690 if component == "t" or component == "r" or component == "z":
691 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3)
693 if component == "U" or component == "e":
694 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3, 3)
696 return output