Coverage for /usr/local/lib/python3.10/site-packages/spam/deformation/deformationField.py: 99%
175 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17: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
19try:
20 multiprocessing.set_start_method("fork")
21except RuntimeError:
22 pass
24import numpy
25import progressbar
27# 2020-02-24 Olga Stamati and Edward Ando
28import spam.deformation
29import spam.label
31# Global number of processes
32nProcessesDefault = multiprocessing.cpu_count()
35def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault, verbose=False):
36 """
37 This function computes the transformation gradient field F from a given displacement field.
38 Please note: the transformation gradient tensor: F = I + du/dx.
40 This function computes du/dx in the centre of an 8-node cell (Q8 in Finite Elements terminology) using order one (linear) shape functions.
42 Parameters
43 ----------
44 displacementField : 4D array of floats
45 The vector field to compute the derivatives.
46 #Its shape is (nz, ny, nx, 3)
48 nodeSpacing : 3-component list of floats
49 Length between two nodes in every direction (*i.e.,* size of a cell)
51 nProcesses : integer, optional
52 Number of processes for multiprocessing
53 Default = number of CPUs in the system
55 verbose : boolean, optional
56 Print progress?
57 Default = True
59 Returns
60 -------
61 F : (nz-1) x (ny-1) x (nx-1) x 3x3 array of n cells
62 The field of the transformation gradient tensors
63 """
64 # import spam.DIC.deformationFunction
65 # import spam.mesh.strain
67 # Define dimensions
68 dims = list(displacementField.shape[0:3])
70 # Q8 has 1 element fewer than the number of displacement points
71 cellDims = [n - 1 for n in dims]
73 # Check if a 2D field is passed
74 if dims[0] == 1:
75 # Add a ficticious layer of nodes and cells in Z direction
76 dims[0] += 1
77 cellDims[0] += 1
78 nodeSpacing[0] += 1
80 # Add a ficticious layer of equal displacements so that the strain in z is null
81 displacementField = numpy.concatenate((displacementField, displacementField))
83 numberOfCells = cellDims[0] * cellDims[1] * cellDims[2]
84 dims = tuple(dims)
85 cellDims = tuple(cellDims)
87 # Transformation gradient tensor F = du/dx +I
88 Ffield = numpy.zeros((cellDims[0], cellDims[1], cellDims[2], 3, 3))
89 FfieldFlat = Ffield.reshape((numberOfCells, 3, 3))
91 # Define the coordinates of the Parent Element
92 # we're using isoparametric Q8 elements
93 lid = numpy.zeros((8, 3)).astype("<u1") # local index
94 lid[0] = [0, 0, 0]
95 lid[1] = [0, 0, 1]
96 lid[2] = [0, 1, 0]
97 lid[3] = [0, 1, 1]
98 lid[4] = [1, 0, 0]
99 lid[5] = [1, 0, 1]
100 lid[6] = [1, 1, 0]
101 lid[7] = [1, 1, 1]
103 # Calculate the derivatives of the shape functions
104 # Since the center is equidistant from all 8 nodes, each one gets equal weighting
105 SFderivative = numpy.zeros((8, 3))
106 for node in range(8):
107 # (local nodes coordinates) / weighting of each node
108 SFderivative[node, 0] = (2.0 * (float(lid[node, 0]) - 0.5)) / 8.0
109 SFderivative[node, 1] = (2.0 * (float(lid[node, 1]) - 0.5)) / 8.0
110 SFderivative[node, 2] = (2.0 * (float(lid[node, 2]) - 0.5)) / 8.0
112 # Compute the jacobian to go from local(Parent Element) to global base
113 jacZ = 2.0 / float(nodeSpacing[0])
114 jacY = 2.0 / float(nodeSpacing[1])
115 jacX = 2.0 / float(nodeSpacing[2])
117 if verbose:
118 pbar = progressbar.ProgressBar(maxval=numberOfCells).start()
119 finishedCells = 0
121 # Loop over the cells
122 global _multiprocessingComputeOneQ8
124 def _multiprocessingComputeOneQ8(cell):
125 zCell, yCell, xCell = numpy.unravel_index(cell, cellDims)
127 # Check for nans in one of the 8 nodes of the cell
128 if not numpy.all(numpy.isfinite(displacementField[zCell : zCell + 2, yCell : yCell + 2, xCell : xCell + 2])):
129 F = numpy.zeros((3, 3)) * numpy.nan
131 # If no nans start the strain calculation
132 else:
133 # Initialise the gradient of the displacement tensor
134 dudx = numpy.zeros((3, 3))
136 # Loop over each node of the cell
137 for node in range(8):
138 # Get the displacement value
139 d = displacementField[
140 int(zCell + lid[node, 0]),
141 int(yCell + lid[node, 1]),
142 int(xCell + lid[node, 2]),
143 :,
144 ]
146 # Compute the influence of each node to the displacement gradient tensor
147 dudx[0, 0] += jacZ * SFderivative[node, 0] * d[0]
148 dudx[1, 1] += jacY * SFderivative[node, 1] * d[1]
149 dudx[2, 2] += jacX * SFderivative[node, 2] * d[2]
150 dudx[1, 0] += jacY * SFderivative[node, 1] * d[0]
151 dudx[0, 1] += jacZ * SFderivative[node, 0] * d[1]
152 dudx[2, 1] += jacX * SFderivative[node, 2] * d[1]
153 dudx[1, 2] += jacY * SFderivative[node, 1] * d[2]
154 dudx[2, 0] += jacX * SFderivative[node, 2] * d[0]
155 dudx[0, 2] += jacZ * SFderivative[node, 0] * d[2]
156 # Adding a transpose to dudx, it's ugly but allows us to pass #142
157 F = numpy.eye(3) + dudx.T
158 return cell, F
160 # Run multiprocessing
161 with multiprocessing.Pool(processes=nProcesses) as pool:
162 for returns in pool.imap_unordered(_multiprocessingComputeOneQ8, range(numberOfCells)):
163 if verbose:
164 finishedCells += 1
165 pbar.update(finishedCells)
166 FfieldFlat[returns[0]] = returns[1]
167 pool.close()
168 pool.join()
170 if verbose:
171 pbar.finish()
173 return Ffield
176def FfieldRegularGeers(
177 displacementField,
178 nodeSpacing,
179 neighbourRadius=1.5,
180 nProcesses=nProcessesDefault,
181 verbose=False,
182):
183 """
184 This function computes the transformation gradient field F from a given displacement field.
185 Please note: the transformation gradient tensor: F = I + du/dx.
187 This function computes du/dx as a weighted function of neighbouring points.
188 Here is implemented the linear model proposed in:
189 "Computing strain fields from discrete displacement fields in 2D-solids", Geers et al., 1996
191 Parameters
192 ----------
193 displacementField : 4D array of floats
194 The vector field to compute the derivatives.
195 Its shape is (nz, ny, nx, 3).
197 nodeSpacing : 3-component list of floats
198 Length between two nodes in every direction (*i.e.,* size of a cell)
200 neighbourRadius : float, optional
201 Distance in nodeSpacings to include neighbours in the strain calcuation.
202 Default = 1.5*nodeSpacing which will give radius = 1.5*min(nodeSpacing)
204 mask : bool, optional
205 Avoid non-correlated NaN points in the displacement field?
206 Default = True
208 nProcesses : integer, optional
209 Number of processes for multiprocessing
210 Default = number of CPUs in the system
212 verbose : boolean, optional
213 Print progress?
214 Default = True
216 Returns
217 -------
218 Ffield: nz x ny x nx x 3x3 array of n cells
219 The field of the transformation gradient tensors
221 Note
222 ----
223 Taken from the implementation in "TomoWarp2: A local digital volume correlation code", Tudisco et al., 2017
224 """
225 import scipy.spatial
227 # Define dimensions
228 dims = displacementField.shape[0:3]
229 nNodes = dims[0] * dims[1] * dims[2]
230 displacementFieldFlat = displacementField.reshape(nNodes, 3)
232 # Check if a 2D field is passed
233 twoD = False
234 if dims[0] == 1:
235 twoD = True
237 # Deformation gradient tensor F = du/dx +I
238 # Ffield = numpy.zeros((dims[0], dims[1], dims[2], 3, 3))
239 FfieldFlat = numpy.zeros((nNodes, 3, 3))
241 if twoD:
242 fieldCoordsFlat = (
243 numpy.mgrid[
244 0:1:1,
245 nodeSpacing[1] : dims[1] * nodeSpacing[1] + nodeSpacing[1] : nodeSpacing[1],
246 nodeSpacing[2] : dims[2] * nodeSpacing[2] + nodeSpacing[2] : nodeSpacing[2],
247 ]
248 .reshape(3, nNodes)
249 .T
250 )
251 else:
252 fieldCoordsFlat = (
253 numpy.mgrid[
254 nodeSpacing[0] : dims[0] * nodeSpacing[0] + nodeSpacing[0] : nodeSpacing[0],
255 nodeSpacing[1] : dims[1] * nodeSpacing[1] + nodeSpacing[1] : nodeSpacing[1],
256 nodeSpacing[2] : dims[2] * nodeSpacing[2] + nodeSpacing[2] : nodeSpacing[2],
257 ]
258 .reshape(3, nNodes)
259 .T
260 )
262 # Get non-nan displacements
263 goodPointsMask = numpy.isfinite(displacementField[:, :, :, 0].reshape(nNodes))
264 badPointsMask = numpy.isnan(displacementField[:, :, :, 0].reshape(nNodes))
265 # Flattened variables
266 fieldCoordsFlatGood = fieldCoordsFlat[goodPointsMask]
267 displacementFieldFlatGood = displacementFieldFlat[goodPointsMask]
268 # set bad points to nan
269 FfieldFlat[badPointsMask] = numpy.eye(3) * numpy.nan
271 # build KD-tree for neighbour identification
272 treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
274 # Output array for good points
275 FfieldFlatGood = numpy.zeros_like(FfieldFlat[goodPointsMask])
277 # Function for parallel mode
278 global _multiprocessingGeersOnePoint
280 def _multiprocessingGeersOnePoint(goodPoint):
281 # This is for the linear model, equation 15 in Geers
282 centralNodePosition = fieldCoordsFlatGood[goodPoint]
283 centralNodeDisplacement = displacementFieldFlatGood[goodPoint]
284 sX0X0 = numpy.zeros((3, 3))
285 sX0Xt = numpy.zeros((3, 3))
286 m0 = numpy.zeros(3)
287 mt = numpy.zeros(3)
289 # Option 2: KDTree on distance
290 # KD-tree will always give the current point as zero-distance
291 ind = treeCoord.query_ball_point(centralNodePosition, neighbourRadius * max(nodeSpacing))
293 # We know that the current point will also be included, so remove it from the index list.
294 ind = numpy.array(ind)
295 ind = ind[ind != goodPoint]
296 nNeighbours = len(ind)
297 nodalRelativePositionsRef = numpy.zeros((nNeighbours, 3)) # Delta_X_0 in paper
298 nodalRelativePositionsDef = numpy.zeros((nNeighbours, 3)) # Delta_X_t in paper
300 for neighbour, i in enumerate(ind):
301 # Relative position in reference configuration
302 # absolute position of this neighbour node
303 # minus abs position of central node
304 nodalRelativePositionsRef[neighbour, :] = fieldCoordsFlatGood[i] - centralNodePosition
306 # Relative position in deformed configuration (i.e., plus displacements)
307 # absolute position of this neighbour node
308 # plus displacement of this neighbour node
309 # minus abs position of central node
310 # minus displacement of central node
311 nodalRelativePositionsDef[neighbour, :] = fieldCoordsFlatGood[i] + displacementFieldFlatGood[i] - centralNodePosition - centralNodeDisplacement
313 for u in range(3):
314 for v in range(3):
315 # sX0X0[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
316 # sX0Xt[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
317 # Proposed solution for #142 for direction of rotation
318 sX0X0[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
319 sX0Xt[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
321 m0 += nodalRelativePositionsRef[neighbour, :]
322 mt += nodalRelativePositionsDef[neighbour, :]
324 sX0X0 = nNeighbours * sX0X0
325 sX0Xt = nNeighbours * sX0Xt
327 A = sX0X0 - numpy.dot(m0, m0)
328 C = sX0Xt - numpy.dot(m0, mt)
329 F = numpy.zeros((3, 3))
331 if twoD:
332 A = A[1:, 1:]
333 C = C[1:, 1:]
334 try:
335 F[1:, 1:] = numpy.dot(numpy.linalg.inv(A), C)
336 F[0, 0] = 1.0
337 except numpy.linalg.LinAlgError:
338 # print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
339 pass
340 else:
341 try:
342 F = numpy.dot(numpy.linalg.inv(A), C)
343 except numpy.linalg.LinAlgError:
344 # print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
345 pass
347 return goodPoint, F
349 # Iterate through flat field of Fs
350 if verbose:
351 pbar = progressbar.ProgressBar(maxval=fieldCoordsFlatGood.shape[0]).start()
352 finishedPoints = 0
354 # Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat
355 with multiprocessing.Pool(processes=nProcesses) as pool:
356 for returns in pool.imap_unordered(_multiprocessingGeersOnePoint, range(fieldCoordsFlatGood.shape[0])):
357 if verbose:
358 finishedPoints += 1
359 pbar.update(finishedPoints)
360 FfieldFlatGood[returns[0]] = returns[1]
362 if verbose:
363 pbar.finish()
365 FfieldFlat[goodPointsMask] = FfieldFlatGood
366 return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
369def FfieldBagi(points, connectivity, displacements, nProcesses=nProcessesDefault, verbose=False):
370 """
371 Calculates transformation gradient function using Bagi's 1996 paper, especially equation 3 on page 174.
372 Required inputs are connectivity matrix for tetrahedra (for example from spam.mesh.triangulate) and
373 nodal positions in reference and deformed configurations.
375 Parameters
376 ----------
377 points : m x 3 numpy array
378 M Particles' points in reference configuration
380 connectivity : n x 4 numpy array
381 Delaunay triangulation connectivity generated by spam.mesh.triangulate for example
383 displacements : m x 3 numpy array
384 M Particles' displacement
386 nProcesses : integer, optional
387 Number of processes for multiprocessing
388 Default = number of CPUs in the system
390 verbose : boolean, optional
391 Print progress?
392 Default = True
394 Returns
395 -------
396 Ffield: nx3x3 array of n cells
397 The field of the transformation gradient tensors
398 """
399 import spam.mesh
401 Ffield = numpy.zeros([connectivity.shape[0], 3, 3], dtype="<f4")
403 connectivity = connectivity.astype(numpy.uint)
405 # define dimension
406 # D = 3.0
408 # Import modules
410 # Construct 4-list of 3-lists of combinations constituting a face of the tet
411 combs = [[0, 1, 2], [1, 2, 3], [2, 3, 0], [0, 1, 3]]
412 unode = [3, 0, 1, 2]
414 # Precompute tetrahedron Volumes
415 tetVolumes = spam.mesh.tetVolumes(points, connectivity)
417 # Initialize arrays for tet strains
418 # print("spam.mesh.bagiStrain(): Constructing strain from Delaunay and Displacements")
420 # Loop through tetrahdra to get avec1, uPos1
421 global _multiprocessingComputeOneTet
423 def _multiprocessingComputeOneTet(tet):
424 # Get the list of IDs, centroids, center of tet
425 tetIDs = connectivity[tet, :]
426 # 2019-10-07 EA: Skip references to missing particles
427 # if max(tetIDs) >= points.shape[0]:
428 # print("spam.mesh.unstructured.bagiStrain(): this tet has node > points.shape[0], skipping")
429 # pass
430 # else:
431 if True:
432 tetCoords = points[tetIDs, :]
433 tetDisp = displacements[tetIDs, :]
434 tetCen = numpy.average(tetCoords, axis=0)
435 if numpy.isfinite(tetCoords).sum() + numpy.isfinite(tetDisp).sum() != 3 * 4 * 2:
436 if verbose:
437 print("spam.mesh.unstructured.bagiStrain(): nans in position or displacement, skipping")
438 # Compute strains
439 F = numpy.zeros((3, 3)) * numpy.nan
440 else:
441 # Loop through each face of tet to get avec, upos (Bagi, 1996, pg. 172)
442 # aVec1 = numpy.zeros([4, 3], dtype='<f4')
443 # uPos1 = numpy.zeros([4, 3], dtype='<f4')
444 # uPos2 = numpy.zeros([4, 3], dtype='<f4')
445 dudx = numpy.zeros((3, 3), dtype="<f4")
447 for face in range(4):
448 faceNorm = numpy.cross(
449 tetCoords[combs[face][0]] - tetCoords[combs[face][1]],
450 tetCoords[combs[face][0]] - tetCoords[combs[face][2]],
451 )
453 # Get a norm vector to face point towards center of tet
454 faceCen = numpy.average(tetCoords[combs[face]], axis=0)
455 tmpnorm = faceNorm / (numpy.linalg.norm(faceNorm))
456 facetocen = tetCen - faceCen
457 if numpy.dot(facetocen, tmpnorm) < 0:
458 tmpnorm = -tmpnorm
460 # Divide by 6 (1/3 for 1/Dimension; 1/2 for area from cross product)
461 # See first eqn., Bagi, 1996, pg. 172.
462 # aVec1[face] = tmpnorm*numpy.linalg.norm(faceNorm)/6
464 # Undeformed positions
465 # uPos1[face] = tetCoords[unode[face]]
466 # Deformed positions
467 # uPos2[face] = tetComs2[unode[face]]
469 dudx += numpy.tensordot(
470 tetDisp[unode[face]],
471 tmpnorm * numpy.linalg.norm(faceNorm) / 6,
472 axes=0,
473 )
475 dudx /= float(tetVolumes[tet])
477 F = numpy.eye(3) + dudx
478 return tet, F
480 # Iterate through flat field of Fs
481 if verbose:
482 pbar = progressbar.ProgressBar(maxval=connectivity.shape[0]).start()
483 finishedTets = 0
485 # Run multiprocessing
486 with multiprocessing.Pool(processes=nProcesses) as pool:
487 for returns in pool.imap_unordered(_multiprocessingComputeOneTet, range(connectivity.shape[0])):
488 if verbose:
489 finishedTets += 1
490 pbar.update(finishedTets)
491 Ffield[returns[0]] = returns[1]
492 pool.close()
493 pool.join()
495 if verbose:
496 pbar.finish()
498 return Ffield
501def decomposeFfield(Ffield, components, twoD=False, nProcesses=nProcessesDefault, verbose=False):
502 """
503 This function takes in an F field (from either FfieldRegularQ8, FfieldRegularGeers, FfieldBagi) and
504 returns fields of desired transformation components.
506 Parameters
507 ----------
508 Ffield : multidimensional x 3 x 3 numpy array of floats
509 Spatial field of Fs
511 components : list of strings
512 These indicate the desired components consistent with spam.deformation.decomposeF or decomposePhi
514 twoD : bool, optional
515 Is the Ffield in 2D? This changes the strain calculation.
516 Default = False
518 nProcesses : integer, optional
519 Number of processes for multiprocessing
520 Default = number of CPUs in the system
522 verbose : boolean, optional
523 Print progress?
524 Default = True
526 Returns
527 -------
528 Dictionary containing appropriately reshaped fields of the transformation components requested.
530 Keys:
531 vol, dev, volss, devss are scalars
532 r, z, Up are 3-component vectors
533 e and U are 3x3 tensors
534 """
535 # The last two are for the 3x3 F field
536 fieldDimensions = Ffield.shape[0:-2]
537 fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
538 FfieldFlat = Ffield.reshape(fieldRavelLength, 3, 3)
540 output = {}
541 for component in components:
542 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
543 output[component] = numpy.zeros(fieldRavelLength)
544 if component == "r" or component == "z" or component == "Up":
545 output[component] = numpy.zeros((fieldRavelLength, 3))
546 if component == "U" or component == "e":
547 output[component] = numpy.zeros((fieldRavelLength, 3, 3))
549 # Function for parallel mode
550 global _multiprocessingDecomposeOneF
552 def _multiprocessingDecomposeOneF(n):
553 F = FfieldFlat[n]
554 if numpy.isfinite(F).sum() == 9:
555 decomposedF = spam.deformation.decomposeF(F, twoD=twoD)
556 return n, decomposedF
557 else:
558 return n, {
559 "r": numpy.array([numpy.nan] * 3),
560 "z": numpy.array([numpy.nan] * 3),
561 "Up": numpy.array([numpy.nan] * 3),
562 "U": numpy.eye(3) * numpy.nan,
563 "e": numpy.eye(3) * numpy.nan,
564 "vol": numpy.nan,
565 "dev": numpy.nan,
566 "volss": numpy.nan,
567 "devss": numpy.nan,
568 }
570 # Iterate through flat field of Fs
571 if verbose:
572 pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
573 finishedPoints = 0
575 # Run multiprocessing
576 with multiprocessing.Pool(processes=nProcesses) as pool:
577 for returns in pool.imap_unordered(_multiprocessingDecomposeOneF, range(fieldRavelLength)):
578 if verbose:
579 finishedPoints += 1
580 pbar.update(finishedPoints)
581 for component in components:
582 output[component][returns[0]] = returns[1][component]
583 pool.close()
584 pool.join()
586 if verbose:
587 pbar.finish()
589 # Reshape on the output
590 for component in components:
591 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
592 output[component] = numpy.array(output[component]).reshape(fieldDimensions)
594 if component == "r" or component == "z" or component == "Up":
595 output[component] = numpy.array(output[component]).reshape(Ffield.shape[0:-1])
597 if component == "U" or component == "e":
598 output[component] = numpy.array(output[component]).reshape(Ffield.shape)
600 return output
603def decomposePhiField(PhiField, components, twoD=False, nProcesses=nProcessesDefault, verbose=False):
604 """
605 This function takes in a Phi field (from readCorrelationTSV?) and
606 returns fields of desired transformation components.
608 Parameters
609 ----------
610 PhiField : multidimensional x 4 x 4 numpy array of floats
611 Spatial field of Phis
613 components : list of strings
614 These indicate the desired components consistent with spam.deformation.decomposePhi
616 twoD : bool, optional
617 Is the PhiField in 2D? This changes the strain calculation.
618 Default = False
620 nProcesses : integer, optional
621 Number of processes for multiprocessing
622 Default = number of CPUs in the system
624 verbose : boolean, optional
625 Print progress?
626 Default = True
628 Returns
629 -------
630 Dictionary containing appropriately reshaped fields of the transformation components requested.
632 Keys:
633 vol, dev, volss, devss are scalars
634 t, r, z and Up are 3-component vectors
635 e and U are 3x3 tensors
636 """
637 # The last two are for the 4x4 Phi field
638 fieldDimensions = PhiField.shape[0:-2]
639 fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
640 PhiFieldFlat = PhiField.reshape(fieldRavelLength, 4, 4)
642 output = {}
643 for component in components:
644 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
645 output[component] = numpy.zeros(fieldRavelLength)
646 if component == "t" or component == "r" or component == "z" or component == "Up":
647 output[component] = numpy.zeros((fieldRavelLength, 3))
648 if component == "U" or component == "e":
649 output[component] = numpy.zeros((fieldRavelLength, 3, 3))
651 # Function for parallel mode
652 global _multiprocessingDecomposeOnePhi
654 def _multiprocessingDecomposeOnePhi(n):
655 Phi = PhiFieldFlat[n]
656 if numpy.isfinite(Phi).sum() == 16:
657 decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
658 return n, decomposedPhi
659 else:
660 return n, {
661 "t": numpy.array([numpy.nan] * 3),
662 "r": numpy.array([numpy.nan] * 3),
663 "z": numpy.array([numpy.nan] * 3),
664 "Up": numpy.array([numpy.nan] * 3),
665 "U": numpy.eye(3) * numpy.nan,
666 "e": numpy.eye(3) * numpy.nan,
667 "vol": numpy.nan,
668 "dev": numpy.nan,
669 "volss": numpy.nan,
670 "devss": numpy.nan,
671 }
673 # Iterate through flat field of Fs
674 if verbose:
675 pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
676 finishedPoints = 0
678 # Run multiprocessing
679 with multiprocessing.Pool(processes=nProcesses) as pool:
680 for returns in pool.imap_unordered(_multiprocessingDecomposeOnePhi, range(fieldRavelLength)):
681 if verbose:
682 finishedPoints += 1
683 pbar.update(finishedPoints)
684 for component in components:
685 output[component][returns[0]] = returns[1][component]
686 pool.close()
687 pool.join()
689 if verbose:
690 pbar.finish()
692 # Reshape on the output
693 for component in components:
694 if component == "vol" or component == "dev" or component == "volss" or component == "devss":
695 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2])
697 if component == "t" or component == "r" or component == "z":
698 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3)
700 if component == "U" or component == "e":
701 output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3, 3)
703 return output