Coverage for /usr/local/lib/python3.10/site-packages/spam/DIC/globalDVC.py: 40%

1# Library of SPAM image correlation functions. 

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/>. 

16 

17""" 

18Set of functions for performing mechanically regularised global DVC. 

19The regularisation implementation is taken from [Mendoza2019]_. 

20 

21Here is how it can be used: 

22 

23.. code-block:: python 

24 

25 import spam.mesh 

26 import spam.DIC 

27 

28 # Mesh 

29 points, connectivity = spam.mesh.createCylinder([-24, 27], 12.5, 97.3, 10, zOrigin=-4.2) 

30 

31 # Regularisation: parameters 

32 myParameters = { 

33 "young": 10, 

34 "poisson": 0.25, 

35 "dirichlet": { 

36 "z_start": {"dof": [0]}, 

37 "z_end": {"dof": [0]}, 

38 }, 

39 } 

40 p = spam.DIC.regularisationParameters(myParameters) 

41 

42 # Regularisation step 1: get labels 

43 labels = spam.DIC.surfaceLabels(points, surfaces=p["surfaces"]) 

44 

45 # Regularisation step 2: build regularisation matrix 

46 regularisationMatrix, regularisationField = spam.DIC.regularisationMatrix( 

47 points, 

48 connectivity, 

49 young=p["young"], 

50 poisson=p["poisson"], 

51 xiBulk=p["xi"], 

52 dirichlet=p["dirichlet"], 

53 periods=p["periods"], 

54 labels=labels, 

55 ) 

56 

57 # Global DVC 

58 globalCorrelation( 

59 im1, 

60 im2, 

61 points, 

62 connectivity, 

63 regularisationMatrix=regularisationMatrix, 

64 regularisationField=regularisationField, 

65 ) 

66 

67.. [Mendoza2019] A. Mendoza, J. Neggers, F. Hild, S Roux (2019). Complete mechanical regularization applied to digital image and volume correlation. 

68 *Computer Methods in Applied Mechanics and Engineering*, Volume 355, 1 October 2019, Pages 27-43 

69 https://doi.org/10.1016/j.cma.2019.06.005 

70 

71.. _RegularisationParameter: _static/gdvc-regularisation-parameters.yaml 

72 

73 Note 

74 ---- 

75 Make a link to the script 

76 

77""" 

78 

79import time # for debug 

80 

81import numpy 

82import spam.DIC 

83import spam.label 

84import spam.mesh 

85import tifffile 

86 

87# 2017-05-29 ER and EA 

88# This is spam's C++ DIC toolkit, but since we're in the tools/ directory we can import it directly 

89from spambind.DIC.DICToolkit import ( 

90 computeDICglobalMatrix, 

91 computeDICglobalVector, 

92 computeGradientPerTet, 

93) 

94 

95 

96def _check_symmetric(a, rtol=1e-05, atol=1e-08): 

97 """Helper function to check if a matrix is symmetric.""" 

98 return numpy.allclose(a, a.T, rtol=rtol, atol=atol) 

99 

100 

101def _errorCalc(im1, im2, im2ref, meshPaddingSlice): 

102 """Helper function to compute the error between two images.""" 

103 errorInitial = numpy.sqrt(numpy.square(im2ref[meshPaddingSlice] - im1[meshPaddingSlice]).sum()) 

104 errorCurrent = numpy.sqrt(numpy.square(im2[meshPaddingSlice] - im1[meshPaddingSlice]).sum()) 

105 return errorCurrent / errorInitial 

106 

107 

108def _computeFunctional(u, K, L=None): 

109 """Helper function to compute global DVC functional""" 

110 u = numpy.ravel(u) 

111 if L is None: 

112 return numpy.matmul(u.T, numpy.matmul(K.T, numpy.matmul(K, u))) 

113 else: 

114 return numpy.matmul(u.T, numpy.matmul(K.T, numpy.matmul(L, numpy.matmul(K, u)))) 

115 

116 

117def _normalisedEnergies(v, Km, Ks, Ls): 

118 """Helper function to compute globale DVC normalised energies""" 

119 Em = _computeFunctional(v, Km) 

120 Es = [_computeFunctional(v, K, L=L) for K, L in zip(Ks, Ls)] 

121 return Em, Es 

122 

123 

124def _computeWeights(kMag, xiBulk, xiDirichlet): 

125 """Helper function to compute global DVC weights""" 

126 print(f"[regularisation] xi bulk = {xiBulk:.2f}") 

127 Wm = (xiBulk * kMag) ** 4 

128 Ws = [] 

129 for i, xi in enumerate(xiDirichlet): 

130 print(f"[regularisation] xi dirichlet {i} = {xi:.2f}") 

131 Ws.append((xi * kMag) ** 4) 

132 

133 return Wm, Ws 

134 

135 

136def surfaceLabels(points, surfaces=[], connectivity=None): 

137 """Creates a label for each points based on a list of keywords that corresponds to surfaces (`ie.` ``["z_start", "z_end"]``). 

138 The label value is based on the position of the surface in the list. 

139 

140 Parameters 

141 ---------- 

142 points: Nx3 array 

143 List of coordinates of the mesh nodes. 

144 

145 surfaces: list of string 

146 A list of keywords corresponding to surfaces. 

147 

148 - ``z_start``: plane at ``z == min(points[:,0])`` 

149 - ``z_end``: plane at ``z == max(points[:,0])`` 

150 - ``y_start``: plane at ``y == min(points[:,1])`` 

151 - ``y_end``: plane at ``y == max(points[:,1])`` 

152 - ``x_start``: plane at ``x == min(points[:,2])`` 

153 - ``x_end``: plane at ``x == max(points[:,2])`` 

154 - ``z_lateral``: lateral surface of a cylinder oriented in the first direction. 

155 - ``y_lateral``: lateral surface of a cylinder oriented in the second direction. 

156 - ``x_lateral``: lateral surface of a cylinder oriented in the third direction. 

157 

158 connectivity: array (only for debug purposes) 

159 Connectivity matrix. If set, creates a VTK file with labels. 

160 

161 Returns 

162 ------- 

163 N x 1 array: 

164 Surface label for each points. 

165 

166 Note 

167 ---- 

168 - Surface labels start at 1, 0 corresponding to bulk or not specified surfaces. 

169 - Points belong to a single surface. The first surface in `surfaces` prevails. 

170 

171 Example 

172 ------- 

173 >>> import spam.mesh 

174 >>> import spam.DIC 

175 >>> 

176 >>> # create mesh 

177 >>> points, connectivity = spam.mesh.createCylinder([-24, 27], 12.5, 97.3, 10, zOrigin=-4.2) 

178 >>> # compute labels for bottom and top surface only 

179 >>> labels = spam.DIC.surfaceLabels(points, surfaces=["z_start", "z_end"], connectivity=connectivity) 

180 

181 """ 

182 

183 def _belongs_to_lateral_surface(point, centre, radius, epsilon=1e-6): 

184 """Returns True if point belongs to the lateral surface of a cylinder""" 

185 d = ((centre[0] - point[1]) ** 2 + (centre[1] - point[2]) ** 2) ** 0.5 

186 return abs(d - radius) < epsilon 

187 

188 def _belongs_to_plane(point, direction, coordinate, epsilon=1e-6): 

189 """Returns True if point belongs to a surface of position `coordinate` in direction `direction`""" 

190 return abs(point[direction] - coordinate) < epsilon 

191 

192 # Get geometrical data from the points coordinates 

193 maxCoord = numpy.max(points, axis=0) 

194 minCoord = numpy.min(points, axis=0) 

195 centre = [0.0] * 3 

196 radii = [0.0] * 3 

197 for dz in range(3): 

198 dy = (dz + 1) % 3 

199 dx = (dz + 2) % 3 

200 centre[dz] = [0.5 * (minCoord[d] + maxCoord[d]) for d in [dy, dx]] 

201 radii[dz] = 0.25 * (maxCoord[dx] + maxCoord[dy] - minCoord[dx] - minCoord[dy]) 

202 

203 labels = numpy.zeros(points.shape[0], dtype=int) 

204 

205 # Loop over the points 

206 for A, point in enumerate(points): 

207 

208 # Loop over the surfaces to enforce order for edges and vertices 

209 for i, surface in enumerate(surfaces): 

210 

211 direction, position = surface.split("_") 

212 direction = {"z": 0, "y": 1, "x": 2}[direction] # direction as integer z: 0, y: 1, x: 2 

213 

214 if position == "start" and _belongs_to_plane(point, direction, minCoord[direction]): 

215 labels[A] = i + 1 

216 break 

217 

218 elif position == "end" and _belongs_to_plane(point, direction, maxCoord[direction]): 

219 labels[A] = i + 1 

220 break 

221 

222 elif position == "lateral" and _belongs_to_lateral_surface(point, centre[direction], radii[direction]): 

223 labels[A] = i + 1 

224 break 

225 

226 if connectivity is not None: # debug 

227 spam.helpers.writeUnstructuredVTK(points, connectivity, pointData={"label": labels}, fileName="gdvc-labels.vtk") 

228 

229 return labels 

230 

231 

232def projectionMatrices(points, connectivity, labels, dirichlet=[]): 

233 """Compute binary diagonal matrices and Laplace-Beltrami operator. 

234 Details on the meaning of these matrices can be found in [Mendoza2019]_ `eq. (12) to (14)` and `eq. (17) to (19)`. 

235 

236 

237 Parameters 

238 ---------- 

239 points: Nx3 array 

240 List of coordinates of the mesh nodes. 

241 

242 connectivity: Mx4 numpay array 

243 Connectivity matrice of the mesh (tetrahedra). 

244 

245 labels: Nx1 array 

246 Surface labels for each points as defined in ``spam.DIC.surfaceLabels()`` 

247 

248 dirichlet: list of (int, list, ) 

249 Each element of the list defines a surface with Dirichlet boundary conditions by a tuple. 

250 

251 - The first element of the tuple is the surface label which should belong to `labels`. 

252 - The second element of the tuple is a list of degrees of freedom directions to consider for the Dirichlet boundary conditions. 

253 - Extra elements in the tuple are ignored. 

254 .. code-block:: python 

255 

256 dirichlet = [ 

257 # surface 1, dof in x, y 

258 ( 

259 1, 

260 [1, 2], 

261 ), 

262 # surface 3, dof in z 

263 ( 

264 3, 

265 [0], 

266 ), 

267 ] 

268 

269 Returns 

270 ------- 

271 3Nx3N array: 

272 :math:`D_m` The binary diagonal matrix corresponding to all dofs of the bulk and neumann surfaces nodes (ie the non dirichlet) 

273 List of 3Nx3N array: 

274 :math:`D_{S_i}` Binary diagonal matrices corresponding to the dofs of each Dirichlet surfaces 

275 List of 3Nx3N array: 

276 :math:`L_{S_i}` List of Laplace-Beltrami operators corresponding to the same Dirichlet surfaces. 

277 

278 """ 

279 if dirichlet is None: 

280 dirichlet = [] 

281 

282 # all matrices out here are of shape ndof x ndof (same as stiffness matrix) 

283 ndof = 3 * points.shape[0] 

284 

285 print(f"[projectionMatrices] number of dirichlet surfaces: {len(dirichlet)} ({dirichlet})") 

286 

287 # initialise list of A and B matrices in order to compute the L operator (size ndof x ndof) 

288 matricesA = [numpy.zeros((ndof, ndof)) for _ in dirichlet] 

289 matricesB = [numpy.zeros((ndof, ndof)) for _ in dirichlet] 

290 matricesD = [numpy.zeros((ndof, ndof)) for _ in dirichlet] 

291 

292 # list of all nodes and dof numbers on dirichlet surfaces 

293 # NOTE: taking label == 0 (ie background) leads to singlular sub_A matrix 

294 dirichletSurfaces = [(numpy.where(labels == d[0])[0], d[1]) for d in dirichlet if d[0]] 

295 

296 DEBUG_dirichlet_connectivity_all = [] 

297 DEBUG_dirichlet_points_all = [] 

298 for si, _ in enumerate(dirichlet): 

299 DEBUG_dirichlet_connectivity_all.append([]) 

300 DEBUG_dirichlet_points_all.append([]) 

301 

302 # # loop over dirichlet matrix 

303 for si, (surfacePoints, dofDirections) in enumerate(dirichletSurfaces): 

304 dofs = [3 * A + d for A in surfacePoints for d in dofDirections] 

305 

306 print(f"[projectionMatrices] Dirichlet #{si} label {dirichlet[si][0]}) #points = {len(surfacePoints)} dofs directions = {dofDirections} {len(dofs)} dofs") 

307 

308 surfacePoints = set(surfacePoints) 

309 for tetra_connectivity in connectivity: 

310 tri_con = list(surfacePoints.intersection(tetra_connectivity)) 

311 tri_nodes = [list(points[i]) for i in tri_con] 

312 if len(tri_con) != 3: 

313 # the element doesn't have a triangle a this surface 

314 continue 

315 

316 # get 4th point of the tetrahedron to compute 3D shape functions 

317 point_4_n = list(set(tri_con).symmetric_difference(set(tetra_connectivity)))[0] 

318 _, tri_coefficients = spam.mesh.shapeFunctions(tri_nodes + [list(points[point_4_n])]) 

319 # we've got a triangle on the surface!!! 

320 DEBUG_dirichlet_connectivity_all[si].append(list(tri_con)) 

321 DEBUG_dirichlet_points_all[si].append(tri_nodes) 

322 

323 # assemble the dirichlet connectivity matrix 

324 a = numpy.array(tri_coefficients[:3, 0]) 

325 bcd = numpy.array(tri_coefficients[:3, 1:]) 

326 B = numpy.array(tri_nodes).T 

327 

328 def phi(i, L): 

329 return a[i] + numpy.matmul(bcd[i], numpy.matmul(B, L.T)) 

330 

331 def dphi(i): 

332 return bcd[i] 

333 

334 # gauss points 

335 L_gp = numpy.array([[0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]]) 

336 

337 # STEP 3: compute the area 

338 area = 0.5 * numpy.linalg.norm( 

339 numpy.cross( 

340 numpy.subtract(tri_nodes[1], tri_nodes[0]), 

341 numpy.subtract(tri_nodes[2], tri_nodes[0]), 

342 ) 

343 ) 

344 

345 # STEP 3: compute inner products of the shape functions 

346 for i in range(3): 

347 for j in range(3): 

348 inner = 0.0 

349 for L in L_gp: 

350 # print(inner, i, j, L, phi(i, L), phi(j, L)) 

351 inner += phi(i, L) * phi(j, L) 

352 inner *= area / 3.0 # the 1/3. comes from the weight 

353 dinner = area * numpy.inner(dphi(i), dphi(j)) 

354 for d in dofDirections: 

355 P = int(3 * tri_con[i] + d) 

356 Q = int(3 * tri_con[j] + d) 

357 matricesA[si][P, Q] += inner 

358 matricesB[si][P, Q] += dinner 

359 

360 # # invert matrix A 

361 # for si, (surfacePoints, dofDirections) in enumerate(dirichletSurfaces): 

362 dofs = [3 * A + d for A in surfacePoints for d in dofDirections] 

363 # 1. extract submatrix A from full size matricesA and invert 

364 sub_A = matricesA[si][:, dofs] 

365 sub_A = sub_A[dofs, :] 

366 

367 sub_A = numpy.linalg.inv(sub_A) 

368 # 2. push back inverted submatrix into full size matricesA 

369 for i, P in enumerate(dofs): 

370 matricesD[si][P, P] = 1 

371 for j, Q in enumerate(dofs): 

372 matricesA[si][P, Q] = sub_A[i, j] 

373 

374 # return also bulk + neumann binary diagonal projection matrices 

375 # Dm = Db + Dn = I - Sum(Dd) 

376 Dm = numpy.eye(ndof) 

377 for Ds in matricesD: 

378 Dm -= Ds 

379 

380 # return a list of D and L = AxB for each dirichlet surface 

381 return Dm, matricesD, [numpy.matmul(A, B) for A, B in zip(matricesA, matricesB)] 

382 

383 

384def isochoricField(points, periods=3, connectivity=None): 

385 r"""Helper function to build the isochoric test function used to normalised the mechanical regularisation functionals. 

386 The function is a shear displacement field with its wave vector in the direction of the longest dimension of the mesh. 

387 

388 .. math:: 

389 

390 \textbf{v}(\textbf{x}) = \sin(2 \pi \textbf{k}\cdot \textbf{x} ) 

391 

392 [Mendoza2019]_ `eq. (25)`. 

393 

394 Parameters 

395 ---------- 

396 points: Nx3 array 

397 List of coordinates of the mesh nodes. 

398 

399 periods: int 

400 The number of periods of the shear wave. 

401 

402 connectivity: array (only for debug purposes) 

403 Connectivity matrix. If set, creates a VTK file with the field. 

404 

405 Returns 

406 ------- 

407 float: 

408 The magnitude of the wave vector :math:`\bf{k}`. 

409 Nx3 array: 

410 The field :math:`\bf{v}` at each nodes. 

411 """ 

412 

413 sizes = [ma - mi for ma, mi in zip(numpy.max(points, axis=0), numpy.min(points, axis=0))] 

414 d = numpy.argmax(sizes) # direction of largest size 

415 size = float(max(sizes)) # largest size 

416 mag = periods / size 

417 v = numpy.zeros_like(points) 

418 for i, point in enumerate(points): 

419 kdotx = mag * point[d] 

420 v[i][(d + 0) % 3] = 0 

421 v[i][(d + 1) % 3] = numpy.sqrt(0.5) * numpy.sin(2 * numpy.pi * kdotx) 

422 v[i][(d + 2) % 3] = numpy.sqrt(0.5) * numpy.sin(2 * numpy.pi * kdotx) 

423 

424 if connectivity is not None: 

425 print("[isochoricField] plot vtk") 

426 spam.helpers.writeUnstructuredVTK( 

427 points, 

428 connectivity, 

429 elementType="tetra", 

430 pointData={"v": v}, 

431 cellData={}, 

432 fileName="gdvc-v.vtk", 

433 ) 

434 

435 return mag, v 

436 

437 

438def regularisationMatrix(points, connectivity, young=1, poisson=0.25, voxelSize=1, xiBulk=None, dirichlet=[], labels=[], periods=3): 

439 r"""Computes the mechanical regularisation matrix :math:`M_{reg}` for the global DVC: 

440 $$(M_{reg} + M_{c}) \\delta\\textbf{u} = \\textbf{b} - M_{reg} \\textbf{u} $$ 

441 where 

442 $$M_{reg} = M_m + \\sum_{i} M_{S_i}$$ 

443 corresponds to both bulk/Neuman and Dirichlet surfaces contributions ([Mendoza2019]_ `eq. (29) and (30)`). 

444 

445 Parameters 

446 ---------- 

447 points: Nx3 array 

448 List of coordinates of the mesh nodes. 

449 

450 connectivity: Mx4 array 

451 Connectivity matrix of the mesh. 

452 

453 young: float (optional) 

454 The Young modulus used to build the stiffness matrix :math:`K` from which :math:`M_m` and :math:`M_{S_i}` are derived. 

455 Default = 1 (it is useless to change this value if you don't impose forces with meaningful units on the Neumann surfaces) 

456 

457 poisson: float (optional) 

458 The Poisson ratio used to build the stiffness matrix :math:`K` from which :math:`M_m` and :math:`M_{S_i}` are derived. 

459 Default = 0.25 

460 

461 voxelSize: float (optional) 

462 The size of a voxel in the same unit used to set the Young modulus assuming that the mesh geometry is in voxels. 

463 This voxel size :math:`v_s` is used to convert the Young modulus in Newton per voxels squared. 

464 Default = 1 (it is useless to change this value if you don't impose forces with meaningful units on the Neumann surfaces). 

465 

466 xiBulk: float (optional) 

467 The regularisation length :math:`\xi_m` of the bulk/Neumann contribution :math:`M_{reg}`. 

468 It has to be compared to the characteristic length of the mesh. 

469 The bigger the regularisation length the more important the regularisation contribution. 

470 If None it is taken to be equal to the mesh characteristic length. 

471 Default = None 

472 

473 dirichlet: list of (int, list, float) (optional) 

474 Each element of the list defines a surface with Dirichlet boundary conditions by a tuple. 

475 

476 - The first element of the tuple is the surface label which should belong to `labels`. 

477 - The second element of the tuple is a list of degrees of freedom directions to consider for the Dirichlet boundary conditions. 

478 - The third element of the tuple correspond to the regularisation length of the surface :math:`\xi_{S_i}` 

479 

480 .. code-block:: python 

481 

482 dirichlet = [ 

483 # surface 1, dof in x, y, xi = 0.1 

484 [1, [1, 2], 0.1], 

485 # surface 3, dof in z, xi not set 

486 [3, [0], None], 

487 ] 

488 

489 labels: Nx1 array (optional) 

490 Surface labels for each points as defined in ``spam.DIC.surfaceLabels()``. Mandatory only if ``dirichlet`` is set. 

491 

492 periods: float (optional) 

493 Number of periods of the isochoric function :math:`v` used to compute the normalized energy ([Mendoza2019]_ `eq. (27)`). 

494 :math:`v` is computed with ``spam.DIC.isochoricField()`` with a default number of periods of 3. 

495 

496 Returns 

497 ------- 

498 3Nx3N array: 

499 :math:`M_{req}` the regularisation matrix. 

500 Nx3 array: 

501 :math:`v` the isochoric field at each nodes of the mesh. 

502 

503 Note 

504 ---- 

505 - The ``dirichlet`` argument is compatible with the one in ``spam.DIC.projectionMatrices()`` 

506 - As long as you don't impose forces to the Neumann surfaces it is usless to set a specific Young modulus and a voxel size. 

507 - Imposing forces is not implemented yet :) 

508 

509 Warning 

510 ------- 

511 This function is not tested yet 

512 """ 

513 

514 def _imshow(matlist, title): 

515 from matplotlib import pyplot as plt 

516 

517 if not isinstance(matlist, list): 

518 matlist = [matlist] 

519 

520 for i, mat in enumerate(matlist): 

521 plt.imshow(mat) 

522 plt.title(f"{title} {i}") 

523 # plt.show() 

524 

525 print(f"[regularisation] Young modulus using si: {young}") 

526 young *= voxelSize**2 

527 print(f"[regularisation] Young modulus using voxels: {young}") 

528 

529 print("[regularisation] Build projection matrices") 

530 Dm, Ds, Ls = projectionMatrices(points, connectivity, labels, dirichlet=dirichlet) 

531 _imshow(Dm, "Dm") 

532 _imshow(Ds, "Ds") 

533 _imshow(Ls, "Ls") 

534 

535 print("[regularisation] Build global stiffness matrix") 

536 K = spam.mesh.globalStiffnessMatrix(points, connectivity, young, poisson) 

537 print("[regularisation] Build bulk stiffness matrix") 

538 Km = numpy.matmul(Dm, K) 

539 print("[regularisation] Build dirichlet stiffness matrices") 

540 Ks = [numpy.matmul(Ds_i, K) for Ds_i in Ds] 

541 _imshow(K, "K") 

542 _imshow(Km, "Km") 

543 _imshow(Ks, "Ks") 

544 del K 

545 

546 print("[regularisation] Build isochoric function") 

547 kMag, v = isochoricField(points, periods=periods) 

548 

549 print("[regularisation] Compute normalised energy and weight") 

550 Em, Es = _normalisedEnergies(v, Km, Ks, Ls) 

551 print(f'[regularisation] Em = {Em:.3e}, Es = {" ".join([f"{_:.3e}" for _ in Es])}') 

552 

553 lc = spam.mesh.getMeshCharacteristicLength(points, connectivity) 

554 print(f"[regularisation] mesh characteristic length lc = {lc}") 

555 

556 xiBulk = lc if xiBulk is None else xiBulk 

557 xiDirichlet = [lc if _[2] is None else _[2] for _ in dirichlet] 

558 Wm, Ws = _computeWeights(kMag, xiBulk, xiDirichlet) 

559 print(f'[regularisation] Wm = {Wm:.2f}, Ws = {" ".join([f"{_:.2f}" for _ in Ws])}') 

560 

561 print(f'[regularisation] Wm/Em = {Wm/Em:.3e} Ws/Es = {" ".join([f"{a/b:.3e}" for a, b in zip(Ws, Es)])}') 

562 

563 # 5.4 compute Mreg 

564 Mreg = numpy.zeros_like(Km) 

565 Mreg = Wm * numpy.matmul(Km.T, Km) / Em 

566 for W, E, K, L in [(W, E, K, L) for W, E, K, L in zip(Ws, Es, Ks, Ls) if E]: 

567 Mreg += W * numpy.matmul(K.T, numpy.matmul(L, K)) / E 

568 

569 _imshow(Mreg, "Mreg") 

570 return Mreg, v 

571 

572 

573def regularisationParameters(userParameters): 

574 """ 

575 Convert a user friendly dictionary of parameters into a dictionary of variables compatible with the regularisation functions of this module. 

576 

577 Parameters 

578 ---------- 

579 userParameters: (str | dict) 

580 The user parameters for the mechanical regularisation scheme. It can be: 

581 

582 - if ``str``: the path to a YAML file. A dummy example can be downloaded here: `RegularisationParameter`_ 

583 - if ``dict``: a dictionary containing the following keys and values: 

584 

585 .. code-block:: python 

586 

587 { 

588 # bulk / Neumann regularisation 

589 "young": 10, # mandatory (the Young modulus) 

590 "poisson": 0.25, # mandatory (the Poisson ratio) 

591 "xiBulk": 30, # optional (the bulk/Neumann regularisation lentgh) 

592 "periods": 3, # the number of periods of the isochoric function (optional) 

593 "voxelSize": 0.01, # the voxel size (optional) 

594 # Information about the Dirichlet surface regularisation 

595 # Each surface is defined by a search keyword among 

596 # z_start, z_end, y_start, y_end, x_start and x_end for plane lookups 

597 # z_lateral, y_lateral, x_lateral for cylinder lateral surface lookups 

598 "dirichlet": { 

599 "z_start": { # surface label 1 

600 "xi": 5, # the regularisation length (optional) 

601 "dof": [0, 1, 2], # mandatory 

602 }, 

603 "z_end": {"dof": [0]}, # surface label 2 # mandatory 

604 }, 

605 } 

606 

607 Returns 

608 ------- 

609 dict: 

610 A dictionary with the variables needed for the regularisation functions. 

611 The dictionary keys are named to match the function's signatures: 

612 

613 - ``surface``: the dirichlet surfaces for ``spam.DIC.surfaceLabels()`` 

614 - ``dirichlet``: the dirichlet surfaces for ``spam.DIC.regularisationMatrix()`` 

615 - ``young``: the Young modulus for ``spam.DIC.regularisationMatrix()`` 

616 - ``poisson``: the Poisson ratio for ``spam.DIC.regularisationMatrix()`` 

617 - ``xiBulk``: the bulk characteristic lentgh for ``spam.DIC.regularisationMatrix()`` 

618 - ``periods``: the Poisson ratio for ``spam.DIC.regularisationMatrix()`` 

619 - ``voxelSize``: the Poisson ratio for ``spam.DIC.regularisationMatrix()`` 

620 """ 

621 

622 surfaces = [k for k in userParameters.get("dirichlet", {})] 

623 dirichlet = [(i + 1, surf["dof"], surf.get("xi")) for i, surf in enumerate(userParameters.get("dirichlet", {}).values())] 

624 young = userParameters["young"] 

625 poisson = userParameters["poisson"] 

626 xiBulk = userParameters.get("xiBulk") 

627 periods = userParameters.get("periods", 3) 

628 voxelSize = userParameters.get("voxelSize", 1.0) 

629 

630 parameters = {"surfaces": surfaces, "young": young, "poisson": poisson, "xiBulk": xiBulk, "dirichlet": dirichlet, "periods": periods, "voxelSize": voxelSize} 

631 

632 for k, v in parameters.items(): 

633 print(f"[regularisation parameters] {k:<10}: {v}") 

634 return parameters 

635 

636 

637def globalCorrelation( 

638 im1, 

639 im2, 

640 points, 

641 connectivity, 

642 regularisationMatrix=None, 

643 regularisationField=None, 

644 initialDisplacements=None, 

645 convergenceCriterion=0.01, 

646 maxIterations=20, 

647 medianFilterEachIteration=False, 

648 debugFiles=False, 

649 prefix="globalCorrelation", 

650 nThreads=None, 

651): 

652 """ 

653 Global DVC (works only with 3D images). 

654 

655 Parameters 

656 ---------- 

657 im1 : 3D array 

658 Reference image in which the mesh is defined 

659 

660 im2 : 3D array 

661 Deformed image, should be same size as im1 

662 

663 points : M x 3 array 

664 M nodal coordinates in reference configuration 

665 

666 connectivity : N x 4 array 

667 connectivityhedral connectivity generated by spam.mesh.triangulate() for example 

668 

669 regularisationMatrix : 3N x 3N array (optional) 

670 Mechanical regularisation stiffness matrix. If None no mechanical regularisation is applied. 

671 First output of `spam.DIC.regularisationMatrix` 

672 

673 

674 regularisationField : N x 3 array (optional) 

675 Isochoric displacement field used to compute the normalisation energies. 

676 Second output of `spam.DIC.regularisationMatrix` 

677 

678 initialDisplacements : M x 3 array of floats (optional) 

679 Initial guess for nodal displacements, must be coherent with input mesh 

680 Default = None 

681 

682 convergenceCriterion : float 

683 Convergence criterion for change in displacements in px 

684 Default = 0.01 

685 

686 maxIterations : int 

687 Number of iterations to stop after if convergence has not been reached 

688 Default = 20 

689 

690 debugFiles : bool 

691 Write temporary results to file for debugging? 

692 Default = 'globalCorrelation' 

693 

694 prefix : string 

695 Output file prefix for debugFiles 

696 Default = None 

697 

698 Returns 

699 ------- 

700 displacements : N x 3 array of floats 

701 (converged?) Nodal displacements 

702 

703 Example 

704 ------- 

705 >>> import spam.DIC 

706 >>> spam.DIC.globalCorrelation( 

707 imRef, 

708 imDef 

709 ) 

710 """ 

711 import multiprocessing 

712 

713 try: 

714 multiprocessing.set_start_method("fork") 

715 except RuntimeError: 

716 pass 

717 

718 import spam.helpers 

719 import spam.mesh 

720 

721 # Global number of processes 

722 nThreads = multiprocessing.cpu_count() if nThreads is None else nThreads 

723 print(f"[globalCorrelation] C++ parallelisation on {nThreads} threads") 

724 

725 print(f"[globalCorrelation] Convergence criterion = {convergenceCriterion}") 

726 print(f"[globalCorrelation] Max iterations = {maxIterations}") 

727 print("[globalCorrelation] Converting im1 to 32-bit float") 

728 im1 = im1.astype("<f4") 

729 

730 points = points.astype("<f8") 

731 connectivity = connectivity.astype("<u4") 

732 

733 maxCoord = numpy.amax(points, axis=0).astype("<u2") 

734 minCoord = numpy.amin(points, axis=0).astype("<u2") 

735 print(f"[globalCorrelation] Mesh box: min = {minCoord} max = {maxCoord}") 

736 

737 meshPaddingSlice = ( 

738 slice(minCoord[0], maxCoord[0]), 

739 slice(minCoord[1], maxCoord[1]), 

740 slice(minCoord[2], maxCoord[2]), 

741 ) 

742 

743 displacements = numpy.zeros((points.shape[0], 3), dtype="<f8") 

744 

745 print(f"[globalCorrelation] Points: {points.shape}") 

746 print(f"[globalCorrelation] Displacements: {displacements.shape}") 

747 print(f"[globalCorrelation] Cells: {connectivity.shape}") 

748 print(f"[globalCorrelation] Padding: {meshPaddingSlice}") 

749 

750 ############################################################### 

751 # Step 2-1 Apply deformation and interpolate pixels 

752 ############################################################### 

753 

754 print("[globalCorrelation] Allocating 3D data (deformed image)") 

755 if initialDisplacements is None: 

756 im1Def = im1.copy() 

757 imTetLabel = spam.label.labelTetrahedra(im1.shape, points, connectivity, nThreads=nThreads) 

758 else: 

759 print("[globalCorrelation] Applying initial deformation to image") 

760 displacements = initialDisplacements.copy() 

761 tic = time.perf_counter() 

762 imTetLabel = spam.label.labelTetrahedra(im1.shape, points + displacements, connectivity, nThreads=nThreads) 

763 print(f"[globalCorrelation] Running labelTetrahedra: {time.perf_counter()-tic:.3f} seconds.") 

764 

765 im1Def = spam.DIC.applyMeshTransformation( 

766 im1, 

767 points, 

768 connectivity, 

769 displacements, 

770 imTetLabel=imTetLabel, 

771 nThreads=nThreads, 

772 ) 

773 if debugFiles: 

774 print("[globalCorrelation] Saving initial images") 

775 for name, image in [ 

776 [f"{prefix}-def-init.tif", im1Def], 

777 [f"{prefix}-imTetLabel-init.tif", imTetLabel], 

778 ]: 

779 print(f"[globalCorrelation]\t{name}: {image.shape}") 

780 tifffile.imwrite(name, image) 

781 

782 # print("[globalCorrelation] Correlating (MF)!") 

783 print("[globalCorrelation] Calculating gradient of IM TWO...") 

784 im2Grad = numpy.array(numpy.gradient(im2), dtype="<f4") 

785 

786 print("[globalCorrelation] Computing global matrix") 

787 # This generates the globalMatrix (big Mc matrix) with imGrad as input 

788 Mc = numpy.zeros((3 * points.shape[0], 3 * points.shape[0]), dtype="<f8") 

789 

790 if debugFiles: 

791 print("[globalCorrelation] Computing debug files fields") 

792 gradientPerTet = numpy.zeros((connectivity.shape[0], 3), dtype="<f8") 

793 IDPerTet = numpy.array([_ for _ in range(connectivity.shape[0])]) 

794 

795 computeGradientPerTet( 

796 imTetLabel.astype("<u4"), 

797 im2Grad.astype("<f4"), 

798 connectivity.astype("<u4"), 

799 (points + displacements).astype("<f8"), 

800 gradientPerTet, 

801 ) 

802 

803 spam.helpers.writeUnstructuredVTK( 

804 (points + displacements), 

805 connectivity, 

806 cellData={"meanGradient": gradientPerTet, "id": IDPerTet}, 

807 fileName=f"{prefix}-gradient.vtk", 

808 ) 

809 del gradientPerTet 

810 

811 computeDICglobalMatrix( 

812 imTetLabel.astype("<u4"), 

813 im2Grad.astype("<f4"), 

814 connectivity.astype("<u4"), 

815 (points + displacements).astype("<f8"), 

816 Mc, 

817 ) 

818 

819 ############################################################### 

820 # Setup left hand vector 

821 ############################################################### 

822 if regularisationMatrix: 

823 regularisationField = isochoricField(points) if regularisationField is None else regularisationField 

824 Ec = _computeFunctional(regularisationField, Mc) 

825 regularisationMatrix *= Ec 

826 print(f"[regularisation] Wc/Ec = {1/Ec:.3e}") 

827 left_hand_inverse = numpy.linalg.inv(Mc + regularisationMatrix) 

828 else: 

829 print("[globalCorrelation] Skip regularisation") 

830 left_hand_inverse = numpy.linalg.inv(Mc) 

831 del Mc 

832 

833 # error = _errorCalc(im2, im1Def, im1, meshPaddingSlice) 

834 # print("\[globalCorrelation] Initial Error (abs) = ", error) 

835 

836 # We try to solve Md=F 

837 # while error > 0.1 and error < errorIn: 

838 # while error > 0.1 and i <= maxIterations and error < errorIn: 

839 dxNorm = numpy.inf 

840 i = 0 

841 while dxNorm > convergenceCriterion and i < maxIterations: 

842 i += 1 

843 

844 # This function returns globalVector (F) taking in im1Def and im2 and the gradients 

845 tic = time.perf_counter() 

846 # print("[globalCorrelation] [newton] run computeDICglobalVector: ", end="") 

847 right_hand_vector = numpy.zeros((3 * points.shape[0]), dtype="<f8") 

848 computeDICglobalVector( 

849 imTetLabel.astype("<u4"), 

850 im2Grad.astype("<f4"), 

851 im1Def.astype("<f4"), 

852 im2.astype("<f4"), 

853 connectivity.astype("<u4"), 

854 (points + displacements).astype("<f8"), 

855 right_hand_vector, 

856 ) 

857 # print(f"{time.perf_counter()-tic:.3f} seconds.") 

858 

859 tic = time.perf_counter() 

860 # print("[globalCorrelation] [newton] run solve: ", end="") 

861 

862 # solve: we can use solve here for sake of precision (over computing 

863 # M^-1). However solve takes quite a lot of time for "small" meshes). 

864 

865 if regularisationMatrix: 

866 right_hand_vector -= numpy.matmul(regularisationMatrix, displacements.ravel()) 

867 dx = numpy.matmul(left_hand_inverse, right_hand_vector).astype("<f8") 

868 # dx_solve = numpy.linalg.solve( 

869 # Mc, 

870 # right_hand_vector 

871 # ).astype('<f8') 

872 # print(numpy.linalg.norm(dx - dx_solve)) 

873 

874 displacements += dx.reshape(points.shape[0], 3) 

875 dxNorm = numpy.linalg.norm(dx) 

876 # print(f"{time.perf_counter()-tic:.3f} seconds.") 

877 

878 if medianFilterEachIteration: 

879 # use connectivity to filter 

880 print("[globalCorrelation] [newton] Median filter of displacements...") 

881 for nodeNumber in range(points.shape[0]): 

882 # get rows of connectivity (i.e., tets) which include this point 

883 connectedTets = numpy.where(connectivity == nodeNumber)[0] 

884 neighbourPoints = numpy.unique(connectivity[connectedTets]) 

885 diff = numpy.median(displacements[neighbourPoints], axis=0) - displacements[nodeNumber] 

886 displacements[nodeNumber] += 0.5 * diff 

887 

888 tic = time.perf_counter() 

889 # print("[globalCorrelation] [newton] run labelTetrahedra: ", end="") 

890 

891 imTetLabel = spam.label.labelTetrahedra(im1.shape, points + displacements, connectivity, nThreads=nThreads) 

892 # print(f"{time.perf_counter()-tic:.3f} seconds.") 

893 

894 tic = time.perf_counter()