Coverage for /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pysagas/geometry/cell.py: 95%
192 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-30 04:27 +0000
1from __future__ import annotations
2import numpy as np
3import pysagas.flow
4from pysagas.geometry import Vector
5from numpy.typing import ArrayLike
6from typing import Union, List, Optional
9class Cell:
10 """A triangular cell object.
12 Attributes
13 -----------
14 p0 : Vector
15 The first vertex of the cell.
17 p1 : Vector
18 The second vertex of the cell.
20 p2 : Vector
21 The third vertex of the cell.
23 A : float
24 The cell face area.
26 n : Vector
27 The cell normal.
29 c : Vector
30 The cell centroid.
32 dndv : np.array
33 The sensitivity of the cell's normal vector to each vertex.
35 dAdv : np.array
36 The sensitivity of the cell's area to each vertex.
38 dvdp : np.array
39 The sensitivity of the cell's vertices to each geometric parameter.
41 dndp : np.array
42 The sensitivity of the cell's normal vector to each geometric parameter.
44 dAdp : np.array
45 The sensitivity of the cell's area to each geometric parameter.
47 dcdp : np.array
48 The sensitivity of the centroid to each geometric parameter.
50 flowstate : FlowState
51 The flowstate associated with the cell.
53 sensitivities : np.array
54 An array containing the [x,y,z] force sensitivities of the cell.
55 """
57 def __init__(
58 self, p0: Vector, p1: Vector, p2: Vector, face_ids: Optional[list[int]] = None
59 ):
60 """Constructs a cell, defined by three points.
62 Parameters
63 ----------
64 p0 : Vector
65 The first point defining the cell.
67 p1 : Vector
68 The second point defining the cell.
70 p2 : Vector
71 The third point defining the cell.
72 """
73 # Save points
74 self.p0 = p0
75 self.p1 = p1
76 self.p2 = p2
78 # Save connectivity information (PyMesh)
79 self._face_ids = face_ids
81 # Calculate cell properites
82 self.n = -self.calc_normal(p0, p1, p2)
83 self.A = self.calc_area(p0, p1, p2)
84 self.c = self.calc_centroid(p0, p1, p2)
86 # Calculate geometric sensitivities
87 self._dndv: ArrayLike = None
88 self._dAdv: ArrayLike = None
89 self._dcdv: ArrayLike = None
91 # Parameter sensitivities
92 self.dvdp: ArrayLike = None # vertex-parameter sensitivities
93 self.dndp: ArrayLike = None # normal-parameter sensitivities
94 self.dAdp: ArrayLike = None # area-parameter sensitivities
95 self.dcdp: ArrayLike = None # centroid-parameter sensitivities
97 # FlowState
98 self.flowstate: pysagas.flow.FlowState = None
100 # Sensitivities
101 self.sensitivities = None
102 self.moment_sensitivities = None
104 # Cell attributes
105 self.attributes = {}
107 @property
108 def dndv(self):
109 if self._dndv is not None:
110 # Already exists, return it
111 return self._dndv
112 else:
113 # Calculate and return it
114 self._dndv = self.n_sensitivity(self.p0, self.p1, self.p2)
115 return self._dndv
117 @property
118 def dAdv(self):
119 if self._dAdv is not None:
120 # Already exists, return it
121 return self._dAdv
122 else:
123 # Calculate and return it
124 self._dAdv = self.A_sensitivity(self.p0, self.p1, self.p2)
125 return self._dAdv
127 @property
128 def dcdv(self):
129 if self._dcdv is not None:
130 # Already exists, return it
131 return self._dcdv
132 else:
133 # Calculate and return it
134 self._dcdv = self.c_sensitivity(self.p0, self.p1, self.p2)
135 return self._dcdv
137 @property
138 def vertices(self):
139 """The cell vertices."""
140 return np.array([getattr(self, p).vec for p in ["p0", "p1", "p2"]])
142 def to_dict(self):
143 """Returns the Cell as a dictionary."""
144 # TODO - need to think about representing the cells
145 # vs. the points
146 # The values at the points will need to be averaged
147 # across cells
148 pass
150 @classmethod
151 def from_points(
152 cls, points: Union[List[Vector], np.array[Vector]], **kwargs
153 ) -> Cell:
154 """Constructs a Vector object from an array of coordinates.
156 Parameters
157 ----------
158 points : Union[List[Vector], np.array[Vector]]
159 The points defining the cell.
161 Returns
162 -------
163 Cell
164 """
165 return cls(*points, **kwargs)
167 def __repr__(self) -> str:
168 return f"Cell({self.p0}, {self.p1}, {self.p2})"
170 def __str__(self) -> str:
171 return "A Cell"
173 def _add_sensitivities(self, dvdp: np.array) -> None:
174 """Adds the cell sensitivities to the cell.
176 Parameters
177 -----------
178 dvdp : array
179 The sensitivity of each vertex x, y and z, with respect to each
180 parameter. The dimensions of dvdp are (9, p), for p parameters.
181 Each 3 rows corresponds to the x, y and z sensitivities for each
182 vertex of the cell. Each column is for the relevant parameter.
183 """
184 self.dvdp = dvdp
185 self.dndp = np.dot(self.dndv, self.dvdp)
186 self.dAdp = np.dot(self.dAdv, self.dvdp)
187 self.dcdp = np.dot(self.dcdv, self.dvdp)
189 @staticmethod
190 def calc_normal(p0: Vector, p1: Vector, p2: Vector) -> Vector:
191 """Calculates the normal vector of a cell defined by three points.
193 Parameters
194 ----------
195 p0 : Vector
196 The first point defining the cell.
198 p1 : Vector
199 The second point defining the cell.
201 p2 : Vector
202 The third point defining the cell.
204 Returns
205 --------
206 normal : Vector
207 The unit normal vector of the cell defined by the points.
209 References
210 -----------
211 https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal
212 """
213 # Define cell edges
214 A = p1 - p0
215 B = p2 - p0
217 # Calculate normal components
218 Nx = A.y * B.z - A.z * B.y
219 Ny = A.z * B.x - A.x * B.z
220 Nz = A.x * B.y - A.y * B.x
222 # Construct normal vector
223 normal = Vector(x=Nx, y=Ny, z=Nz)
225 if normal.norm == 0:
226 # Degenerate cell
227 raise DegenerateCell()
229 # Convert to unit vector
230 unit_normal = normal / np.linalg.norm(normal.vec)
232 return unit_normal
234 @staticmethod
235 def calc_area(p0: Vector, p1: Vector, p2: Vector) -> float:
236 """Calculates the area of a cell defined by three points.
238 Parameters
239 ----------
240 p0 : Vector
241 The first point defining the cell.
243 p1 : Vector
244 The second point defining the cell.
246 p2 : Vector
247 The third point defining the cell.
249 Returns
250 --------
251 area : float
252 The area of the cell defined by the points.
254 References
255 -----------
256 https://en.wikipedia.org/wiki/Cross_product
257 """
258 # Define cell edges
259 A = p1 - p0
260 B = p2 - p0
262 # Calculate area
263 S = 0.5 * np.linalg.norm(np.cross(A.vec, B.vec))
264 return S
266 @staticmethod
267 def calc_centroid(p0: Vector, p1: Vector, p2: Vector) -> Vector:
268 """Calculates the centroid of a cell defined by three points.
270 Parameters
271 ----------
272 p0 : Vector
273 The first point defining the cell.
275 p1 : Vector
276 The second point defining the cell.
278 p2 : Vector
279 The third point defining the cell.
281 Returns
282 --------
283 c : Vector
284 The centroid of the cell defined by the points.
286 References
287 -----------
288 https://en.wikipedia.org/wiki/Centroid
289 """
290 cx = (p0.x + p1.x + p2.x) / 3
291 cy = (p0.y + p1.y + p2.y) / 3
292 cz = (p0.z + p1.z + p2.z) / 3
293 return Vector(cx, cy, cz)
295 @staticmethod
296 def n_sensitivity(p0: Vector, p1: Vector, p2: Vector) -> np.array:
297 """Calculates the sensitivity of a cell's normal vector
298 to the points defining the cell analytically.
300 Parameters
301 ----------
302 p0 : Vector
303 The first point defining the cell.
305 p1 : Vector
306 The second point defining the cell.
308 p2 : Vector
309 The third point defining the cell.
311 Returns
312 -------
313 sensitivity : np.array
314 The sensitivity matrix with size m x n x p. Rows m refer to
315 the vertices, columns n refer to the vertex coordinates, and
316 slices p refer to the components of the normal vector.
317 """
318 g = (
319 ((p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)) ** 2
320 + ((p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)) ** 2
321 + ((p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)) ** 2
322 )
324 # n = np.empty(3)
325 # A = p0
326 # B = p1
327 # C = p2
328 # n[0] = ((B.y-A.y)*(C.z-A.z) - (B.z-A.z)*(C.y-A.y)) / g**0.5
329 # n[1] = ((B.z-A.z)*(C.x-A.x) - (B.x-A.x)*(C.z-A.z)) / g**0.5
330 # n[2] = ((B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x)) / g**0.5
332 M_sense = np.empty((3, 9))
333 for row in range(3):
334 if row == 0:
335 f = (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y) # f_x
336 elif row == 1:
337 f = (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z) # f_y
338 elif row == 2:
339 f = (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x) # f_z
341 for col in range(9):
342 if col == 0: # d/dp0.x
343 g_dash = 2 * (
344 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
345 ) * (-(p1.z - p0.z) + (p2.z - p0.z)) + 2 * (
346 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
347 ) * (
348 -(p2.y - p0.y) + (p1.y - p0.y)
349 )
350 elif col == 1: # d/dp0.y
351 g_dash = 2 * (
352 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
353 ) * (-(p2.z - p0.z) + (p1.z - p0.z)) + 2 * (
354 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
355 ) * (
356 -(p1.x - p0.x) + (p2.x - p0.x)
357 )
358 elif col == 2: # d/dp0.z
359 g_dash = 2 * (
360 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
361 ) * (-(p1.y - p0.y) + (p2.y - p0.y)) + 2 * (
362 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
363 ) * (
364 -(p2.x - p0.x) + (p1.x - p0.x)
365 )
366 elif col == 3: # d/dp1.x
367 g_dash = 2 * (
368 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
369 ) * (-(p2.z - p0.z)) + 2 * (
370 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
371 ) * (
372 (p2.y - p0.y)
373 )
374 elif col == 4: # d/dp1.y
375 g_dash = 2 * (
376 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
377 ) * ((p2.z - p0.z)) + 2 * (
378 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
379 ) * (
380 -(p2.x - p0.x)
381 )
382 elif col == 5: # d/dp1.z
383 g_dash = 2 * (
384 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
385 ) * (-(p2.y - p0.y)) + 2 * (
386 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
387 ) * (
388 (p2.x - p0.x)
389 )
390 elif col == 6: # d/dp2.x
391 g_dash = 2 * (
392 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
393 ) * ((p1.z - p0.z)) + 2 * (
394 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
395 ) * (
396 -(p1.y - p0.y)
397 )
398 elif col == 7: # d/dp2.y
399 g_dash = 2 * (
400 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
401 ) * (-(p1.z - p0.z)) + 2 * (
402 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
403 ) * (
404 (p1.x - p0.x)
405 )
406 elif col == 8: # d/dp2.z
407 g_dash = 2 * (
408 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
409 ) * ((p1.y - p0.y)) + 2 * (
410 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
411 ) * (
412 -(p1.x - p0.x)
413 )
415 if row == 0: # for n..x
416 # ((p1.y-p0.y)*(p2.z-p0.z) - (p1.z-p0.z)*(p2.y-p0.y)) # f_x
417 if col == 1: # d/dp0.y
418 f_dash = -(p2.z - p0.z) + (p1.z - p0.z)
419 elif col == 2: # d/dp0.z
420 f_dash = -(p1.y - p0.y) + (p2.y - p0.y)
421 elif col == 4: # d/dp1.y
422 f_dash = p2.z - p0.z
423 elif col == 5: # d/dp1.z
424 f_dash = -(p2.y - p0.y)
425 elif col == 7: # d/dp2.y
426 f_dash = -(p1.z - p0.z)
427 elif col == 8: # d/dp2.z
428 f_dash = p1.y - p0.y
429 else:
430 f_dash = 0
431 if row == 1: # for n.y # CHECK
432 # f = ((p1.z-p0.z)*(p2.x-p0.x) - (p1.x-p0.x)*(p2.z-p0.z)) # f_y
433 if col == 0: # d/dp0.x
434 f_dash = -(p1.z - p0.z) + (p2.z - p0.z)
435 elif col == 2: # d/dp0.z
436 f_dash = -(p2.x - p0.x) + (p1.x - p0.x) # CHECK
437 elif col == 3: # d/dp1.x
438 f_dash = -(p2.z - p0.z)
439 elif col == 5: # d/dp1.z
440 f_dash = p2.x - p0.x
441 elif col == 6: # d/dp2.x
442 f_dash = p1.z - p0.z
443 elif col == 8: # d/dp2.z
444 f_dash = -(p1.x - p0.x)
445 else:
446 f_dash = 0
447 if row == 2: # for n.z # CHECK
448 # f = ((p1.x-p0.x)*(p2.y-p0.y) - (p1.y-p0.y)*(p2.x-p0.x)) # f_z
449 if col == 0: # d/dp0.x
450 f_dash = -(p2.y - p0.y) + (p1.y - p0.y)
451 elif col == 1: # d/dp0.y
452 f_dash = -(p1.x - p0.x) + (p2.x - p0.x)
453 elif col == 3: # d/dp1.x
454 f_dash = p2.y - p0.y
455 elif col == 4: # d/dp1.y
456 f_dash = -(p2.x - p0.x)
457 elif col == 6: # d/dp2.x
458 f_dash = -(p1.y - p0.y)
459 elif col == 7: # d/dp2.y
460 f_dash = p1.x - p0.x
461 else:
462 f_dash = 0
464 # Assign sensitivity
465 M_sense[row, col] = (
466 f_dash * (g ** (-0.5)) + f * (-1 / 2) * g ** (-3 / 2) * g_dash
467 )
468 return M_sense
470 @staticmethod
471 def A_sensitivity(p0: Vector, p1: Vector, p2: Vector) -> np.array:
472 """Calculates the sensitivity of a cell's area to the points
473 defining the cell analytically.
475 Parameters
476 ----------
477 p0 : Vector
478 The first point defining the cell.
480 p1 : Vector
481 The second point defining the cell.
483 p2 : Vector
484 The third point defining the cell.
486 Returns
487 -------
488 sensitivity : np.array
489 The sensitivity matrix with size m x n. Rows m refer to
490 the vertices, columns n refer to the vertex coordinates.
491 """
492 g = (
493 ((p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)) ** 2
494 + ((p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)) ** 2
495 + ((p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)) ** 2
496 )
498 A_sense = np.empty(9)
499 for col in range(9):
500 if col == 0: # d/dp0.x
501 g_dash = 2 * (
502 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
503 ) * (-(p1.z - p0.z) + (p2.z - p0.z)) + 2 * (
504 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
505 ) * (
506 -(p2.y - p0.y) + (p1.y - p0.y)
507 )
508 elif col == 1: # d/dp0.y
509 g_dash = 2 * (
510 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
511 ) * (-(p2.z - p0.z) + (p1.z - p0.z)) + 2 * (
512 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
513 ) * (
514 -(p1.x - p0.x) + (p2.x - p0.x)
515 )
516 elif col == 2: # d/dp0.z
517 g_dash = 2 * (
518 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
519 ) * (-(p1.y - p0.y) + (p2.y - p0.y)) + 2 * (
520 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
521 ) * (
522 -(p2.x - p0.x) + (p1.x - p0.x)
523 )
524 elif col == 3: # d/dp1.x
525 g_dash = 2 * (
526 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
527 ) * (-(p2.z - p0.z)) + 2 * (
528 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
529 ) * (
530 (p2.y - p0.y)
531 )
532 elif col == 4: # d/dp1.y
533 g_dash = 2 * (
534 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
535 ) * ((p2.z - p0.z)) + 2 * (
536 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
537 ) * (
538 -(p2.x - p0.x)
539 )
540 elif col == 5: # d/dp1.z
541 g_dash = 2 * (
542 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
543 ) * (-(p2.y - p0.y)) + 2 * (
544 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
545 ) * (
546 (p2.x - p0.x)
547 )
548 elif col == 6: # d/dp2.x
549 g_dash = 2 * (
550 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
551 ) * ((p1.z - p0.z)) + 2 * (
552 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
553 ) * (
554 -(p1.y - p0.y)
555 )
556 elif col == 7: # d/dp2.y
557 g_dash = 2 * (
558 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
559 ) * (-(p1.z - p0.z)) + 2 * (
560 (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x)
561 ) * (
562 (p1.x - p0.x)
563 )
564 elif col == 8: # d/dp2.z
565 g_dash = 2 * (
566 (p1.y - p0.y) * (p2.z - p0.z) - (p1.z - p0.z) * (p2.y - p0.y)
567 ) * ((p1.y - p0.y)) + 2 * (
568 (p1.z - p0.z) * (p2.x - p0.x) - (p1.x - p0.x) * (p2.z - p0.z)
569 ) * (
570 -(p1.x - p0.x)
571 )
573 A_sense[col] = 0.5 * (1 / 2) * g ** (-1 / 2) * g_dash
575 return A_sense
577 @staticmethod
578 def c_sensitivity(p0: Vector, p1: Vector, p2: Vector) -> np.array:
579 """Calculates the sensitivity of a cell's centroid to the
580 points defining the cell.
582 Parameters
583 ----------
584 p0 : Vector
585 The first point defining the cell.
587 p1 : Vector
588 The second point defining the cell.
590 p2 : Vector
591 The third point defining the cell.
593 Returns
594 -------
595 sensitivity : np.array
596 The sensitivity matrix with size m x n x p. Rows m refer to
597 the vertices, columns n refer to the vertex coordinates, and
598 slices p refer to the components of the centroid point.
599 """
600 c_sens = np.array(
601 [
602 [1 / 3, 0, 0, 1 / 3, 0, 0, 1 / 3, 0, 0],
603 [0, 1 / 3, 0, 0, 1 / 3, 0, 0, 1 / 3, 0],
604 [0, 0, 1 / 3, 0, 0, 1 / 3, 0, 0, 1 / 3],
605 ]
606 )
607 return c_sens
610class DegenerateCell(Exception):
611 """Exception raised for degenerate cells."""
613 def __init__(self, message="Degenerate cell"):
614 self.message = message
615 super().__init__(self.message)