Coverage for src/dictk/element.py: 100%

1"""A single four-noded quadrilateral (Q4) finite element: shape functions,

2Jacobian, deformation gradient, and strain at its Gauss points.

4The shape-function/Jacobian/deformation-gradient chain and the

5Green-Lagrange strain measure are ported from the `hdic` codebase's

6`~/hdic/src/hdic/types/fea.py` (see [[project_dictk_hdic_provenance]]),

7which was one hdic FEA implementation actually wired into a

8real pipeline there, and validated end to end against a hand-worked

9example (a unit square stretched 5% in $x$ gives $E_{11} = 0.05125$

10exactly), reproduced here as this module's own regression test.

12The logarithmic (Hencky) strain measure has no hdic counterpart --

13hdic never implemented it -- so it's dictk-original, following the

14[Seth-Hill Strain

15Family](../getting_started/continuum_mechanics.html#seth-hill-strain-family)'s

16$m=0$ case already documented in `continuum_mechanics.md`.

18See [Finite Element

19Method](../getting_started/finite_element_method.html) for the full

20derivation this module implements -- each function's docstring links to

21the specific subsection its math comes from.

22"""

24from collections.abc import Sequence

25from typing import Final

27import numpy as np

29from dictk.image import PixelCoordinate, SubpixelCoordinate

31#: The 2x2 Gauss quadrature rule's local coordinate value -- see

32#: [Gauss Points](../getting_started/finite_element_method.html#gauss-points).

33GAUSS_POINT_COORDINATE: Final[float] = 1.0 / np.sqrt(3.0)

36def gauss_points() -> list[tuple[float, float]]:

37 r"""The four $(\xi, \eta)$ locations of the standard 2x2 Gauss rule.

39 See [Gauss

40 Points](../getting_started/finite_element_method.html#gauss-points).

41 In the standard FEA local coordinate system, the y-axis points up,

42 different from the image coordinate system where the y-axis points down.

43 For both coordinate systems, the x-axis points right.

44 Ordered to match `shape_functions`' own $N_1$..$N_4$ corner

45 convention (bottom-left, bottom-right, top-right, top-left in local

46 coordinates), for a y-axis-up convention:

48 ```

49 eta

50 ^

51 |

52 N4 (-1, 1) o------+------o N3 (1, 1)

53 | | |

54 | g4 | g3 |

55 | | |

56 |------+------|---> xi

57 | | |

58 | g1 | g2 |

59 | | |

60 N1 (-1,-1) o------+------o N2 (1,-1)

61 ```

63 `N1`..`N4` are the element's 4 corner nodes, at the local coordinates

64 labeled above. `g1`..`g4` are this function's own 4 returned

65 `(xi, eta)` pairs, in return order -- each at

66 $(\pm 1/\sqrt{3}, \pm 1/\sqrt{3})$, inset from its nearest corner

67 node.

69 Returns:

70 A length-4 list of `(xi, eta)` pairs.

71 """

72 g = GAUSS_POINT_COORDINATE

73 return [(-g, -g), (g, -g), (g, g), (-g, g)]

76def gauss_point_coordinates(

77 *, points: Sequence[PixelCoordinate | SubpixelCoordinate]

78) -> list[tuple[float, float]]:

79 r"""An element's 4 Gauss points' own global $(X, Y)$ position.

81 dictk-original -- no hdic counterpart. Computed as

82 `shape_functions(xi, eta) @ points` at each of `gauss_points()`'s 4

83 locations, so each Gauss point's position is a bilinear interpolation

84 of the element's 4 corner nodes, not the corners themselves.

86 Returns `(float, float)` pairs, not `PixelCoordinate` -- a Gauss

87 point's position is generally sub-pixel (bilinear interpolation of

88 integer corners rarely lands back on an integer), and

89 `PixelCoordinate.x`/`.y` are `int`, which would silently truncate it.

91 Args:

92 points: The element's 4 corner nodes' positions, in $N_1$..$N_4$

93 order (see [Shape

94 Functions](../getting_started/finite_element_method.html#shape-functions)).

95 Pass the *current* (deformed) positions to place each Gauss

96 point at its current-configuration location -- e.g. for

97 overlaying on a deformed image, matching the reference

98 hdic implementation this pattern is modeled on

99 (`~/hdic/src/hdic/types/fea_vis.py`'s

100 `plot_strain_at_gauss_points`, which always uses deformed

101 node positions).

102

103 Returns:

104 A length-4 list of `(X, Y)` pairs, in `gauss_points()`'s order.

105

106 Raises:

107 ValueError: If `points` is not length 4.

108 """

109 if len(points) != 4:

110 raise ValueError(f"points has {len(points)} points, must be exactly 4")

111

112 coordinates = np.array([[p.x, p.y] for p in points], dtype=float)

113 return [

114 tuple(shape_functions(xi=xi, eta=eta) @ coordinates)

115 for xi, eta in gauss_points()

116 ]

117

118

119def shape_functions(*, xi: float, eta: float) -> np.ndarray:

120 r"""The 4 Q4 shape functions $N_1$..$N_4$ at local coordinate $(\xi, \eta)$.

121

122 See [Shape

123 Functions](../getting_started/finite_element_method.html#shape-functions).

124

125 Args:

126 xi: Local coordinate along $\xi$, in $[-1, 1]$.

127 eta: Local coordinate along $\eta$, in $[-1, 1]$.

128

129 Returns:

130 `N`, shape `(4,)`: $N_1, N_2, N_3, N_4$.

131 """

132 return 0.25 * np.array(

133 [

134 (1 - xi) * (1 - eta),

135 (1 + xi) * (1 - eta),

136 (1 + xi) * (1 + eta),

137 (1 - xi) * (1 + eta),

138 ]

139 )

140

141

142def shape_function_derivatives(*, xi: float, eta: float) -> np.ndarray:

143 r"""The shape functions' derivatives with respect to local coordinates $(\xi, \eta)$.

144

145 See [Shape Function Derivatives in Local

146 Coordinates](../getting_started/finite_element_method.html#shape-function-derivatives-in-local-coordinates).

147

148 Args:

149 xi: Local coordinate along $\xi$, in $[-1, 1]$.

150 eta: Local coordinate along $\eta$, in $[-1, 1]$.

151

152 Returns:

153 `dN/d(xi, eta)`, shape `(2, 4)`: row 0 is $\partial N_a/\partial

154 \xi$ for $a=1..4$, row 1 is $\partial N_a/\partial \eta$.

155 """

156 return 0.25 * np.array(

157 [

158 [-(1 - eta), (1 - eta), (1 + eta), -(1 + eta)],

159 [-(1 - xi), -(1 + xi), (1 + xi), (1 - xi)],

160 ]

161 )

162

163

164def jacobian(*, derivatives: np.ndarray, coordinates: np.ndarray) -> np.ndarray:

165 r"""The Jacobian matrix $\boldsymbol{j}_0$ mapping local to global coordinate derivatives.

166

167 See [Jacobian

168 Matrix](../getting_started/finite_element_method.html#jacobian-matrix).

169

170 Args:

171 derivatives: `shape_function_derivatives`' own output, shape `(2, 4)`.

172 coordinates: The element's 4 corner nodes' positions, shape

173 `(4, 2)`, each row `[X_a, Y_a]`, in $N_1$..$N_4$ order.

174

175 Returns:

176 $\boldsymbol{j}_0$ = `derivatives @ coordinates`, shape `(2, 2)`.

177 """

178 return derivatives @ coordinates

179

180

181def shape_function_gradients(

182 *, derivatives: np.ndarray, jacobian: np.ndarray

183) -> np.ndarray:

184 r"""The shape functions' derivatives with respect to global coordinates $(X, Y)$.

185

186 See [Shape Function Derivatives in Global

187 Coordinates](../getting_started/finite_element_method.html#shape-function-derivatives-in-global-coordinates).

188

189 Args:

190 derivatives: `shape_function_derivatives`' own output, shape `(2, 4)`.

191 jacobian: `jacobian()`'s own output for this same element and

192 $(\xi, \eta)$, shape `(2, 2)`.

193

194 Returns:

195 $\partial N_a/\partial X$ = $\boldsymbol{j}_0^{-1}$ `@ derivatives`,

196 shape `(2, 4)`.

197

198 Raises:

199 ValueError: If `jacobian` is singular (a degenerate element, e.g.

200 two coincident corners) -- hdic's own `types/fea.py` doesn't

201 guard this; this port does.

202 """

203 if np.isclose(np.linalg.det(jacobian), 0.0):

204 raise ValueError(

205 f"jacobian {jacobian.tolist()} is singular -- the element is degenerate "

206 "(e.g. two coincident corners), so shape function gradients are undefined"

207 )

208 return np.linalg.inv(jacobian) @ derivatives

209

210

211def displacement_gradient(

212 *, gradients: np.ndarray, displacements: np.ndarray

213) -> np.ndarray:

214 r"""The displacement field's gradient $\boldsymbol{\nabla}_0\boldsymbol{u}$ at a Gauss point.

215

216 See [Displacement

217 Gradient](../getting_started/finite_element_method.html#displacement-gradient).

218

219 Args:

220 gradients: `shape_function_gradients`' own output, shape `(2, 4)`.

221 displacements: The element's 4 corner nodes' displacements

222 (current position minus reference position), shape `(4, 2)`,

223 each row `[u_a, v_a]`, same $N_1$..$N_4$ order as `coordinates`.

224

225 Returns:

226 $\boldsymbol{\nabla}_0\boldsymbol{u}$ = `(gradients @

227 displacements).T`, shape `(2, 2)`.

228 """

229 return (gradients @ displacements).T

230

231

232def deformation_gradient(*, displacement_gradient: np.ndarray) -> np.ndarray:

233 r"""The deformation gradient $\boldsymbol{F} = \boldsymbol{I} + \boldsymbol{\nabla}_0\boldsymbol{u}$.

234

235 See [Deformation

236 Gradient](../getting_started/finite_element_method.html#deformation-gradient).

237

238 Args:

239 displacement_gradient: `displacement_gradient()`'s own output, shape `(2, 2)`.

240

241 Returns:

242 $\boldsymbol{F}$, shape `(2, 2)`.

243 """

244 return np.eye(2) + displacement_gradient

245

246

247def green_lagrange_strain(*, deformation_gradient: np.ndarray) -> np.ndarray:

248 r"""The Green-Lagrange strain tensor $\boldsymbol{E} = \frac{1}{2}(\boldsymbol{F}^T\boldsymbol{F} - \boldsymbol{I})$.

249

250 See [Green-Lagrange

251 Strain](../getting_started/continuum_mechanics.html#green-lagrange-strain).

252

253 Args:

254 deformation_gradient: `deformation_gradient()`'s own output, shape `(2, 2)`.

255

256 Returns:

257 $\boldsymbol{E}$, shape `(2, 2)`.

258 """

259 return 0.5 * (deformation_gradient.T @ deformation_gradient - np.eye(2))

260

261

262def log_strain(*, deformation_gradient: np.ndarray) -> np.ndarray:

263 r"""The logarithmic (Hencky, natural) strain tensor $\boldsymbol{E}^{(0)} = \ln\boldsymbol{U}$.

264

265 dictk-original -- hdic's `types/fea.py` only implements

266 Green-Lagrange strain, so there is no hdic source to port here.

267 Added to compare against VIC-2D, which reports logarithmic (Euler)

268 strain, not Green-Lagrange (see [Verification Against

269 VIC-2D](../getting_started/simple_stretch.html#verification-against-vic-2d)).

270

271 Computed via the [Seth-Hill Strain

272 Family](../getting_started/continuum_mechanics.html#seth-hill-strain-family)'s

273 $m=0$ case, using its [spectral

274 representation](../getting_started/continuum_mechanics.html#spectral-representation):

275 eigendecompose the right Cauchy-Green tensor

276 $\boldsymbol{C} = \boldsymbol{F}^T\boldsymbol{F}$ for principal

277 stretches $\lambda_\alpha$ and principal directions

278 $\boldsymbol{N}_\alpha$, then

279 $\boldsymbol{E}^{(0)} = \sum_\alpha \ln(\lambda_\alpha)\,

280 \boldsymbol{N}_\alpha \otimes \boldsymbol{N}_\alpha$.

281

282 Args:

283 deformation_gradient: `deformation_gradient()`'s own output, shape `(2, 2)`.

284

285 Returns:

286 $\boldsymbol{E}^{(0)}$, shape `(2, 2)`.

287 """

288 right_cauchy_green = deformation_gradient.T @ deformation_gradient

289 eigenvalues, eigenvectors = np.linalg.eigh(right_cauchy_green)

290 stretches = np.sqrt(eigenvalues)

291 return eigenvectors @ np.diag(np.log(stretches)) @ eigenvectors.T

292

293

294def _gauss_point_deformation_gradients(

295 *,

296 reference_points: Sequence[PixelCoordinate | SubpixelCoordinate],

297 current_points: Sequence[PixelCoordinate | SubpixelCoordinate],

298) -> list[np.ndarray]:

299 r"""The deformation gradient $\boldsymbol{F}$ at an element's 4 Gauss points.

300

301 Shared plumbing for `gauss_point_green_lagrange_strains` and `gauss_point_log_strains`:

302 composes `shape_function_derivatives` -> `jacobian` ->

303 `shape_function_gradients` -> `displacement_gradient` ->

304 `deformation_gradient` at each of `gauss_points()`'s 4 locations,

305 given the element's 4 corner nodes' reference and current positions

306 directly -- unlike the hdic chain this is ported from, which requires

307 the caller to pre-subtract raw numpy arrays into a `displacements`

308 matrix before calling in.

309

310 Args:

311 reference_points: The element's 4 corner nodes' reference

312 positions, in $N_1$..$N_4$ order (see [Shape

313 Functions](../getting_started/finite_element_method.html#shape-functions)).

314 current_points: The same 4 nodes' current positions, same order

315 and indexing as `reference_points`.

316

317 Returns:

318 A length-4 list of `(2, 2)` deformation gradients, one per Gauss

319 point, in `gauss_points()`'s order.

320

321 Raises:

322 ValueError: If `reference_points` or `current_points` is not

323 length 4, or the Jacobian is singular at some Gauss point

324 (see `shape_function_gradients`).

325 """

326 if len(reference_points) != 4:

327 raise ValueError(

328 f"reference_points has {len(reference_points)} points, must be exactly 4"

329 )

330 if len(current_points) != 4:

331 raise ValueError(

332 f"current_points has {len(current_points)} points, must be exactly 4"

333 )

334

335 coordinates = np.array([[p.x, p.y] for p in reference_points], dtype=float)

336 current = np.array([[p.x, p.y] for p in current_points], dtype=float)

337 displacements = current - coordinates

338

339 gradients_at_gauss_points = []

340 for xi, eta in gauss_points():

341 derivatives = shape_function_derivatives(xi=xi, eta=eta)

342 j = jacobian(derivatives=derivatives, coordinates=coordinates)

343 gradients = shape_function_gradients(derivatives=derivatives, jacobian=j)

344 grad_u = displacement_gradient(gradients=gradients, displacements=displacements)

345 gradients_at_gauss_points.append(

346 deformation_gradient(displacement_gradient=grad_u)

347 )

348 return gradients_at_gauss_points

349

350

351def gauss_point_green_lagrange_strains(

352 *,

353 reference_points: Sequence[PixelCoordinate | SubpixelCoordinate],

354 current_points: Sequence[PixelCoordinate | SubpixelCoordinate],

355) -> list[np.ndarray]:

356 r"""Green-Lagrange strain at an element's 4 Gauss points.

357

358 Args:

359 reference_points: The element's 4 corner nodes' reference

360 positions, in $N_1$..$N_4$ order (see [Shape

361 Functions](../getting_started/finite_element_method.html#shape-functions)).

362 current_points: The same 4 nodes' current positions, same order

363 and indexing as `reference_points`.

364

365 Returns:

366 A length-4 list of `(2, 2)` Green-Lagrange strain tensors, one

367 per Gauss point, in `gauss_points()`'s order.

368

369 Raises:

370 ValueError: If `reference_points` or `current_points` is not

371 length 4, or the Jacobian is singular at some Gauss point

372 (see `shape_function_gradients`).

373 """

374 return [

375 green_lagrange_strain(deformation_gradient=f)

376 for f in _gauss_point_deformation_gradients(

377 reference_points=reference_points, current_points=current_points

378 )

379 ]

380

381

382def gauss_point_log_strains(

383 *,

384 reference_points: Sequence[PixelCoordinate | SubpixelCoordinate],

385 current_points: Sequence[PixelCoordinate | SubpixelCoordinate],

386) -> list[np.ndarray]:

387 r"""Logarithmic (Hencky) strain at an element's 4 Gauss points.

388

389 Same composition as `gauss_point_green_lagrange_strains`, substituting `log_strain`

390 for `green_lagrange_strain` as the final step.

391

392 Args:

393 reference_points: The element's 4 corner nodes' reference

394 positions, in $N_1$..$N_4$ order (see [Shape

395 Functions](../getting_started/finite_element_method.html#shape-functions)).

396 current_points: The same 4 nodes' current positions, same order

397 and indexing as `reference_points`.

398

399 Returns:

400 A length-4 list of `(2, 2)` logarithmic strain tensors, one per

401 Gauss point, in `gauss_points()`'s order.

402

403 Raises:

404 ValueError: If `reference_points` or `current_points` is not

405 length 4, or the Jacobian is singular at some Gauss point

406 (see `shape_function_gradients`).

407 """

408 return [

409 log_strain(deformation_gradient=f)

410 for f in _gauss_point_deformation_gradients(

411 reference_points=reference_points, current_points=current_points

412 )

413 ]