From 4b6fceb95813e5b40e84477a7c2a3ee164dcca15 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:19:04 +0100 Subject: [PATCH 1/7] Refactor NumericalData class Relocate and refactor with base class and FD specific class --- src/compas_fd/numdata/__init__.py | 33 +++++++++++++++++++ .../{fd => numdata}/fd_numerical_data.py | 8 +++-- src/compas_fd/numdata/numerical_data.py | 26 +++++++++++++++ 3 files changed, 64 insertions(+), 3 deletions(-) create mode 100644 src/compas_fd/numdata/__init__.py rename src/compas_fd/{fd => numdata}/fd_numerical_data.py (93%) create mode 100644 src/compas_fd/numdata/numerical_data.py diff --git a/src/compas_fd/numdata/__init__.py b/src/compas_fd/numdata/__init__.py new file mode 100644 index 00000000..036c517b --- /dev/null +++ b/src/compas_fd/numdata/__init__.py @@ -0,0 +1,33 @@ +""" +****************** +compas_fd.numdata +****************** + +.. currentmodule:: compas_fd.numdata + + +Classes +======= + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + NumericalData + FDNumericalData + +""" +from __future__ import print_function +from __future__ import absolute_import +from __future__ import division + +import compas +from .numerical_data import NumericalData + +if not compas.IPY: + from .fd_numerical_data import FDNumericalData + +__all__ = ['NumericalData'] + +if not compas.IPY: + __all__ += ['FDNumericalData'] diff --git a/src/compas_fd/fd/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py similarity index 93% rename from src/compas_fd/fd/fd_numerical_data.py rename to src/compas_fd/numdata/fd_numerical_data.py index fb061df5..e3fd4689 100644 --- a/src/compas_fd/fd/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -17,11 +17,12 @@ from compas.numerical import connectivity_matrix -from .result import Result +from .numerical_data import NumericalData +from ..result import Result @dataclass -class FDNumericalData: +class FDNumericalData(NumericalData): """Stores numerical data used by the force density algorithms.""" free: int fixed: int @@ -48,7 +49,8 @@ def from_params(cls, fixed: List[int], edges: List[Tuple[int, int]], forcedensities: List[float], - loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None): + loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None + ): """Construct numerical arrays from force density solver input parameters.""" free = list(set(range(len(vertices))) - set(fixed)) xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py new file mode 100644 index 00000000..752cac39 --- /dev/null +++ b/src/compas_fd/numdata/numerical_data.py @@ -0,0 +1,26 @@ +from dataclasses import dataclass +from dataclasses import astuple + +from ..result import Result + + +@dataclass +class NumericalData: + """Storage class for numerical arrays used by solver algorithms.""" + + def __iter__(self): + return iter(astuple(self)) + + @classmethod + def from_params(cls, *args, **kwargs): + """Construct nuerical arrays from algorithm input parameters.""" + raise NotImplementedError + + @classmethod + def from_mesh(cls, mesh): + """Construct numerical arrays from input mesh.""" + raise NotImplementedError + + def to_result(self) -> Result: + """Parse relevant numerical data into a Result object.""" + raise NotImplementedError From 5cbb2be40eb518f975345d022da511f45bf37de2 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 14:19:00 +0100 Subject: [PATCH 2/7] Cleanup --- src/compas_fd/numdata/fd_numerical_data.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index e3fd4689..be33c9ac 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -8,7 +8,6 @@ from nptyping import NDArray from dataclasses import dataclass -from dataclasses import astuple from numpy import asarray from numpy import float64 @@ -40,9 +39,6 @@ class FDNumericalData(NumericalData): tangent_residuals: NDArray[(Any, 3), float64] = None normal_residuals: NDArray[(Any, 1), float64] = None - def __iter__(self): - return iter(astuple(self)) - @classmethod def from_params(cls, vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], From e525695057ee5983e37d1b42fb118192b16a9bdb Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 14:25:42 +0100 Subject: [PATCH 3/7] Add FaceDataMixin FaceData is a mixin for NumericalData. Its responsibility is to keep track of all numerical mesh face data, and lazily compute face data arrays- only when demanded. Its numerical arrays can be used for general purposes: face load computations, face stresses, Natural Force Density solvers, ... --- src/compas_fd/numdata/face_data.py | 123 ++++++++++++++++ src/compas_fd/numdata/fd_numerical_data.py | 159 ++++++++++++++++++++- 2 files changed, 277 insertions(+), 5 deletions(-) create mode 100644 src/compas_fd/numdata/face_data.py diff --git a/src/compas_fd/numdata/face_data.py b/src/compas_fd/numdata/face_data.py new file mode 100644 index 00000000..741f5a25 --- /dev/null +++ b/src/compas_fd/numdata/face_data.py @@ -0,0 +1,123 @@ +from numpy import add +from numpy import cross +from numpy import roll +from numpy import zeros +from numpy.linalg import norm +from scipy.sparse import coo_matrix + +from compas.datastructures import Mesh +from compas.numerical import face_matrix + + +class FaceDataMixin: + """Stores numerical mesh face data. + Use as mixin for NumericalData classes. + """ + + def __init__(self, mesh: Mesh): + self.faces_vertices = mesh + self.reset_face_data() + + def reset_face_data(self): + self._tributary_area_matrix = None + self._face_areas = None + self._face_normals = None + self._face_centroids = None + self._face_cross_products = None + + @property + def faces_vertices(self): + return self._faces_vertices + + @faces_vertices.setter + def faces_vertices(self, mesh: Mesh): + v_i = mesh.vertex_index() + f_i = {fkey: index for index, fkey in enumerate(mesh.faces())} + faces_vertices = [None] * mesh.number_of_faces() + for f in mesh.faces(): + faces_vertices[f_i[f]] = [v_i[v] for v in mesh.face_vertices(f)] + self._faces_vertices = faces_vertices + self._face_matrix = None + + @property + def face_matrix(self): + if self._face_matrix is None: + self._face_matrix = face_matrix(self.faces_vertices, rtype='csr', normalize=True) + return self._face_matrix + + @property + def tributary_area_matrix(self): + """Tributary areas matrix as a sparse (vertices x faces) array. + Entry a_ij holds tributary area of face j for vertex i. + """ + if self._tributary_area_matrix is None: + self._compute_tributary_area_matrix() + return self._tributary_area_matrix + + def _compute_tributary_area_matrix(self): + v_count = self.xyz.shape[0] + f_count = len(self.faces_vertices) + data = [] + rows = [] + cols = [] + + for face, fcp in enumerate(self.face_cross_products): + face_vertices = self.faces_vertices[face] + partial_areas = norm(fcp, axis=1) + tributary_area = partial_areas + roll(partial_areas, -1) + data.extend(tributary_area) + rows.extend(face_vertices) + cols.extend([face] * len(face_vertices)) + + self._tributary_area_matrix = coo_matrix((data, (rows, cols)), + (v_count, f_count)).tocsr() * 0.25 + + @property + def face_areas(self): + if self._face_areas is None: + self._face_areas = self.tributary_area_matrix.sum(axis=0).reshape((-1, 1)) + return self._face_areas + + @property + def face_normals(self): + if self._face_normals is None: + self._compute_face_normals() + return self._face_normals + + def _compute_face_normals(self): + f_count = len(self.faces_vertices) + self._face_normals = zeros((f_count, 3), dtype=float) + for face, fcp in enumerate(self.face_cross_products): + self._face_normals[face, :] = add.reduce(fcp) + self._face_normals /= norm(self._face_normals, axis=1)[:, None] + + @property + def face_centroids(self): + if self._face_centroids is None: + self._face_centroids = self.F.dot(self.xyz) + return self._face_centroids + + @property + def face_cross_products(self): + """Cross products of the consecutive (centroid -> vertex) + vectors for each of the vertices of the face. + """ + if self._face_cross_products is None: + self._compute_face_cross_products() + return self._face_cross_products + + def _compute_face_cross_products(self): + self._face_cross_products = [] + centroids = self.face_centroids + for face, face_vertices in enumerate(self.faces_vertices): + vecs = self.xyz[face_vertices] - centroids[face] + vecs_shift = roll(vecs, -1, axis=0) + self._face_cross_products.append(cross(vecs, vecs_shift)) + + # aliases + F = face_matrix + TA = tributary_area_matrix + fa = face_areas + fn = face_normals + fc = face_centroids + fcp = face_cross_products diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index be33c9ac..11252885 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -14,14 +14,14 @@ from numpy import zeros_like from scipy.sparse import diags +from compas.datastructures import Mesh from compas.numerical import connectivity_matrix from .numerical_data import NumericalData -from ..result import Result +from .face_data import FaceDataMixin -@dataclass -class FDNumericalData(NumericalData): +class FDNumericalData(FaceDataMixin, NumericalData): """Stores numerical data used by the force density algorithms.""" free: int fixed: int @@ -63,10 +63,159 @@ def from_params(cls, return cls(free, fixed, xyz, C, q, Q, p, A, Ai, Af) @classmethod - def from_mesh(cls, mesh): + def from_mesh(cls, mesh: Mesh): """Construct numerical arrays from input mesh.""" - raise NotImplementedError + v_i = mesh.vertex_index() + vertices = asarray(mesh.vertices_attributes('xyz')) + fixed = [v_i[v] for v in mesh.vertices_where({'is_anchor': True})] + edges = [(v_i[u], v_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] + force_densities = asarray([attr['q'] for key, attr in + mesh.edges_where({'_is_edge': True}, True)]) + loads = asarray(mesh.vertices_attributes(('px', 'py', 'pz'))) + + numdata = cls(vertices, fixed, edges, force_densities, loads) + super(FDNumericalData, numdata).__init__(mesh) + numdata.has_face_data = True + return numdata def to_result(self) -> Result: """Parse relevant numerical data into a Result object.""" return Result(self.xyz, self.residuals, self.forces, self.lengths) + + @property + def xyz(self): + return self._xyz + + @xyz.setter + def xyz(self, vertices): + self._xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) + self.reset_xyz() + + def reset_xyz(self): + self._lengths = None + self._residuals = None + if getattr(self, 'has_face_data', False): + FaceDataMixin.reset_face_data(self) + + @property + def fixed(self): + return self._fixed + + @fixed.setter + def fixed(self, fixed): + self._fixed = fixed + self._free = list(set(range(len(self.xyz))) - set(fixed)) + + @property + def free(self): + return self._free + + @property + def edges(self): + return self._edges + + @edges.setter + def edges(self, edges): + self._edges = edges + self._connectivity_matrix = None + self._connectivity_matrix_free = None + self._connectivity_matrix_fixed = None + + @property + def connectivity_matrix(self): + if self._connectivity_matrix is None: + self._connectivity_matrix = connectivity_matrix(self.edges, 'csr') + return self._connectivity_matrix + + @property + def connectivity_matrix_free(self): + if self._connectivity_matrix_free is None: + self._connectivity_matrix_free = self.connectivity_matrix[:, self.free] + return self._connectivity_matrix_free + + @property + def connectivity_matrix_fixed(self): + if self._connectivity_matrix_fixed is None: + self._connectivity_matrix_fixed = self.connectivity_matrix[:, self.fixed] + return self._connectivity_matrix_fixed + + @property + def force_densities(self): + return self._force_densities + + @force_densities.setter + def force_densities(self, force_densities): + self._force_densities = asarray(force_densities, dtype=float64).reshape((-1, 1)) + self._force_densities_matrix = None + self._stifness_matrix = None + self._stifness_matrix_free = None + self._stifness_matrix_fixed = None + self._forces = None + self._residuals = None + + @property + def force_densities_matrix(self): + if self._force_densities_matrix is None: + self._force_densities_matrix = diags( + [self.force_densities.flatten()], [0]) + return self._force_densities_matrix + + @property + def loads(self): + return self._loads + + @loads.setter + def loads(self, loads): + self._loads = (zeros_like(self.xyz) if loads is None else + asarray(loads, dtype=float64).reshape((-1, 3))) + self._residuals = None + + @property + def stifness_matrix(self): + if self._stifness_matrix is None: + self._stifness_matrix = self.C.T.dot(self.Q).dot(self.C) + return self._stifness_matrix + + @property + def stifness_matrix_free(self): + if self._stifness_matrix_free is None: + self._stifness_matrix_free = self.Ci.T.dot(self.Q).dot(self.Ci) + return self._stifness_matrix_free + + @property + def stifness_matrix_fixed(self): + if self._stifness_matrix_fixed is None: + self._stifness_matrix_fixed = self.Ci.T.dot(self.Q).dot(self.Cf) + return self._stifness_matrix_fixed + + @property + def forces(self): + if self._forces is None: + self._forces = self.q * self.ls + return self._forces + + @property + def lengths(self): + if self._lengths is None: + self._lengths = normrow(self.C.dot(self.xyz)) + return self._lengths + + @property + def residuals(self): + if self._residuals is None: + self._residuals = self.p - self.D.dot(self.xyz) + return self._residuals + + # aliases + C = connectivity_matrix + Ci = connectivity_matrix_free + Cf = connectivity_matrix_fixed + q = force_densities + Q = force_densities_matrix + p = loads + D = stifness_matrix + Di = stifness_matrix_free + Df = stifness_matrix_fixed + f = forces + ls = lengths + r = residuals From 850a0e4d928a0040361ba7376e2a9c2665d7fc0b Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 14:31:01 +0100 Subject: [PATCH 4/7] Lint --- src/compas_fd/numdata/fd_numerical_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index 11252885..45fa3958 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -7,16 +7,16 @@ from typing_extensions import Annotated from nptyping import NDArray -from dataclasses import dataclass - from numpy import asarray from numpy import float64 +from numpy import normrow from numpy import zeros_like from scipy.sparse import diags from compas.datastructures import Mesh from compas.numerical import connectivity_matrix +from compas_fd.result import Result from .numerical_data import NumericalData from .face_data import FaceDataMixin From a11c6f39730016eaa4f789fe0af525951741d7b7 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 14:35:21 +0100 Subject: [PATCH 5/7] Refactor FDNumericalData Change FDNumericalData from dataclass into regular class with lazy calculation of property variables. Implicit recalculation of dependent variables. Add shorthand aliases to numerical arrays. --- src/compas_fd/numdata/fd_numerical_data.py | 51 ++++++---------------- 1 file changed, 13 insertions(+), 38 deletions(-) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index 45fa3958..6e8ae38b 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -9,12 +9,12 @@ from numpy import asarray from numpy import float64 -from numpy import normrow from numpy import zeros_like from scipy.sparse import diags from compas.datastructures import Mesh from compas.numerical import connectivity_matrix +from compas.numerical import normrow from compas_fd.result import Result from .numerical_data import NumericalData @@ -23,44 +23,19 @@ class FDNumericalData(FaceDataMixin, NumericalData): """Stores numerical data used by the force density algorithms.""" - free: int - fixed: int - xyz: NDArray[(Any, 3), float64] - C: NDArray[(Any, Any), int] - q: NDArray[(Any, 1), float64] - Q: NDArray[(Any, Any), float64] - p: NDArray[(Any, 1), float64] - A: NDArray[(Any, Any), float64] - Ai: NDArray[(Any, Any), float64] - Af: NDArray[(Any, Any), float64] - forces: NDArray[(Any, 1), float64] = None - lengths: NDArray[(Any, 1), float64] = None - residuals: NDArray[(Any, 3), float64] = None - tangent_residuals: NDArray[(Any, 3), float64] = None - normal_residuals: NDArray[(Any, 1), float64] = None - @classmethod - def from_params(cls, - vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], - fixed: List[int], - edges: List[Tuple[int, int]], - forcedensities: List[float], - loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None - ): - """Construct numerical arrays from force density solver input parameters.""" - free = list(set(range(len(vertices))) - set(fixed)) - xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) - C = connectivity_matrix(edges, 'csr') - Ci = C[:, free] - Cf = C[:, fixed] - q = asarray(forcedensities, dtype=float64).reshape((-1, 1)) - Q = diags([q.flatten()], [0]) - p = (zeros_like(xyz) if loads is None else - asarray(loads, dtype=float64).reshape((-1, 3))) - A = C.T.dot(Q).dot(C) - Ai = Ci.T.dot(Q).dot(Ci) - Af = Ci.T.dot(Q).dot(Cf) - return cls(free, fixed, xyz, C, q, Q, p, A, Ai, Af) + def __init__(self, + vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], + fixed: List[int], + edges: List[Tuple[int, int]], + force_densities: List[float], + loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None + ): + self.xyz = vertices + self.fixed = fixed + self.edges = edges + self.force_densities = force_densities + self.loads = loads @classmethod def from_mesh(cls, mesh: Mesh): From 4d7f89d0557a70519abf3ebad468ff1c92880572 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 18:39:55 +0100 Subject: [PATCH 6/7] Refactor NumericalData Formalize lazy evaluation into a decorator. Add matrices to FDNumericalData. --- src/compas_fd/numdata/face_data.py | 79 +++++++------- src/compas_fd/numdata/fd_numerical_data.py | 120 ++++++++++++--------- src/compas_fd/numdata/numerical_data.py | 35 ++++-- 3 files changed, 127 insertions(+), 107 deletions(-) diff --git a/src/compas_fd/numdata/face_data.py b/src/compas_fd/numdata/face_data.py index 741f5a25..d5135cd5 100644 --- a/src/compas_fd/numdata/face_data.py +++ b/src/compas_fd/numdata/face_data.py @@ -1,5 +1,11 @@ +from typing import Any +from typing import List +from typing import Tuple +from nptyping import NDArray + from numpy import add from numpy import cross +from numpy import float64 from numpy import roll from numpy import zeros from numpy.linalg import norm @@ -7,6 +13,7 @@ from compas.datastructures import Mesh from compas.numerical import face_matrix +from .numerical_data import lazy_eval class FaceDataMixin: @@ -14,11 +21,11 @@ class FaceDataMixin: Use as mixin for NumericalData classes. """ - def __init__(self, mesh: Mesh): + def __init__(self, mesh: Mesh) -> None: self.faces_vertices = mesh self.reset_face_data() - def reset_face_data(self): + def reset_face_data(self) -> None: self._tributary_area_matrix = None self._face_areas = None self._face_normals = None @@ -26,11 +33,11 @@ def reset_face_data(self): self._face_cross_products = None @property - def faces_vertices(self): + def faces_vertices(self) -> List[Tuple[int]]: return self._faces_vertices @faces_vertices.setter - def faces_vertices(self, mesh: Mesh): + def faces_vertices(self, mesh: Mesh) -> None: v_i = mesh.vertex_index() f_i = {fkey: index for index, fkey in enumerate(mesh.faces())} faces_vertices = [None] * mesh.number_of_faces() @@ -40,21 +47,16 @@ def faces_vertices(self, mesh: Mesh): self._face_matrix = None @property - def face_matrix(self): - if self._face_matrix is None: - self._face_matrix = face_matrix(self.faces_vertices, rtype='csr', normalize=True) - return self._face_matrix + @lazy_eval + def face_matrix(self) -> NDArray[(Any, Any), float64]: + return face_matrix(self.faces_vertices, rtype='csr', normalize=True) @property - def tributary_area_matrix(self): + @lazy_eval + def tributary_area_matrix(self) -> NDArray[(Any, Any), float64]: """Tributary areas matrix as a sparse (vertices x faces) array. Entry a_ij holds tributary area of face j for vertex i. """ - if self._tributary_area_matrix is None: - self._compute_tributary_area_matrix() - return self._tributary_area_matrix - - def _compute_tributary_area_matrix(self): v_count = self.xyz.shape[0] f_count = len(self.faces_vertices) data = [] @@ -64,55 +66,46 @@ def _compute_tributary_area_matrix(self): for face, fcp in enumerate(self.face_cross_products): face_vertices = self.faces_vertices[face] partial_areas = norm(fcp, axis=1) - tributary_area = partial_areas + roll(partial_areas, -1) - data.extend(tributary_area) + tributary_areas = partial_areas + roll(partial_areas, -1) + data.extend(tributary_areas) rows.extend(face_vertices) cols.extend([face] * len(face_vertices)) - self._tributary_area_matrix = coo_matrix((data, (rows, cols)), - (v_count, f_count)).tocsr() * 0.25 + return coo_matrix((data, (rows, cols)), (v_count, f_count)).tocsr() * 0.25 @property - def face_areas(self): - if self._face_areas is None: - self._face_areas = self.tributary_area_matrix.sum(axis=0).reshape((-1, 1)) - return self._face_areas + @lazy_eval + def face_areas(self) -> NDArray[(Any, 1), float64]: + return self.tributary_area_matrix.sum(axis=0).reshape((-1, 1)) @property - def face_normals(self): - if self._face_normals is None: - self._compute_face_normals() - return self._face_normals - - def _compute_face_normals(self): + @lazy_eval + def face_normals(self) -> NDArray[(Any, 1), float64]: f_count = len(self.faces_vertices) - self._face_normals = zeros((f_count, 3), dtype=float) + face_normals = zeros((f_count, 3), dtype=float) for face, fcp in enumerate(self.face_cross_products): - self._face_normals[face, :] = add.reduce(fcp) - self._face_normals /= norm(self._face_normals, axis=1)[:, None] + face_normals[face, :] = add.reduce(fcp) + face_normals /= norm(face_normals, axis=1)[:, None] + return face_normals @property - def face_centroids(self): - if self._face_centroids is None: - self._face_centroids = self.F.dot(self.xyz) - return self._face_centroids + @lazy_eval + def face_centroids(self) -> NDArray[(Any, 1), float64]: + return self.F.dot(self.xyz) @property - def face_cross_products(self): + @lazy_eval + def face_cross_products(self) -> List[NDArray[(Any, 3), float64]]: """Cross products of the consecutive (centroid -> vertex) vectors for each of the vertices of the face. """ - if self._face_cross_products is None: - self._compute_face_cross_products() - return self._face_cross_products - - def _compute_face_cross_products(self): - self._face_cross_products = [] + face_cross_products = [] centroids = self.face_centroids for face, face_vertices in enumerate(self.faces_vertices): vecs = self.xyz[face_vertices] - centroids[face] vecs_shift = roll(vecs, -1, axis=0) - self._face_cross_products.append(cross(vecs, vecs_shift)) + face_cross_products.append(cross(vecs, vecs_shift)) + return face_cross_products # aliases F = face_matrix diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index 6e8ae38b..58c4d63f 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -11,6 +11,7 @@ from numpy import float64 from numpy import zeros_like from scipy.sparse import diags +from scipy.sparse import vstack as svstack from compas.datastructures import Mesh from compas.numerical import connectivity_matrix @@ -18,6 +19,7 @@ from compas_fd.result import Result from .numerical_data import NumericalData +from .numerical_data import lazy_eval from .face_data import FaceDataMixin @@ -58,68 +60,84 @@ def to_result(self) -> Result: return Result(self.xyz, self.residuals, self.forces, self.lengths) @property - def xyz(self): + def xyz(self) -> NDArray[(Any, 3), float64]: return self._xyz @xyz.setter - def xyz(self, vertices): + def xyz(self, vertices) -> None: self._xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) self.reset_xyz() - def reset_xyz(self): + def reset_xyz(self) -> None: self._lengths = None self._residuals = None + self._coordinate_difference_matrices = None + self._equilibrium_matrix = None if getattr(self, 'has_face_data', False): FaceDataMixin.reset_face_data(self) @property - def fixed(self): + def fixed(self) -> List[int]: return self._fixed @fixed.setter - def fixed(self, fixed): + def fixed(self, fixed) -> None: self._fixed = fixed self._free = list(set(range(len(self.xyz))) - set(fixed)) @property - def free(self): + def free(self) -> List[int]: return self._free @property - def edges(self): + def edges(self) -> List[Tuple[int, int]]: return self._edges @edges.setter - def edges(self, edges): + def edges(self, edges: List[List[int]]) -> None: self._edges = edges self._connectivity_matrix = None self._connectivity_matrix_free = None self._connectivity_matrix_fixed = None @property - def connectivity_matrix(self): - if self._connectivity_matrix is None: - self._connectivity_matrix = connectivity_matrix(self.edges, 'csr') - return self._connectivity_matrix + @lazy_eval + def connectivity_matrix(self) -> NDArray[(Any, Any), float64]: + return connectivity_matrix(self.edges, 'csr') @property - def connectivity_matrix_free(self): - if self._connectivity_matrix_free is None: - self._connectivity_matrix_free = self.connectivity_matrix[:, self.free] - return self._connectivity_matrix_free + @lazy_eval + def connectivity_matrix_free(self) -> NDArray[(Any, Any), float64]: + return self.connectivity_matrix[:, self.free] @property - def connectivity_matrix_fixed(self): - if self._connectivity_matrix_fixed is None: - self._connectivity_matrix_fixed = self.connectivity_matrix[:, self.fixed] - return self._connectivity_matrix_fixed + @lazy_eval + def connectivity_matrix_fixed(self) -> NDArray[(Any, Any), float64]: + return self.connectivity_matrix[:, self.fixed] @property - def force_densities(self): + @lazy_eval + def coordinate_difference_matrices(self) -> Tuple[NDArray[(Any, Any), float64]]: + uvw = self.C.dot(self.xyz) + U = diags([uvw[:, 0].flatten()], [0]) + V = diags([uvw[:, 1].flatten()], [0]) + W = diags([uvw[:, 2].flatten()], [0]) + return (U, V, W) + + @property + @lazy_eval + def equilibrium_matrix(self) -> NDArray[(Any, Any), float64]: + U, V, W = self.coordinate_difference_matrices + return svstack((self.C.T.dot(U), self.C.T.dot(V), self.C.T.dot(W))) + + @property + def force_densities(self) -> NDArray[(Any, 1), float64]: return self._force_densities @force_densities.setter - def force_densities(self, force_densities): + def force_densities(self, + force_densities: Union[List[float], NDArray[(Any, 3), float64]] + ) -> None: self._force_densities = asarray(force_densities, dtype=float64).reshape((-1, 1)) self._force_densities_matrix = None self._stifness_matrix = None @@ -129,62 +147,58 @@ def force_densities(self, force_densities): self._residuals = None @property - def force_densities_matrix(self): - if self._force_densities_matrix is None: - self._force_densities_matrix = diags( - [self.force_densities.flatten()], [0]) - return self._force_densities_matrix + @lazy_eval + def force_densities_matrix(self) -> NDArray[(Any, Any), float64]: + return diags([self.force_densities.flatten()], [0]) @property - def loads(self): + def loads(self) -> NDArray[(Any, 3), float64]: return self._loads @loads.setter - def loads(self, loads): + def loads(self, + loads: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]] + ) -> None: self._loads = (zeros_like(self.xyz) if loads is None else asarray(loads, dtype=float64).reshape((-1, 3))) self._residuals = None @property - def stifness_matrix(self): - if self._stifness_matrix is None: - self._stifness_matrix = self.C.T.dot(self.Q).dot(self.C) - return self._stifness_matrix + @lazy_eval + def stifness_matrix(self) -> NDArray[(Any, Any), float64]: + return self.C.T.dot(self.Q).dot(self.C) @property - def stifness_matrix_free(self): - if self._stifness_matrix_free is None: - self._stifness_matrix_free = self.Ci.T.dot(self.Q).dot(self.Ci) - return self._stifness_matrix_free + @lazy_eval + def stifness_matrix_free(self) -> NDArray[(Any, Any), float64]: + return self.Ci.T.dot(self.Q).dot(self.Ci) @property - def stifness_matrix_fixed(self): - if self._stifness_matrix_fixed is None: - self._stifness_matrix_fixed = self.Ci.T.dot(self.Q).dot(self.Cf) - return self._stifness_matrix_fixed + @lazy_eval + def stifness_matrix_fixed(self) -> NDArray[(Any, Any), float64]: + return self.Ci.T.dot(self.Q).dot(self.Cf) @property - def forces(self): - if self._forces is None: - self._forces = self.q * self.ls - return self._forces + @lazy_eval + def forces(self) -> NDArray[(Any, 1), float64]: + return self.q * self.ls @property - def lengths(self): - if self._lengths is None: - self._lengths = normrow(self.C.dot(self.xyz)) - return self._lengths + @lazy_eval + def lengths(self) -> NDArray[(Any, 1), float64]: + return normrow(self.C.dot(self.xyz)) @property - def residuals(self): - if self._residuals is None: - self._residuals = self.p - self.D.dot(self.xyz) - return self._residuals + @lazy_eval + def residuals(self) -> NDArray[(Any, 3), float64]: + return self.p - self.D.dot(self.xyz) # aliases C = connectivity_matrix Ci = connectivity_matrix_free Cf = connectivity_matrix_fixed + UVW = coordinate_difference_matrices + E = equilibrium_matrix q = force_densities Q = force_densities_matrix p = loads diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py index 752cac39..49622794 100644 --- a/src/compas_fd/numdata/numerical_data.py +++ b/src/compas_fd/numdata/numerical_data.py @@ -1,21 +1,12 @@ -from dataclasses import dataclass -from dataclasses import astuple +from typing import Callable +from functools import wraps from ..result import Result -@dataclass class NumericalData: """Storage class for numerical arrays used by solver algorithms.""" - def __iter__(self): - return iter(astuple(self)) - - @classmethod - def from_params(cls, *args, **kwargs): - """Construct nuerical arrays from algorithm input parameters.""" - raise NotImplementedError - @classmethod def from_mesh(cls, mesh): """Construct numerical arrays from input mesh.""" @@ -24,3 +15,25 @@ def from_mesh(cls, mesh): def to_result(self) -> Result: """Parse relevant numerical data into a Result object.""" raise NotImplementedError + + +def lazy_eval(getter: Callable) -> Callable: + """Decorator for lazy evaluation on property getters. + + A private attribute name is assumed as ('_' + getter name). + The computation in the getter function is done only if: + - the private attribute does not exist in the instance dictionary. + - or the private attribute is None. + After computation by the getter function, the private attribute is stored + in the instance attribute dictionary. + """ + private_attr = '_' + getter.__name__ + + @wraps(getter) + def wrapper(self, *args, **kwargs): + attr = getattr(self, private_attr, None) + if attr is None: + attr = getter(self, *args, **kwargs) + self.__dict__[private_attr] = attr + return attr + return wrapper From 98038f8d6b8a912890d5edc0c7f640ca2d230982 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 21:38:14 +0100 Subject: [PATCH 7/7] Change result location Change result location --- src/compas_fd/numdata/numerical_data.py | 2 +- src/compas_fd/result/__init__.py | 25 +++++++++++++++++++++++++ src/compas_fd/{fd => result}/result.py | 0 3 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 src/compas_fd/result/__init__.py rename src/compas_fd/{fd => result}/result.py (100%) diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py index 49622794..71499790 100644 --- a/src/compas_fd/numdata/numerical_data.py +++ b/src/compas_fd/numdata/numerical_data.py @@ -1,7 +1,7 @@ from typing import Callable from functools import wraps -from ..result import Result +from compas_fd.result import Result class NumericalData: diff --git a/src/compas_fd/result/__init__.py b/src/compas_fd/result/__init__.py new file mode 100644 index 00000000..fa376bd3 --- /dev/null +++ b/src/compas_fd/result/__init__.py @@ -0,0 +1,25 @@ +""" +****************** +compas_fd.result +****************** + +.. currentmodule:: compas_fd.result + + +Classes +======= + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + Result + +""" +from __future__ import print_function +from __future__ import absolute_import +from __future__ import division + +from .result import Result + +__all__ = ['Result'] diff --git a/src/compas_fd/fd/result.py b/src/compas_fd/result/result.py similarity index 100% rename from src/compas_fd/fd/result.py rename to src/compas_fd/result/result.py