From 0441bac8b0c1d0159c23bea69aeb3e5513fad764 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Wed, 20 Oct 2021 10:15:00 +0200 Subject: [PATCH 01/13] Add LoadCalculator Updated load processor for self-weight, wind (normal), and snow (projected) loads. Lazy calculation, only loaded faces are processed. --- src/compas_fd/datastructures/cablemesh.py | 8 +- src/compas_fd/loads/__init__.py | 6 +- src/compas_fd/loads/load_calculator.py | 210 ++++++++++++++++++++++ 3 files changed, 220 insertions(+), 4 deletions(-) create mode 100644 src/compas_fd/loads/load_calculator.py diff --git a/src/compas_fd/datastructures/cablemesh.py b/src/compas_fd/datastructures/cablemesh.py index a114b3a8..05fb1c9f 100644 --- a/src/compas_fd/datastructures/cablemesh.py +++ b/src/compas_fd/datastructures/cablemesh.py @@ -36,6 +36,9 @@ def __init__(self): 'is_fixed': False, 'constraint': None, 'param': None, + '_px': 0.0, + '_py': 0.0, + '_pz': 0.0, '_rx': 0.0, '_ry': 0.0, '_rz': 0.0, @@ -48,5 +51,8 @@ def __init__(self): '_r': 0.0, }) self.default_face_attributes.update({ - 'is_loaded': True, + 'is_loaded': False, + 't': 0.0, + 'wind': 0.0, + 'snow': 0.0 }) diff --git a/src/compas_fd/loads/__init__.py b/src/compas_fd/loads/__init__.py index c8e1cc9a..ffdeb8a4 100644 --- a/src/compas_fd/loads/__init__.py +++ b/src/compas_fd/loads/__init__.py @@ -13,7 +13,7 @@ :toctree: generated/ :nosignatures: - SelfweightCalculator + LoadCalculator """ from __future__ import print_function @@ -23,9 +23,9 @@ import compas if not compas.IPY: - from .selfweight import SelfweightCalculator + from .load_calculator import LoadCalculator __all__ = [] if not compas.IPY: - __all__ += ['SelfweightCalculator'] + __all__ += ['LoadCalculator'] diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py new file mode 100644 index 00000000..3c6be7b3 --- /dev/null +++ b/src/compas_fd/loads/load_calculator.py @@ -0,0 +1,210 @@ +from numpy import add +from numpy import array +from numpy import asarray +from numpy import average +from numpy import cross +from numpy import empty +from numpy import float64 +from numpy import roll + +from numpy.linalg import norm +from scipy.sparse import coo_matrix + +from compas.datastructures import Mesh + +from typing import Any +from typing import List +from typing import Sequence +from typing import Union +from typing_extensions import Annotated +from nptyping import NDArray + + +class LoadCalculator: + """Represents assembly of all input loads into a global vertex load matrix. + + By calling an instance of LoadCalculator with vertex coordinates, + the global load matrix is calculated for the current geometry. + This calculation is lazy, face data (tributary areas or normals) + is only calculated if face loads or oriented loads are applied. + + Parameters + ---------- + mesh : :class:`Mesh` + Instance of Mesh datastructure. + + Attributes + ---------- + loads_mat : ndarray of float + global vertex load matrix as (v x 3) array. + face_normals : ndarray of float + Face normals for current geometry as (f x 3) array. + trib_areas : ndarray of float + Tributary areas matrix as sparse (v x f) array. + Entry a_ij holds tributary area of face j for vertex i. + + Examples + -------- + >>> import compas + >>> from compas.datastructures import Mesh + >>> from compas_fd import LoadCalculator + >>> mesh = Mesh.from_obj(compas.get('hypar.obj')) + >>> dva = {'px': .0, 'py': .0, 'pz': .1} + >>> dfa = {'density': 2.4, 't': .1, + 'is_loaded': True, 'wind': 1.0} + >>> mesh.update_default_vertex_attributes(dva) + >>> mesh.update_default_face_attributes(dfa) + >>> lc = LoadCalculator(mesh) + >>> xyz = [mesh.vertex_coordinates(v) + for v in mesh.vertices()] + >>> load_mat = lc(xyz) + >>> lc.update_mesh() + """ + + # attribute names in mesh data + PX, PY, PZ = 'px', 'py', 'pz' # input vertex loads + _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads + THICKNESS = 't' + RHO = 'density' + SW = 'is_loaded' + NORMAL = 'wind' + PROJECT = 'snow' + + def __init__(self, mesh: Mesh) -> None: + self.mesh = mesh + self._set_mesh_maps() + self._preprocess_loads() + self._set_load_flags() + + def __call__(self, xyz: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], + all_faces: bool = False) -> NDArray[(Any, 3), float64]: + """Assemble global load matrix for current geometry. + + Parameters + ---------- + xyz : ndarray of float or list of lists + Coordinates of the vertices as (v x 3) array. + all_faces : bool + Switch to calculate all face tributary areas + and normals, regardless of the applied loads. + Default is ``False``. + + Returns + ------- + ndarray of float + Resultant global vertex loads as (v x 3) array. + """ + if not self._load_flags['any_face_load'] and not all_faces: + return self._vertex_loads + self.loads_mat = array(self._vertex_loads, copy=True) + self._process_faces(asarray(xyz).reshape(-1, 3), all_faces) + self._add_face_loads() + return self.loads_mat + + def update_mesh(self) -> None: + """Set the resultant vertex loads in mesh data by + using the latest calculated internal load matrix.""" + try: + lmat = self.loads_mat + except AttributeError: + lmat = self._vertex_loads + for v, vkey in enumerate(self.mesh.vertices()): + self.mesh.vertex_attributes(vkey, [self._PX, self._PY, self._PZ], + [lmat[v, 0], lmat[v, 1], lmat[v, 2]]) + + def _set_mesh_maps(self) -> None: + """Set mesh counts and maps.""" + self._v_count = self.mesh.number_of_vertices() + self._f_count = self.mesh.number_of_faces() + self._v_id = self.mesh.key_index() + self.face_normals = None + self.trib_area = None + + def _preprocess_loads(self) -> None: + """Preprocess loads into arrays.""" + self._vertex_loads = asarray(self.mesh.vertices_attributes( + names=(self.PX, self.PY, self.PZ)) + ).reshape((self._v_count, 3)) + self._weight = asarray([-t * rho * w for t, rho, w in + self.mesh.faces_attributes( + [self.THICKNESS, self.RHO, self.SW])] + ).reshape((self._f_count, 1)) + self._normal = asarray(self.mesh.faces_attribute(self.NORMAL) + ).reshape((self._f_count, 1)) + self._project = asarray(self.mesh.faces_attribute(self.PROJECT) + ).reshape((self._f_count, 1)) + + def _set_load_flags(self) -> None: + """Set flags for which load types are applied.""" + bool_weight = array(self._weight, dtype=bool) + bool_normal = array(self._normal, dtype=bool) + bool_project = array(self._project, dtype=bool) + self._oriented_loads = bool_normal + bool_project + self._face_loads = self._oriented_loads + bool_weight + self._load_flags = {'any_face_load': self._face_loads.any(), + 'any_weight': bool_weight.any(), + 'any_normal': bool_normal.any(), + 'any_projected': bool_project.any()} + + def _add_face_loads(self) -> None: + """Add active face loads to the global load matrix.""" + if self._load_flags['any_weight']: + self._add_weight() + if self._load_flags['any_normal']: + self._add_normal_loads() + if self._load_flags['any_projected']: + self._add_projected_loads() + + def _process_faces(self, xyz: NDArray[(Any, 3), float64], + all_faces: bool) -> None: + """Calculate normals and tributary areas per face, + for the vertex coordinates of the current geometry.""" + self.face_normals = empty((self._f_count, 3), dtype=float) + trib_data = []; trib_rows = []; trib_cols = [] # noqa + + for f, fkey in enumerate(self.mesh.faces()): + if not (self._face_loads[f] or all_faces): + continue + vts = [self._v_id[v] for v in + self.mesh.face_vertices(fkey)] + f_xyz = xyz[vts, :] + cps = self.__face_cross_products(f_xyz) + # face tributary areas per vertex + part_areas = 0.25 * norm(cps, axis=1) + trib_data.extend(part_areas + roll(part_areas, -1)) + trib_rows.extend(vts) + trib_cols.extend([f] * len(vts)) + # face normal vector + if self._oriented_loads[f] or all_faces: + self.face_normals[f, :] = add.reduce(cps) + + self.face_normals /= norm(self.face_normals, axis=1)[:, None] + self.trib_areas = coo_matrix((trib_data, (trib_rows, trib_cols)), + (self._v_count, self._f_count)).tocsr() + + def __face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] + ) -> NDArray[(Any, 3), float64]: + """Utility to calculate cross products of the consecutive + (centroid -> vertex) vectors for each of the vertices of the face.""" + centroid = average(f_xyz, axis=0) + v = f_xyz - centroid + vr = roll(v, -1, axis=0) + return cross(v, vr) + + def _add_weight(self) -> None: + """Convert all face self-weights into global vertex loads + and add them into the global load matrix.""" + self.loads_mat[:, 2:] += self.trib_areas.dot(self._weight) + + def _add_normal_loads(self) -> None: + """Convert all normal face loads into global vertex loads + and add them into the global load matrix.""" + normal_loads = self.face_normals * self._normal + self.loads_mat += self.trib_areas.dot(normal_loads) + + def _add_projected_loads(self) -> None: + """Convert all Z-projected face loads into global vertex loads + and add them into the global load matrix.""" + z = asarray([0, 0, -1]).reshape((3, 1)) + proj_loads = self.face_normals.dot(z) * self._project + self.loads_mat[:, 2:] += self.trib_areas.dot(proj_loads) From e084ce9485ea4838cc1c2307c59e927de00a9968 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Wed, 20 Oct 2021 10:27:58 +0200 Subject: [PATCH 02/13] Update test files --- src/compas_fd/fd/mesh_fd_numpy.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/compas_fd/fd/mesh_fd_numpy.py b/src/compas_fd/fd/mesh_fd_numpy.py index 316fad97..8b8fc53c 100644 --- a/src/compas_fd/fd/mesh_fd_numpy.py +++ b/src/compas_fd/fd/mesh_fd_numpy.py @@ -2,7 +2,7 @@ from numpy import float64 import compas_fd -from compas_fd.loads import SelfweightCalculator +from compas_fd.loads import LoadCalculator from .fd_numpy import fd_numpy @@ -29,9 +29,8 @@ def mesh_fd_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.data p = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) edges = [(k_i[u], k_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] q = array([attr['q'] for key, attr in mesh.edges_where({'_is_edge': True}, True)], dtype=float64).reshape((-1, 1)) - density = mesh.attributes['density'] - calculate_sw = SelfweightCalculator(mesh, density=density) - p[:, 2] -= calculate_sw(xyz)[:, 0] + lc = LoadCalculator(mesh) + p[:, 2] = lc(xyz)[:, 2] result = fd_numpy(vertices=xyz, fixed=fixed, edges=edges, forcedensities=q, loads=p) From e7fddb19f82056ce7216898f5bfb263da63c80c7 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 25 Oct 2021 10:12:02 +0200 Subject: [PATCH 03/13] Density as global variable --- src/compas_fd/loads/load_calculator.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 3c6be7b3..a375c001 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -50,8 +50,8 @@ class LoadCalculator: >>> from compas_fd import LoadCalculator >>> mesh = Mesh.from_obj(compas.get('hypar.obj')) >>> dva = {'px': .0, 'py': .0, 'pz': .1} - >>> dfa = {'density': 2.4, 't': .1, - 'is_loaded': True, 'wind': 1.0} + >>> dfa = {'t': .1, 'is_loaded': True, 'wind': 1.0} + >>> mesh.attributes.update({'density': 22.0}) >>> mesh.update_default_vertex_attributes(dva) >>> mesh.update_default_face_attributes(dfa) >>> lc = LoadCalculator(mesh) @@ -125,9 +125,10 @@ def _preprocess_loads(self) -> None: self._vertex_loads = asarray(self.mesh.vertices_attributes( names=(self.PX, self.PY, self.PZ)) ).reshape((self._v_count, 3)) - self._weight = asarray([-t * rho * w for t, rho, w in + rho = self.mesh.attributes[self.RHO] + self._weight = asarray([-t * w * rho for t, w in self.mesh.faces_attributes( - [self.THICKNESS, self.RHO, self.SW])] + [self.THICKNESS, self.SW])] ).reshape((self._f_count, 1)) self._normal = asarray(self.mesh.faces_attribute(self.NORMAL) ).reshape((self._f_count, 1)) From b80cc09c09494914970e65258cff4a0429b8a77f Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 25 Oct 2021 11:33:24 +0200 Subject: [PATCH 04/13] Update naming for clarity --- src/compas_fd/loads/load_calculator.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index a375c001..962eb99b 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -64,8 +64,8 @@ class LoadCalculator: # attribute names in mesh data PX, PY, PZ = 'px', 'py', 'pz' # input vertex loads _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads - THICKNESS = 't' RHO = 'density' + THICKNESS = 't' SW = 'is_loaded' NORMAL = 'wind' PROJECT = 'snow' @@ -127,19 +127,19 @@ def _preprocess_loads(self) -> None: ).reshape((self._v_count, 3)) rho = self.mesh.attributes[self.RHO] self._weight = asarray([-t * w * rho for t, w in - self.mesh.faces_attributes( - [self.THICKNESS, self.SW])] - ).reshape((self._f_count, 1)) - self._normal = asarray(self.mesh.faces_attribute(self.NORMAL) + self.mesh.faces_attributes([self.THICKNESS, self.SW])] ).reshape((self._f_count, 1)) - self._project = asarray(self.mesh.faces_attribute(self.PROJECT) - ).reshape((self._f_count, 1)) + self._normal_load = asarray(self.mesh.faces_attribute(self.NORMAL) + ).reshape((self._f_count, 1)) + self._project_load = asarray(self.mesh.faces_attribute(self.PROJECT) + ).reshape((self._f_count, 1)) + print(self._weight) def _set_load_flags(self) -> None: """Set flags for which load types are applied.""" bool_weight = array(self._weight, dtype=bool) - bool_normal = array(self._normal, dtype=bool) - bool_project = array(self._project, dtype=bool) + bool_normal = array(self._normal_load, dtype=bool) + bool_project = array(self._project_load, dtype=bool) self._oriented_loads = bool_normal + bool_project self._face_loads = self._oriented_loads + bool_weight self._load_flags = {'any_face_load': self._face_loads.any(), @@ -200,12 +200,12 @@ def _add_weight(self) -> None: def _add_normal_loads(self) -> None: """Convert all normal face loads into global vertex loads and add them into the global load matrix.""" - normal_loads = self.face_normals * self._normal + normal_loads = self.face_normals * self._normal_load self.loads_mat += self.trib_areas.dot(normal_loads) def _add_projected_loads(self) -> None: """Convert all Z-projected face loads into global vertex loads and add them into the global load matrix.""" z = asarray([0, 0, -1]).reshape((3, 1)) - proj_loads = self.face_normals.dot(z) * self._project + proj_loads = self.face_normals.dot(z) * self._project_load self.loads_mat[:, 2:] += self.trib_areas.dot(proj_loads) From ef5c965f5a1457b89ef07704f648b321d5c2339b Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 25 Oct 2021 11:34:00 +0200 Subject: [PATCH 05/13] Cleanup --- src/compas_fd/loads/load_calculator.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 962eb99b..2a579504 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -133,7 +133,6 @@ def _preprocess_loads(self) -> None: ).reshape((self._f_count, 1)) self._project_load = asarray(self.mesh.faces_attribute(self.PROJECT) ).reshape((self._f_count, 1)) - print(self._weight) def _set_load_flags(self) -> None: """Set flags for which load types are applied.""" From d28473973ea90f91ac66d6d0b527b8f14e8b708e Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 25 Oct 2021 16:40:57 +0200 Subject: [PATCH 06/13] Update naming matrix and other naming updated for clarity --- src/compas_fd/loads/load_calculator.py | 74 +++++++++++++------------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 2a579504..92971076 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -35,11 +35,11 @@ class LoadCalculator: Attributes ---------- - loads_mat : ndarray of float + LM : ndarray of float global vertex load matrix as (v x 3) array. - face_normals : ndarray of float + FN : ndarray of float Face normals for current geometry as (f x 3) array. - trib_areas : ndarray of float + TA : ndarray of float Tributary areas matrix as sparse (v x f) array. Entry a_ij holds tributary area of face j for vertex i. @@ -66,7 +66,7 @@ class LoadCalculator: _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads RHO = 'density' THICKNESS = 't' - SW = 'is_loaded' + HAS_SW = 'is_loaded' NORMAL = 'wind' PROJECT = 'snow' @@ -77,14 +77,14 @@ def __init__(self, mesh: Mesh) -> None: self._set_load_flags() def __call__(self, xyz: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], - all_faces: bool = False) -> NDArray[(Any, 3), float64]: + process_all_faces: bool = False) -> NDArray[(Any, 3), float64]: """Assemble global load matrix for current geometry. Parameters ---------- xyz : ndarray of float or list of lists Coordinates of the vertices as (v x 3) array. - all_faces : bool + process_all_faces : bool Switch to calculate all face tributary areas and normals, regardless of the applied loads. Default is ``False``. @@ -94,31 +94,31 @@ def __call__(self, xyz: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, ndarray of float Resultant global vertex loads as (v x 3) array. """ - if not self._load_flags['any_face_load'] and not all_faces: + if not self._load_flags['any_face_load'] and not process_all_faces: return self._vertex_loads - self.loads_mat = array(self._vertex_loads, copy=True) - self._process_faces(asarray(xyz).reshape(-1, 3), all_faces) + self.LM = array(self._vertex_loads, copy=True) + self._process_faces(asarray(xyz).reshape(-1, 3), process_all_faces) self._add_face_loads() - return self.loads_mat + return self.LM def update_mesh(self) -> None: """Set the resultant vertex loads in mesh data by using the latest calculated internal load matrix.""" try: - lmat = self.loads_mat + LM = self.LM except AttributeError: - lmat = self._vertex_loads + LM = self._vertex_loads for v, vkey in enumerate(self.mesh.vertices()): self.mesh.vertex_attributes(vkey, [self._PX, self._PY, self._PZ], - [lmat[v, 0], lmat[v, 1], lmat[v, 2]]) + [LM[v, 0], LM[v, 1], LM[v, 2]]) def _set_mesh_maps(self) -> None: """Set mesh counts and maps.""" self._v_count = self.mesh.number_of_vertices() self._f_count = self.mesh.number_of_faces() self._v_id = self.mesh.key_index() - self.face_normals = None - self.trib_area = None + self.FN = None + self.TA = None def _preprocess_loads(self) -> None: """Preprocess loads into arrays.""" @@ -127,18 +127,18 @@ def _preprocess_loads(self) -> None: ).reshape((self._v_count, 3)) rho = self.mesh.attributes[self.RHO] self._weight = asarray([-t * w * rho for t, w in - self.mesh.faces_attributes([self.THICKNESS, self.SW])] + self.mesh.faces_attributes([self.THICKNESS, self.HAS_SW])] ).reshape((self._f_count, 1)) - self._normal_load = asarray(self.mesh.faces_attribute(self.NORMAL) - ).reshape((self._f_count, 1)) - self._project_load = asarray(self.mesh.faces_attribute(self.PROJECT) + self._normal_loads = asarray(self.mesh.faces_attribute(self.NORMAL) ).reshape((self._f_count, 1)) + self._project_loads = asarray(self.mesh.faces_attribute(self.PROJECT) + ).reshape((self._f_count, 1)) def _set_load_flags(self) -> None: """Set flags for which load types are applied.""" bool_weight = array(self._weight, dtype=bool) - bool_normal = array(self._normal_load, dtype=bool) - bool_project = array(self._project_load, dtype=bool) + bool_normal = array(self._normal_loads, dtype=bool) + bool_project = array(self._project_loads, dtype=bool) self._oriented_loads = bool_normal + bool_project self._face_loads = self._oriented_loads + bool_weight self._load_flags = {'any_face_load': self._face_loads.any(), @@ -156,34 +156,34 @@ def _add_face_loads(self) -> None: self._add_projected_loads() def _process_faces(self, xyz: NDArray[(Any, 3), float64], - all_faces: bool) -> None: + process_all_faces: bool) -> None: """Calculate normals and tributary areas per face, for the vertex coordinates of the current geometry.""" - self.face_normals = empty((self._f_count, 3), dtype=float) + self.FN = empty((self._f_count, 3), dtype=float) trib_data = []; trib_rows = []; trib_cols = [] # noqa for f, fkey in enumerate(self.mesh.faces()): - if not (self._face_loads[f] or all_faces): + if not (self._face_loads[f] or process_all_faces): continue vts = [self._v_id[v] for v in self.mesh.face_vertices(fkey)] f_xyz = xyz[vts, :] - cps = self.__face_cross_products(f_xyz) + cps = self._face_cross_products(f_xyz) # face tributary areas per vertex part_areas = 0.25 * norm(cps, axis=1) trib_data.extend(part_areas + roll(part_areas, -1)) trib_rows.extend(vts) trib_cols.extend([f] * len(vts)) # face normal vector - if self._oriented_loads[f] or all_faces: - self.face_normals[f, :] = add.reduce(cps) + if self._oriented_loads[f] or process_all_faces: + self.FN[f, :] = add.reduce(cps) - self.face_normals /= norm(self.face_normals, axis=1)[:, None] - self.trib_areas = coo_matrix((trib_data, (trib_rows, trib_cols)), - (self._v_count, self._f_count)).tocsr() + self.FN /= norm(self.FN, axis=1)[:, None] + self.TA = coo_matrix((trib_data, (trib_rows, trib_cols)), + (self._v_count, self._f_count)).tocsr() - def __face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] - ) -> NDArray[(Any, 3), float64]: + def _face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] + ) -> NDArray[(Any, 3), float64]: """Utility to calculate cross products of the consecutive (centroid -> vertex) vectors for each of the vertices of the face.""" centroid = average(f_xyz, axis=0) @@ -194,17 +194,17 @@ def __face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] def _add_weight(self) -> None: """Convert all face self-weights into global vertex loads and add them into the global load matrix.""" - self.loads_mat[:, 2:] += self.trib_areas.dot(self._weight) + self.LM[:, 2:] += self.TA.dot(self._weight) def _add_normal_loads(self) -> None: """Convert all normal face loads into global vertex loads and add them into the global load matrix.""" - normal_loads = self.face_normals * self._normal_load - self.loads_mat += self.trib_areas.dot(normal_loads) + NL = self.FN * self._normal_loads + self.LM += self.TA.dot(NL) def _add_projected_loads(self) -> None: """Convert all Z-projected face loads into global vertex loads and add them into the global load matrix.""" z = asarray([0, 0, -1]).reshape((3, 1)) - proj_loads = self.face_normals.dot(z) * self._project_load - self.loads_mat[:, 2:] += self.trib_areas.dot(proj_loads) + PL = self.FN.dot(z) * self._project_loads + self.LM[:, 2:] += self.TA.dot(PL) From 3ee178b1064bc2cb7255d00254019adaadb71c12 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Tue, 26 Oct 2021 10:19:54 +0200 Subject: [PATCH 07/13] Naming again for user clarity --- src/compas_fd/loads/load_calculator.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 92971076..b1cde6de 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -35,8 +35,8 @@ class LoadCalculator: Attributes ---------- - LM : ndarray of float - global vertex load matrix as (v x 3) array. + RL : ndarray of float + global resultant vertex load matrix as (v x 3) array. FN : ndarray of float Face normals for current geometry as (f x 3) array. TA : ndarray of float @@ -96,21 +96,21 @@ def __call__(self, xyz: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, """ if not self._load_flags['any_face_load'] and not process_all_faces: return self._vertex_loads - self.LM = array(self._vertex_loads, copy=True) + self.RL = array(self._vertex_loads, copy=True) self._process_faces(asarray(xyz).reshape(-1, 3), process_all_faces) self._add_face_loads() - return self.LM + return self.RL def update_mesh(self) -> None: """Set the resultant vertex loads in mesh data by using the latest calculated internal load matrix.""" try: - LM = self.LM + RL = self.RL except AttributeError: - LM = self._vertex_loads + RL = self._vertex_loads for v, vkey in enumerate(self.mesh.vertices()): self.mesh.vertex_attributes(vkey, [self._PX, self._PY, self._PZ], - [LM[v, 0], LM[v, 1], LM[v, 2]]) + [RL[v, 0], RL[v, 1], RL[v, 2]]) def _set_mesh_maps(self) -> None: """Set mesh counts and maps.""" @@ -194,17 +194,17 @@ def _face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] def _add_weight(self) -> None: """Convert all face self-weights into global vertex loads and add them into the global load matrix.""" - self.LM[:, 2:] += self.TA.dot(self._weight) + self.RL[:, 2:] += self.TA.dot(self._weight) def _add_normal_loads(self) -> None: """Convert all normal face loads into global vertex loads and add them into the global load matrix.""" NL = self.FN * self._normal_loads - self.LM += self.TA.dot(NL) + self.RL += self.TA.dot(NL) def _add_projected_loads(self) -> None: """Convert all Z-projected face loads into global vertex loads and add them into the global load matrix.""" z = asarray([0, 0, -1]).reshape((3, 1)) PL = self.FN.dot(z) * self._project_loads - self.LM[:, 2:] += self.TA.dot(PL) + self.RL[:, 2:] += self.TA.dot(PL) From 52634451eadf9f0504a0b4ee8c8e117df36ab14a Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Tue, 26 Oct 2021 11:53:34 +0200 Subject: [PATCH 08/13] Add boolean toggles loads --- src/compas_fd/loads/load_calculator.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index b1cde6de..4b871973 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -50,7 +50,8 @@ class LoadCalculator: >>> from compas_fd import LoadCalculator >>> mesh = Mesh.from_obj(compas.get('hypar.obj')) >>> dva = {'px': .0, 'py': .0, 'pz': .1} - >>> dfa = {'t': .1, 'is_loaded': True, 'wind': 1.0} + >>> dfa = {'t': .1, 'wind': 1.0, 'snow': .0, + 'has_weight': True, 'has_wind': True, 'has_snow': False} >>> mesh.attributes.update({'density': 22.0}) >>> mesh.update_default_vertex_attributes(dva) >>> mesh.update_default_face_attributes(dfa) @@ -66,9 +67,11 @@ class LoadCalculator: _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads RHO = 'density' THICKNESS = 't' - HAS_SW = 'is_loaded' NORMAL = 'wind' PROJECT = 'snow' + HAS_SW = 'has_weight' + HAS_NORMAL = 'has_wind' + HAS_PROJ = 'has_snow' def __init__(self, mesh: Mesh) -> None: self.mesh = mesh @@ -126,17 +129,16 @@ def _preprocess_loads(self) -> None: names=(self.PX, self.PY, self.PZ)) ).reshape((self._v_count, 3)) rho = self.mesh.attributes[self.RHO] - self._weight = asarray([-t * w * rho for t, w in - self.mesh.faces_attributes([self.THICKNESS, self.HAS_SW])] - ).reshape((self._f_count, 1)) - self._normal_loads = asarray(self.mesh.faces_attribute(self.NORMAL) - ).reshape((self._f_count, 1)) - self._project_loads = asarray(self.mesh.faces_attribute(self.PROJECT) - ).reshape((self._f_count, 1)) + self._weights = asarray([-t * hw * rho for t, hw in self.mesh.faces_attributes( + [self.THICKNESS, self.HAS_SW])]).reshape((self._f_count, 1)) + self._normal_loads = asarray([n * hn for n, hn in self.mesh.faces_attributes( + [self.NORMAL, self.HAS_NORMAL])]).reshape((self._f_count, 1)) + self._project_loads = asarray([p * hp for p, hp in self.mesh.faces_attributes( + [self.PROJECT, self.HAS_PROJ])]).reshape((self._f_count, 1)) def _set_load_flags(self) -> None: """Set flags for which load types are applied.""" - bool_weight = array(self._weight, dtype=bool) + bool_weight = array(self._weights, dtype=bool) bool_normal = array(self._normal_loads, dtype=bool) bool_project = array(self._project_loads, dtype=bool) self._oriented_loads = bool_normal + bool_project @@ -194,7 +196,7 @@ def _face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] def _add_weight(self) -> None: """Convert all face self-weights into global vertex loads and add them into the global load matrix.""" - self.RL[:, 2:] += self.TA.dot(self._weight) + self.RL[:, 2:] += self.TA.dot(self._weights) def _add_normal_loads(self) -> None: """Convert all normal face loads into global vertex loads From 4394891d17a36cd7c2574ea327e2847cdcfb3879 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Thu, 28 Oct 2021 14:17:50 +0200 Subject: [PATCH 09/13] Add weight by vertex --- src/compas_fd/loads/load_calculator.py | 106 ++++++++++++++----------- 1 file changed, 60 insertions(+), 46 deletions(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 4b871973..c48cadda 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -3,9 +3,9 @@ from numpy import asarray from numpy import average from numpy import cross -from numpy import empty from numpy import float64 from numpy import roll +from numpy import zeros from numpy.linalg import norm from scipy.sparse import coo_matrix @@ -37,11 +37,11 @@ class LoadCalculator: ---------- RL : ndarray of float global resultant vertex load matrix as (v x 3) array. - FN : ndarray of float - Face normals for current geometry as (f x 3) array. TA : ndarray of float Tributary areas matrix as sparse (v x f) array. Entry a_ij holds tributary area of face j for vertex i. + FN : ndarray of float + Face normals for current geometry as (f x 3) array. Examples -------- @@ -49,9 +49,9 @@ class LoadCalculator: >>> from compas.datastructures import Mesh >>> from compas_fd import LoadCalculator >>> mesh = Mesh.from_obj(compas.get('hypar.obj')) - >>> dva = {'px': .0, 'py': .0, 'pz': .1} - >>> dfa = {'t': .1, 'wind': 1.0, 'snow': .0, - 'has_weight': True, 'has_wind': True, 'has_snow': False} + >>> dva = {'px': .0, 'py': .0, 'pz': -.1, 't': .1} + >>> dfa = {'wind': 1.0, 'snow': .0, 'has_weight': True, + 'has_wind': True, 'has_snow': False} >>> mesh.attributes.update({'density': 22.0}) >>> mesh.update_default_vertex_attributes(dva) >>> mesh.update_default_face_attributes(dfa) @@ -62,7 +62,7 @@ class LoadCalculator: >>> lc.update_mesh() """ - # attribute names in mesh data + # hardcoded attribute names in mesh data map PX, PY, PZ = 'px', 'py', 'pz' # input vertex loads _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads RHO = 'density' @@ -97,75 +97,86 @@ def __call__(self, xyz: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, ndarray of float Resultant global vertex loads as (v x 3) array. """ - if not self._load_flags['any_face_load'] and not process_all_faces: + if not self._any_face_load and not process_all_faces: return self._vertex_loads - self.RL = array(self._vertex_loads, copy=True) + self._RL = array(self._vertex_loads, copy=True) self._process_faces(asarray(xyz).reshape(-1, 3), process_all_faces) self._add_face_loads() - return self.RL + return self._RL def update_mesh(self) -> None: - """Set the resultant vertex loads in mesh data by + """Set the resultant vertex loads in the mesh data map by using the latest calculated internal load matrix.""" - try: - RL = self.RL - except AttributeError: - RL = self._vertex_loads + RL = self.RL for v, vkey in enumerate(self.mesh.vertices()): self.mesh.vertex_attributes(vkey, [self._PX, self._PY, self._PZ], [RL[v, 0], RL[v, 1], RL[v, 2]]) + @property + def RL(self): + """global resultant vertex load matrix as (v x 3) array.""" + return getattr(self, '_RL', self._vertex_loads) + + @property + def TA(self): + """Tributary areas matrix as sparse (v x f) array. + Entry a_ij holds tributary area of face j for vertex i.""" + return getattr(self, '_TA', zeros((self._v_count, self._f_count))) + + @property + def FN(self): + """Face normals for current geometry as (f x 3) array.""" + return getattr(self, '_FN', zeros((self._f_count, 3))) + def _set_mesh_maps(self) -> None: """Set mesh counts and maps.""" self._v_count = self.mesh.number_of_vertices() self._f_count = self.mesh.number_of_faces() self._v_id = self.mesh.key_index() - self.FN = None - self.TA = None def _preprocess_loads(self) -> None: """Preprocess loads into arrays.""" self._vertex_loads = asarray(self.mesh.vertices_attributes( names=(self.PX, self.PY, self.PZ)) ).reshape((self._v_count, 3)) - rho = self.mesh.attributes[self.RHO] - self._weights = asarray([-t * hw * rho for t, hw in self.mesh.faces_attributes( - [self.THICKNESS, self.HAS_SW])]).reshape((self._f_count, 1)) + self._v_weights = self.mesh.attributes[self.RHO] * asarray( + self.mesh.vertices_attribute(self.THICKNESS)) self._normal_loads = asarray([n * hn for n, hn in self.mesh.faces_attributes( [self.NORMAL, self.HAS_NORMAL])]).reshape((self._f_count, 1)) self._project_loads = asarray([p * hp for p, hp in self.mesh.faces_attributes( [self.PROJECT, self.HAS_PROJ])]).reshape((self._f_count, 1)) def _set_load_flags(self) -> None: - """Set flags for which load types are applied.""" - bool_weight = array(self._weights, dtype=bool) - bool_normal = array(self._normal_loads, dtype=bool) - bool_project = array(self._project_loads, dtype=bool) - self._oriented_loads = bool_normal + bool_project - self._face_loads = self._oriented_loads + bool_weight - self._load_flags = {'any_face_load': self._face_loads.any(), - 'any_weight': bool_weight.any(), - 'any_normal': bool_normal.any(), - 'any_projected': bool_project.any()} + """Set flags for which face load types are applied.""" + have_normal = array(self._normal_loads, dtype=bool) + have_projected = array(self._project_loads, dtype=bool) + self._have_weight = array(self.mesh.faces_attribute(self.HAS_SW), + dtype=bool).reshape((self._f_count, 1)) + self._have_oriented_load = have_normal + have_projected + self._have_face_load = self._have_oriented_load + self._have_weight + self._any_face_load = self._have_face_load.any() + self._any_weight = self._have_weight.any() + self._any_normal = have_normal.any() + self._any_projected = have_projected.any() def _add_face_loads(self) -> None: """Add active face loads to the global load matrix.""" - if self._load_flags['any_weight']: + if self._any_weight: self._add_weight() - if self._load_flags['any_normal']: + if self._any_normal: self._add_normal_loads() - if self._load_flags['any_projected']: + if self._any_projected: self._add_projected_loads() def _process_faces(self, xyz: NDArray[(Any, 3), float64], process_all_faces: bool) -> None: """Calculate normals and tributary areas per face, for the vertex coordinates of the current geometry.""" - self.FN = empty((self._f_count, 3), dtype=float) + self._FN = zeros((self._f_count, 3), dtype=float) trib_data = []; trib_rows = []; trib_cols = [] # noqa for f, fkey in enumerate(self.mesh.faces()): - if not (self._face_loads[f] or process_all_faces): + if not (self._have_face_load[f] or process_all_faces): continue vts = [self._v_id[v] for v in self.mesh.face_vertices(fkey)] @@ -177,12 +188,12 @@ def _process_faces(self, xyz: NDArray[(Any, 3), float64], trib_rows.extend(vts) trib_cols.extend([f] * len(vts)) # face normal vector - if self._oriented_loads[f] or process_all_faces: - self.FN[f, :] = add.reduce(cps) + if self._have_oriented_load[f] or process_all_faces: + self._FN[f, :] = add.reduce(cps) - self.FN /= norm(self.FN, axis=1)[:, None] - self.TA = coo_matrix((trib_data, (trib_rows, trib_cols)), - (self._v_count, self._f_count)).tocsr() + self._FN /= norm(self._FN, axis=1)[:, None] + self._TA = coo_matrix((trib_data, (trib_rows, trib_cols)), + (self._v_count, self._f_count)).tocsr() def _face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] ) -> NDArray[(Any, 3), float64]: @@ -195,18 +206,21 @@ def _face_cross_products(self, f_xyz: NDArray[(Any, 3), float64] def _add_weight(self) -> None: """Convert all face self-weights into global vertex loads - and add them into the global load matrix.""" - self.RL[:, 2:] += self.TA.dot(self._weights) + and add them into the global load matrix. + Loads are stepped per vertex tributary area.""" + self._RL[:, 2:] += (self._v_weights * add.reduce( + self._TA * self._have_weight, axis=1) + ).reshape((-1, 1)) def _add_normal_loads(self) -> None: """Convert all normal face loads into global vertex loads and add them into the global load matrix.""" - NL = self.FN * self._normal_loads - self.RL += self.TA.dot(NL) + NL = self._FN * self._normal_loads + self._RL += self._TA.dot(NL) def _add_projected_loads(self) -> None: """Convert all Z-projected face loads into global vertex loads and add them into the global load matrix.""" z = asarray([0, 0, -1]).reshape((3, 1)) - PL = self.FN.dot(z) * self._project_loads - self.RL[:, 2:] += self.TA.dot(PL) + PL = self._FN.dot(z) * self._project_loads + self._RL[:, 2:] += self._TA.dot(PL) From 49a4ac744b1562b41601276c2a0fab166be8ac0b Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 1 Nov 2021 09:34:29 +0100 Subject: [PATCH 10/13] Weight minus Z --- src/compas_fd/loads/load_calculator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index c48cadda..1a6f10c2 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -208,7 +208,7 @@ def _add_weight(self) -> None: """Convert all face self-weights into global vertex loads and add them into the global load matrix. Loads are stepped per vertex tributary area.""" - self._RL[:, 2:] += (self._v_weights * add.reduce( + self._RL[:, 2:] -= (self._v_weights * add.reduce( self._TA * self._have_weight, axis=1) ).reshape((-1, 1)) From e24d3f92564ac242d24e96256eb4f107f0dafaa3 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 1 Nov 2021 10:30:36 +0100 Subject: [PATCH 11/13] Include global load factors --- src/compas_fd/datastructures/cablemesh.py | 8 ++++++- src/compas_fd/loads/load_calculator.py | 26 +++++++++++++++++------ 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/src/compas_fd/datastructures/cablemesh.py b/src/compas_fd/datastructures/cablemesh.py index 05fb1c9f..4df861c5 100644 --- a/src/compas_fd/datastructures/cablemesh.py +++ b/src/compas_fd/datastructures/cablemesh.py @@ -23,6 +23,10 @@ def __init__(self): super(CableMesh, self).__init__() self.attributes.update({ 'density': 22.0, + 'point_lc': 1.0, + 'weight_lc': 1.0, + 'wind_lc': 1.0, + 'snow_lc': 1.0 }) self.default_vertex_attributes.update({ 'x': 0.0, @@ -51,7 +55,9 @@ def __init__(self): '_r': 0.0, }) self.default_face_attributes.update({ - 'is_loaded': False, + 'has_weight': False, + 'has_snow': False, + 'has_wind': False, 't': 0.0, 'wind': 0.0, 'snow': 0.0 diff --git a/src/compas_fd/loads/load_calculator.py b/src/compas_fd/loads/load_calculator.py index 1a6f10c2..8c1c9f2f 100644 --- a/src/compas_fd/loads/load_calculator.py +++ b/src/compas_fd/loads/load_calculator.py @@ -49,10 +49,12 @@ class LoadCalculator: >>> from compas.datastructures import Mesh >>> from compas_fd import LoadCalculator >>> mesh = Mesh.from_obj(compas.get('hypar.obj')) + >>> ma = {'density': 22.0 'point_lc': 1, 'weight_lc': 1, + 'wind_lc': 1, 'snow_lc': 1} >>> dva = {'px': .0, 'py': .0, 'pz': -.1, 't': .1} >>> dfa = {'wind': 1.0, 'snow': .0, 'has_weight': True, 'has_wind': True, 'has_snow': False} - >>> mesh.attributes.update({'density': 22.0}) + >>> mesh.attributes.update(ma) >>> mesh.update_default_vertex_attributes(dva) >>> mesh.update_default_face_attributes(dfa) >>> lc = LoadCalculator(mesh) @@ -65,8 +67,14 @@ class LoadCalculator: # hardcoded attribute names in mesh data map PX, PY, PZ = 'px', 'py', 'pz' # input vertex loads _PX, _PY, _PZ = '_px', '_py', '_pz' # resultant vertex loads + RHO = 'density' - THICKNESS = 't' + VTS_LC = 'point_lc' + WEIGHT_LC = 'weight_lc' + NORMAL_LC = 'wind_lc' + PROJ_LC = 'snow_lc' + + TH = 't' NORMAL = 'wind' PROJECT = 'snow' HAS_SW = 'has_weight' @@ -136,15 +144,19 @@ def _set_mesh_maps(self) -> None: def _preprocess_loads(self) -> None: """Preprocess loads into arrays.""" + vts_lc, weight_lc, normal_lc, proj_lc, rho = (self.mesh.attributes[key] for key in ( + self.VTS_LC, self.WEIGHT_LC, self.NORMAL_LC, self.PROJ_LC, self.RHO)) + self._vertex_loads = asarray(self.mesh.vertices_attributes( names=(self.PX, self.PY, self.PZ)) - ).reshape((self._v_count, 3)) - self._v_weights = self.mesh.attributes[self.RHO] * asarray( - self.mesh.vertices_attribute(self.THICKNESS)) + ).reshape((self._v_count, 3)) * vts_lc + self._v_weights = asarray(self.mesh.vertices_attribute(self.TH)) * rho * weight_lc self._normal_loads = asarray([n * hn for n, hn in self.mesh.faces_attributes( - [self.NORMAL, self.HAS_NORMAL])]).reshape((self._f_count, 1)) + [self.NORMAL, self.HAS_NORMAL])] + ).reshape((self._f_count, 1)) * normal_lc self._project_loads = asarray([p * hp for p, hp in self.mesh.faces_attributes( - [self.PROJECT, self.HAS_PROJ])]).reshape((self._f_count, 1)) + [self.PROJECT, self.HAS_PROJ])] + ).reshape((self._f_count, 1)) * proj_lc def _set_load_flags(self) -> None: """Set flags for which face load types are applied.""" From 59486b7d1f1ed33821cb66b0c1a06b9b6447b5f2 Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 1 Nov 2021 10:52:27 +0100 Subject: [PATCH 12/13] Move thickness to vertices --- src/compas_fd/datastructures/cablemesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compas_fd/datastructures/cablemesh.py b/src/compas_fd/datastructures/cablemesh.py index 4df861c5..bdf01026 100644 --- a/src/compas_fd/datastructures/cablemesh.py +++ b/src/compas_fd/datastructures/cablemesh.py @@ -46,6 +46,7 @@ def __init__(self): '_rx': 0.0, '_ry': 0.0, '_rz': 0.0, + 't': 0.0 }) self.default_edge_attributes.update({ 'q': 1.0, @@ -58,7 +59,6 @@ def __init__(self): 'has_weight': False, 'has_snow': False, 'has_wind': False, - 't': 0.0, 'wind': 0.0, 'snow': 0.0 }) From a9fd86fe0dd3b07ed7161107a562ad4ab2d0a5df Mon Sep 17 00:00:00 2001 From: Sam <Sam_Bouten@hotmail.com> Date: Mon, 1 Nov 2021 10:57:40 +0100 Subject: [PATCH 13/13] Correct thickness --- src/compas_fd/datastructures/cablemesh.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/compas_fd/datastructures/cablemesh.py b/src/compas_fd/datastructures/cablemesh.py index bdf01026..ec470f62 100644 --- a/src/compas_fd/datastructures/cablemesh.py +++ b/src/compas_fd/datastructures/cablemesh.py @@ -45,8 +45,7 @@ def __init__(self): '_pz': 0.0, '_rx': 0.0, '_ry': 0.0, - '_rz': 0.0, - 't': 0.0 + '_rz': 0.0 }) self.default_edge_attributes.update({ 'q': 1.0,