diff --git a/CHANGELOG.md b/CHANGELOG.md index b53bdef6..7b64eece 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +* Added `Envelope` base class for masonry structure boundaries with intrados, extrados, and middle meshes +* Added especialized `MeshEnvelope`, `ParametricEnvelope`, and `BrepEnvelope` to deal with different evelope inputs. +* Added direct parametric envelope creation methods inheriting from `ParametricEnvelope` class: + * `CrossVaultEnvelope` - Creates cross vault envelopes with configurable spans and thickness + * `DomeEnvelope` - Creates dome envelopes with configurable radius, center, and oculus + * `PavillionVaultEnvelope` - Creates pavillion vault envelopes with spring angle support + * `PointedVaultEnvelope` - Creates pointed vault envelopes with height control parameters +* The infrastructure around these classes has been updated to enable assigning constraints to the form diagrams. * Added comprehensive form diagram creation methods to `FormDiagram` class: * `create_cross()` - Creates cross discretisation with orthogonal arrangement and quad diagonals * `create_fan()` - Creates fan discretisation with straight lines to corners @@ -25,30 +33,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Added corner detection functionality: * `corner_vertices()` - Identifies corner vertices on boundary with configurable angle threshold * Added comprehensive mesh creation functions in `diagram_rectangular.py`: - * `create_cross_mesh()` - Creates cross pattern meshes - * `create_fan_mesh()` - Creates fan pattern meshes - * `create_parametric_fan_mesh()` - Creates parametric fan meshes with lambda interpolation - * `create_cross_with_diagonal_mesh()` - Creates cross meshes with diagonals - * `create_ortho_mesh()` - Creates orthogonal grid meshes * Added circular mesh creation functions in `diagram_circular.py`: - * `create_circular_radial_mesh()` - Creates circular radial meshes - * `create_circular_radial_spaced_mesh()` - Creates hemispherically spaced circular meshes - * `create_circular_spiral_mesh()` - Creates circular spiral meshes * Added arch mesh creation functions in `diagram_arch.py`: - * `create_arch_linear_mesh()` - Creates arch meshes with semicircular projection - * `create_arch_linear_equally_spaced_mesh()` - Creates arch meshes with equally spaced nodes ### Changed -* Refactored parameter passing from `xy_span=[[x0, x1], [y0, y1]]` to `x_span=(x0, x1), y_span=(y0, y1)` for better API consistency -* Changed `center` parameter from list to tuple in circular diagram functions for immutability and consistency -* Updated function signatures to use single-line format for Black compatibility -* Improved code formatting and linting across all diagram generation files -* Fixed various spelling errors and documentation formatting issues -* Renamed `create_delta_form()` to `create_cross_with_diagonal()` for better clarity -* Updated support assignment to use `is_support` attribute instead of deprecated `is_fixed` -* Removed deprecated `form.parameters` usage from all form creation methods - ### Removed diff --git a/src/compas_tna/envelope/__init__.py b/src/compas_tna/envelope/__init__.py new file mode 100644 index 00000000..bab07f13 --- /dev/null +++ b/src/compas_tna/envelope/__init__.py @@ -0,0 +1,20 @@ +from .envelope import Envelope +from .brepenvelope import BrepEnvelope +from .meshenvelope import MeshEnvelope +from .parametricenvelope import ParametricEnvelope + +from .crossvault import CrossVaultEnvelope +from .dome import DomeEnvelope +from .pavillionvault import PavillionVaultEnvelope +from .pointedvault import PointedVaultEnvelope + +__all__ = [ + "Envelope", + "BrepEnvelope", + "MeshEnvelope", + "ParametricEnvelope", + "PavillionVaultEnvelope", + "PointedVaultEnvelope", + "DomeEnvelope", + "CrossVaultEnvelope", +] diff --git a/src/compas_tna/envelope/brepenvelope.py b/src/compas_tna/envelope/brepenvelope.py new file mode 100644 index 00000000..a94268a1 --- /dev/null +++ b/src/compas_tna/envelope/brepenvelope.py @@ -0,0 +1,12 @@ +from compas.geometry import Brep +from compas_tna.envelope import Envelope + + +class BrepEnvelope(Envelope): + """An Envelope defined by a BRep at intrados, extrados, and middle.""" + + def __init__(self, intrados: Brep = None, extrados: Brep = None, middle: Brep = None, **kwargs): + super().__init__(**kwargs) + self.intrados = intrados + self.extrados = extrados + self.middle = middle diff --git a/src/compas_tna/envelope/crossvault.py b/src/compas_tna/envelope/crossvault.py new file mode 100644 index 00000000..a35f6838 --- /dev/null +++ b/src/compas_tna/envelope/crossvault.py @@ -0,0 +1,458 @@ +import math + +from numpy import array +from numpy import ones +from numpy import zeros + +from compas.datastructures import Mesh +from compas_tna.diagrams.diagram_rectangular import create_cross_mesh +from compas_tna.envelope.parametricenvelope import ParametricEnvelope + + +def crossvault_envelope( + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, +): + """Create an envelope for a cross vault geometry with given parameters. + + Parameters + ---------- + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + thickness : float, optional + Thickness of the vault, by default 0.50 + min_lb : float, optional + Parameter for lower bound in nodes in the boundary, by default 0.0 + n : int, optional + Number of vertices for the mesh, by default 100 + + Returns + ------- + middle : Mesh + Middle mesh + intrados : Mesh + Intrados mesh + extrados : Mesh + Extrados mesh + """ + # Create base topology + base_topology = create_cross_mesh(x_span=x_span, y_span=y_span, n=n) + xyz0, faces_i = base_topology.to_vertices_and_faces() + xi, yi, _ = array(xyz0).transpose() + + # Create middle surface + zt = crossvault_middle(xi, yi, min_lb, x_span=x_span, y_span=y_span, tol=1e-6) + xyzt = array([xi, yi, zt.flatten()]).transpose() + middle = Mesh.from_vertices_and_faces(xyzt, faces_i) + middle.update_default_vertex_attributes(thickness=thickness) + + # Create upper and lower bounds + zub, zlb = crossvault_bounds(xi, yi, thickness, min_lb, x_span=x_span, y_span=y_span, tol=1e-6) + xyzub = array([xi, yi, zub.flatten()]).transpose() + xyzlb = array([xi, yi, zlb.flatten()]).transpose() + + extrados = Mesh.from_vertices_and_faces(xyzub, faces_i) + intrados = Mesh.from_vertices_and_faces(xyzlb, faces_i) + + return intrados, extrados, middle + + +def crossvault_middle(x, y, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Compute middle of a crossvault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + span of the vault in y direction, by default (0.0, 10.0) + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + z : array + Values of the middle surface in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + + rx = (x1 - x0) / 2 + ry = (y1 - y0) / 2 + hc = max(rx, ry) + + z = zeros((len(x), 1)) + + for i in range(len(x)): + xi, yi = x[i], y[i] + if yi > y1: + yi = y1 + if yi < y0: + yi = y0 + if xi > x1: + xi = x1 + if xi < x0: + xi = x0 + xd = x0 + (x1 - x0) / (y1 - y0) * (yi - y0) + yd = y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + hxd = math.sqrt(abs((rx) ** 2 - ((xd - x0) - rx) ** 2)) + hyd = math.sqrt(abs((ry) ** 2 - ((yd - y0) - ry) ** 2)) + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + z[i] = hc * (hxd + math.sqrt((ry) ** 2 - ((yi - y0) - ry) ** 2)) / (rx + ry) + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + z[i] = hc * (hyd + math.sqrt((rx) ** 2 - ((xi - x0) - rx) ** 2)) / (rx + ry) + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + z[i] = hc * (hxd + math.sqrt((ry) ** 2 - ((yi - y0) - ry) ** 2)) / (rx + ry) + elif yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + z[i] = hc * (hyd + math.sqrt((rx) ** 2 - ((xi - x0) - rx) ** 2)) / (rx + ry) + else: + print("Vertex did not belong to any Q. (x,y) = ({0},{1})".format(xi, yi)) + z[i] = -min_lb + + return z + + +def crossvault_bounds(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Compute upper and lower bounds of an crossvault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the arch + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + span of the vault in y direction, by default (0.0, 10.0) + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + ub : array + Values of the upper bound in the points + lb : array + Values of the lower bound in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + + y1_ub = y1 + thk / 2 + y0_ub = y0 - thk / 2 + x1_ub = x1 + thk / 2 + x0_ub = x0 - thk / 2 + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + rx_ub = (x1_ub - x0_ub) / 2 + ry_ub = (y1_ub - y0_ub) / 2 + rx_lb = (x1_lb - x0_lb) / 2 + ry_lb = (y1_lb - y0_lb) / 2 + + hc_ub = max(rx_ub, ry_ub) + hc_lb = max(rx_lb, ry_lb) + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + + for i in range(len(x)): + xi, yi = x[i], y[i] + xd_ub = x0_ub + (x1_ub - x0_ub) / (y1_ub - y0_ub) * (yi - y0_ub) + yd_ub = y0_ub + (y1_ub - y0_ub) / (x1_ub - x0_ub) * (xi - x0_ub) + hxd_ub = math.sqrt((rx_ub) ** 2 - ((xd_ub - x0_ub) - rx_ub) ** 2) + hyd_ub = math.sqrt((ry_ub) ** 2 - ((yd_ub - y0_ub) - ry_ub) ** 2) + + intrados_null = False + + if (yi > y1_lb and (xi > x1_lb or xi < x0_lb)) or (yi < y0_lb and (xi > x1_lb or xi < x0_lb)): + intrados_null = True + else: + yi_intra = yi + xi_intra = xi + if yi > y1_lb: + yi_intra = y1_lb + if yi < y0_lb: + yi_intra = y0_lb + elif xi > x1_lb: + xi_intra = x1_lb + elif xi < x0_lb: + xi_intra = x0_lb + + xd_lb = x0_lb + (x1_lb - x0_lb) / (y1_lb - y0_lb) * (yi_intra - y0_lb) + yd_lb = y0_lb + (y1_lb - y0_lb) / (x1_lb - x0_lb) * (xi_intra - x0_lb) + hxd_lb = _sqrt(((rx_lb) ** 2 - ((xd_lb - x0_lb) - rx_lb) ** 2)) + hyd_lb = _sqrt(((ry_lb) ** 2 - ((yd_lb - y0_lb) - ry_lb) ** 2)) + + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + ub[i] = hc_ub * (hxd_ub + math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2)) / (rx_ub + ry_ub) + if not intrados_null: + lb[i] = hc_lb * (hxd_lb + math.sqrt((ry_lb) ** 2 - ((yi_intra - y0_lb) - ry_lb) ** 2)) / (rx_lb + ry_lb) + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + ub[i] = hc_ub * (hyd_ub + math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2)) / (rx_ub + ry_ub) + if not intrados_null: + lb[i] = hc_lb * (hyd_lb + math.sqrt((rx_lb) ** 2 - ((xi_intra - x0_lb) - rx_lb) ** 2)) / (rx_lb + ry_lb) + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + ub[i] = hc_ub * (hxd_ub + math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2)) / (rx_ub + ry_ub) + if not intrados_null: + lb[i] = hc_lb * (hxd_lb + math.sqrt((ry_lb) ** 2 - ((yi_intra - y0_lb) - ry_lb) ** 2)) / (rx_lb + ry_lb) + elif yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + ub[i] = hc_ub * (hyd_ub + math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2)) / (rx_ub + ry_ub) + if not intrados_null: + lb[i] = hc_lb * (hyd_lb + math.sqrt((rx_lb) ** 2 - ((xi_intra - x0_lb) - rx_lb) ** 2)) / (rx_lb + ry_lb) + else: + print("Error Q. (x,y) = ({0},{1})".format(xi, yi)) + + return ub, lb + + +def crossvault_bounds_derivatives(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Computes the sensitivities of upper and lower bounds in the x, y coordinates and thickness specified. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the arch + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + dub : array + Values of the sensitivities for the upper bound in the points + dlb : array + Values of the sensitivities for the lower bound in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + + y1_ub = y1 + thk / 2 + y0_ub = y0 - thk / 2 + x1_ub = x1 + thk / 2 + x0_ub = x0 - thk / 2 + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + rx_ub = (x1_ub - x0_ub) / 2 + ry_ub = (y1_ub - y0_ub) / 2 + rx_lb = (x1_lb - x0_lb) / 2 + ry_lb = (y1_lb - y0_lb) / 2 + + hc_ub = max(rx_ub, ry_ub) + hc_lb = max(rx_lb, ry_lb) + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + dub = zeros((len(x), 1)) # dzub / dt + dlb = zeros((len(x), 1)) # dzlb / dt + + dubdx = zeros((len(x), len(x))) + dubdy = zeros((len(x), len(x))) + dlbdx = zeros((len(x), len(x))) + dlbdy = zeros((len(x), len(x))) + + yc = ry_ub + y0_ub # Only works for square + xc = rx_ub + x0_ub + + for i in range(len(x)): + xi, yi = x[i], y[i] + xd_ub = x0_ub + (x1_ub - x0_ub) / (y1_ub - y0_ub) * (yi - y0_ub) + yd_ub = y0_ub + (y1_ub - y0_ub) / (x1_ub - x0_ub) * (xi - x0_ub) + hxd_ub = math.sqrt((rx_ub) ** 2 - ((xd_ub - x0_ub) - rx_ub) ** 2) + hyd_ub = math.sqrt((ry_ub) ** 2 - ((yd_ub - y0_ub) - ry_ub) ** 2) + + intrados_null = False + + if (yi > y1_lb and (xi > x1_lb or xi < x0_lb)) or (yi < y0_lb and (xi > x1_lb or xi < x0_lb)): + intrados_null = True + else: + yi_intra = yi + xi_intra = xi + if yi > y1_lb: + yi_intra = y1_lb + elif yi < y0_lb: + yi_intra = y0_lb + if xi > x1_lb: + xi_intra = x1_lb + elif xi < x0_lb: + xi_intra = x0_lb + + xd_lb = x0_lb + (x1_lb - x0_lb) / (y1_lb - y0_lb) * (yi_intra - y0_lb) + yd_lb = y0_lb + (y1_lb - y0_lb) / (x1_lb - x0_lb) * (xi_intra - x0_lb) + hxd_lb = _sqrt(((rx_lb) ** 2 - ((xd_lb - x0_lb) - rx_lb) ** 2)) + hyd_lb = _sqrt(((ry_lb) ** 2 - ((yd_lb - y0_lb) - ry_lb) ** 2)) + + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + ub[i] = hc_ub * (hxd_ub + math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2)) / (rx_ub + ry_ub) + dub[i] = 1 / 2 * ry_ub / ub[i] * hc_ub / ((rx_ub + ry_ub) / 2) + # dubdx[i, i] += 0.0 + dubdy[i, i] += -(yi - yc) / ub[i] + if not intrados_null: + lb[i] = hc_lb * (hxd_lb + math.sqrt((ry_lb) ** 2 - ((yi_intra - y0_lb) - ry_lb) ** 2)) / (rx_lb + ry_lb) + dlb[i] = -1 / 2 * ry_lb / lb[i] * hc_lb / ((rx_lb + ry_lb) / 2) + # dlbdx[i, i] += 0.0 + dlbdy[i, i] += -(yi - yc) / lb[i] + if yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + ub[i] = hc_ub * (hyd_ub + math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2)) / (rx_ub + ry_ub) + dub[i] = 1 / 2 * rx_ub / ub[i] * hc_ub / ((rx_ub + ry_ub) / 2) + # dubdy[i, i] += 0.0 + dubdx[i, i] += -(xi - xc) / ub[i] + if not intrados_null: + lb[i] = hc_lb * (hyd_lb + math.sqrt((rx_lb) ** 2 - ((xi_intra - x0_lb) - rx_lb) ** 2)) / (rx_lb + ry_lb) + dlb[i] = -1 / 2 * rx_lb / lb[i] * hc_lb / ((rx_lb + ry_lb) / 2) + # dlbdy[i, i] += 0.0 + dlbdx[i, i] += -(xi - xc) / lb[i] + if yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + ub[i] = hc_ub * (hxd_ub + math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2)) / (rx_ub + ry_ub) + dub[i] = 1 / 2 * ry_ub / ub[i] * hc_ub / ((rx_ub + ry_ub) / 2) + # dubdx[i, i] += 0.0 + dubdy[i, i] += -(yi - yc) / ub[i] + if not intrados_null: + lb[i] = hc_lb * (hxd_lb + math.sqrt((ry_lb) ** 2 - ((yi_intra - y0_lb) - ry_lb) ** 2)) / (rx_lb + ry_lb) + dlb[i] = -1 / 2 * ry_lb / lb[i] * hc_lb / ((rx_lb + ry_lb) / 2) + # dlbdx[i, i] += 0.0 + dlbdy[i, i] += -(yi - yc) / lb[i] + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + ub[i] = hc_ub * (hyd_ub + math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2)) / (rx_ub + ry_ub) + dub[i] = 1 / 2 * rx_ub / ub[i] * hc_ub / ((rx_ub + ry_ub) / 2) + # dubdy[i, i] += 0.0 + dubdx[i, i] += -(xi - xc) / ub[i] + if not intrados_null: + lb[i] = hc_lb * (hyd_lb + math.sqrt((rx_lb) ** 2 - ((xi_intra - x0_lb) - rx_lb) ** 2)) / (rx_lb + ry_lb) + dlb[i] = -1 / 2 * rx_lb / lb[i] * hc_lb / ((rx_lb + ry_lb) / 2) + # dlbdy[i, i] += 0.0 + dlbdx[i, i] += -(xi - xc) / lb[i] + # else: + # print('Error Q. (x,y) = ({0},{1})'.format(xi, yi)) + + return dub, dlb, dubdx, dubdy, dlbdx, dlbdy + + +def crossvault_bound_react(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Compute the bounds on the reaction vector of the crossvault.""" + pass + + +def crossvault_bound_react_derivatives(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Compute the sensitivities of the bounds on the reaction vector of the crossvault.""" + pass + + +def _sqrt(x): + try: + sqrt_x = math.sqrt(x) + except BaseException: + if x > -10e4: + sqrt_x = math.sqrt(abs(x)) + else: + sqrt_x = 0.0 + print("Problems to sqrt: ", x) + return sqrt_x + + +class CrossVaultEnvelope(ParametricEnvelope): + def __init__( + self, + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, + **kwargs, + ): + super().__init__(thickness=thickness, **kwargs) + self.x_span = x_span + self.y_span = y_span + self.min_lb = min_lb + self.n = n + + self.update_envelope() # Generate the intra/extra/middle meshes + + @property + def __data__(self): + data = super().__data__ + data["x_span"] = self.x_span + data["y_span"] = self.y_span + data["min_lb"] = self.min_lb + data["n"] = self.n + return data + + def __str__(self): + return f"CrossVaultEnvelope(name={self.name})" + + def update_envelope(self): + intrados, extrados, middle = crossvault_envelope(x_span=self.x_span, y_span=self.y_span, thickness=self.thickness, min_lb=self.min_lb, n=self.n) + self.intrados = intrados + self.extrados = extrados + self.middle = middle + + def compute_middle(self, x, y): + return crossvault_middle(x, y, self.min_lb, self.x_span, self.y_span, tol=1e-6) + + def compute_bounds(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return crossvault_bounds(x, y, thickness, self.min_lb, self.x_span, self.y_span, tol=1e-6) + + def compute_bounds_derivatives(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return crossvault_bounds_derivatives(x, y, thickness, self.min_lb, self.x_span, self.y_span, tol=1e-6) + + def compute_bound_react(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return crossvault_bound_react(x, y, thickness, self.min_lb, self.x_span, self.y_span, tol=1e-6) + + def compute_bound_react_derivatives(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return crossvault_bound_react_derivatives(x, y, thickness, self.min_lb, self.x_span, self.y_span, tol=1e-6) diff --git a/src/compas_tna/envelope/dome.py b/src/compas_tna/envelope/dome.py new file mode 100644 index 00000000..19ce1b0e --- /dev/null +++ b/src/compas_tna/envelope/dome.py @@ -0,0 +1,364 @@ +import math + +from numpy import array +from numpy import ones +from numpy import zeros + +from compas.datastructures import Mesh +from compas_tna.diagrams.diagram_circular import create_circular_radial_spaced_mesh +from compas_tna.envelope.parametricenvelope import ParametricEnvelope + + +def dome_envelope( + center: tuple = (5.0, 5.0), + radius: float = 5.0, + thickness: float = 0.50, + min_lb: float = 0.0, + n_hoops: int = 24, + n_parallels: int = 40, + r_oculus: float = 0.0, +): + """Create an envelope for a dome geometry with given parameters. + + Parameters + ---------- + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + radius : float, optional + The radius of the dome, by default 5.0 + thickness : float, optional + Thickness of the dome, by default 0.50 + min_lb : float, optional + Parameter for lower bound in nodes in the boundary, by default 0.0 + n_hoops : int, optional + Number of hoops for the mesh, by default 24 + n_parallels : int, optional + Number of parallels for the mesh, by default 40 + r_oculus : float, optional + Radius of the oculus (opening at the top), by default 0.0 + + Returns + ------- + middle : Mesh + Middle mesh + intrados : Mesh + Intrados mesh + extrados : Mesh + Extrados mesh + """ + # Create meshes for different radii + for radius_current in [radius, radius - thickness / 2, radius + thickness / 2]: + base_topology = create_circular_radial_spaced_mesh( + center=center, + radius=radius_current, + n_hoops=n_hoops, + n_parallels=n_parallels, + r_oculus=r_oculus, + ) + xyz0, faces_i = base_topology.to_vertices_and_faces() + xi, yi, _ = array(xyz0).transpose() + zt = dome_middle(xi, yi, radius_current, min_lb, center=center) + xyzt = array([xi, yi, zt.flatten()]).transpose() + + if radius_current == radius: + middle = Mesh.from_vertices_and_faces(xyzt, faces_i) + elif radius_current == radius - thickness / 2: + intrados = Mesh.from_vertices_and_faces(xyzt, faces_i) + elif radius_current == radius + thickness / 2: + extrados = Mesh.from_vertices_and_faces(xyzt, faces_i) + + # Set thickness attributes + middle.update_default_vertex_attributes(thickness=thickness) + intrados.update_default_vertex_attributes(thickness=thickness) + extrados.update_default_vertex_attributes(thickness=thickness) + + return intrados, extrados, middle + + +def dome_middle(x, y, radius, min_lb, center=(5.0, 5.0)): + """Compute middle of the dome based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + radius : float, optional + The radius of the dome, by default 5.0 + min_lb : float + Parameter for lower bound in nodes in the boundary + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + + Returns + ------- + zt : array + Values of the middle surface in the points + """ + + xc = center[0] + yc = center[1] + zt = ones((len(x), 1)) + + for i in range(len(x)): + zt2 = radius**2 - (x[i] - xc) ** 2 - (y[i] - yc) ** 2 + if zt2 > 0: + zt[i] = math.sqrt(zt2) + else: + zt[i] = 0.0 + + return zt + + +def dome_bounds(x, y, thk, min_lb, center=(5.0, 5.0), radius=5.0): + """Compute upper and lower bounds of the dome based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the dome + min_lb : float + Parameter for lower bound in nodes in the boundary + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + radius : float, optional + The radius of the dome, by default 5.0 + + Returns + ------- + ub : array + Values of the upper bound in the points + lb : array + Values of the lower bound in the points + """ + + xc = center[0] + yc = center[1] + ri = radius - thk / 2 + re = radius + thk / 2 + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + + for i in range(len(x)): + zi2 = ri**2 - (x[i] - xc) ** 2 - (y[i] - yc) ** 2 + ze2 = re**2 - (x[i] - xc) ** 2 - (y[i] - yc) ** 2 + ub[i] = math.sqrt(ze2) + if zi2 > 0.0: + lb[i] = math.sqrt(zi2) + + return ub, lb + + +def dome_bounds_derivatives(x, y, thk, min_lb, center=(5.0, 5.0), radius=5.0): + """Compute sensitivities of upper and lower bounds of the dome based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the dome + min_lb : float + Parameter for lower bound in nodes in the boundary + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + radius : float, optional + The radius of the dome, by default 5.0 + + Returns + ------- + dub : array + Values of the sensitivities of upper bound in the points + dlb : array + Values of the sensitivities of lower bound in the points + """ + + xc = center[0] + yc = center[1] + ri = radius - thk / 2 + re = radius + thk / 2 + dub = zeros((len(x), 1)) + dlb = zeros((len(x), 1)) + dubdx = zeros((len(x), len(x))) + dubdy = zeros((len(x), len(x))) + dlbdx = zeros((len(x), len(x))) + dlbdy = zeros((len(x), len(x))) + + for i in range(len(x)): + zi2 = ri**2 - (x[i] - xc) ** 2 - (y[i] - yc) ** 2 + ze2 = re**2 - (x[i] - xc) ** 2 - (y[i] - yc) ** 2 + ze = math.sqrt(ze2) + dub[i] = 1 / 2 * re / ze + dubdx[i, i] = 1 / 2 / ze * -2 * (x[i] - xc) + dubdy[i, i] = 1 / 2 / ze * -2 * (y[i] - yc) + if zi2 > 0.0: + zi = math.sqrt(zi2) + dlb[i] = -1 / 2 * ri / zi + dlbdx[i, i] = 1 / 2 / zi * -2 * (x[i] - xc) + dlbdy[i, i] = 1 / 2 / zi * -2 * (y[i] - yc) + + return dub, dlb + + +def dome_bound_react(x, y, thk, fixed, center=(5.0, 5.0), radius=5.0): + """Computes the reaction bounds of a dome for a given thickness + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the dome + fixed : list + List with indexes of the fixed vertices + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + radius : float, optional + The radius of the dome, by default 5.0 + + Returns + ------- + b : array + The reaction bounds + """ + + [xc, yc] = center[:2] + b = zeros((len(fixed), 2)) + + for i in range(len(fixed)): + i_ = fixed[i] + theta = math.atan2((y[i_] - yc), (x[i_] - xc)) + x_ = abs(thk / 2 * math.cos(theta)) + y_ = abs(thk / 2 * math.sin(theta)) + b[i, 0] = x_ + b[i, 1] = y_ + + return b + + +def dome_bound_react_derivatives(x, y, thk, fixed, center=(5.0, 5.0), radius=5.0): + """Computes the reaction bounds derivatives of a dome for a given thickness + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the dome + fixed : list + List with indexes of the fixed vertices + center : tuple, optional + x, y coordinates of the center of the dome, by default (5.0, 5.0) + radius : float, optional + The radius of the dome, by default 5.0 + + Returns + ------- + db : array + The sensitivity of the reaction bounds + """ + + [xc, yc] = center[:2] + db = zeros((len(fixed), 2)) + + for i in range(len(fixed)): + i_ = fixed[i] + theta = math.atan2((y[i_] - yc), (x[i_] - xc)) + x_ = abs(1 / 2 * math.cos(theta)) + y_ = abs(1 / 2 * math.sin(theta)) + db[i, :] = [x_, y_] + + return db + + +class DomeEnvelope(ParametricEnvelope): + def __init__( + self, + center: tuple = (5.0, 5.0), + radius: float = 5.0, + thickness: float = 0.50, + min_lb: float = 0.0, + n_hoops: int = 24, + n_parallels: int = 40, + r_oculus: float = 0.0, + **kwargs, + ): + super().__init__(thickness=thickness, **kwargs) + self.center = center + self.radius = radius + self.min_lb = min_lb + self.n_hoops = n_hoops + self.n_parallels = n_parallels + self.r_oculus = r_oculus + + self.update_envelope() # Generate the intra/extra/middle meshes + + @property + def __data__(self): + data = super().__data__ + data["center"] = self.center + data["radius"] = self.radius + data["min_lb"] = self.min_lb + data["n_hoops"] = self.n_hoops + data["n_parallels"] = self.n_parallels + data["r_oculus"] = self.r_oculus + return data + + def __str__(self): + return f"DomeEnvelope(name={self.name})" + + def update_envelope(self): + intrados, extrados, middle = dome_envelope( + center=self.center, + radius=self.radius, + thickness=self.thickness, + min_lb=self.min_lb, + n_hoops=self.n_hoops, + n_parallels=self.n_parallels, + r_oculus=self.r_oculus, + ) + self.intrados = intrados + self.extrados = extrados + self.middle = middle + + def compute_middle(self, x, y): + return dome_middle(x, y, self.radius, self.min_lb, self.center) + + def compute_bounds(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return dome_bounds(x, y, thickness, self.min_lb, self.center, self.radius) + + def compute_bounds_derivatives(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return dome_bounds_derivatives(x, y, thickness, self.min_lb, self.center, self.radius) + + def compute_bound_react(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return dome_bound_react(x, y, thickness, fixed, self.center, self.radius) + + def compute_bound_react_derivatives(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return dome_bound_react_derivatives(x, y, thickness, fixed, self.center, self.radius) diff --git a/src/compas_tna/envelope/envelope.py b/src/compas_tna/envelope/envelope.py new file mode 100644 index 00000000..d39afc3e --- /dev/null +++ b/src/compas_tna/envelope/envelope.py @@ -0,0 +1,126 @@ +from typing import Optional +from typing import Tuple + +import numpy as np + +from compas.data import Data +from compas_tna.diagrams import FormDiagram + + +class Envelope(Data): + """Pure geometric envelope representing masonry structure boundaries.""" + + def __init__(self, rho: Optional[float] = 20.0, is_parametric: bool = False, **kwargs): + super().__init__(**kwargs) + + self.rho = rho + self.is_parametric = is_parametric + + # Computed properties (cached) + self._area = 0.0 + self._volume = 0.0 + self._total_selfweight = 0.0 + + @property + def __data__(self): + data = {} + data["rho"] = self.rho + data["is_parametric"] = self.is_parametric + return data + + def __str__(self): + return f"Envelope(name={self.name})" + + # ============================================================================= + # Properties + # ============================================================================= + + @property + def area(self): + if not self._area: + self._area = self.compute_area() + return self._area + + @property + def volume(self): + if not self._volume: + self._volume = self.compute_volume() + return self._volume + + @property + def selfweight(self): + if not self._total_selfweight: + self._total_selfweight = self.compute_selfweight() + return self._total_selfweight + + # ============================================================================= + # Geometry operations + # ============================================================================= + + def compute_area(self) -> float: + """Compute and returns the total area of the structure based on the appropriate method.""" + + raise NotImplementedError("Implement compute_area for specific envelope type.") + + def compute_volume(self) -> float: + """Compute and returns the total volume of the structure based on the appropriate method.""" + + raise NotImplementedError("Implement compute_volume for specific envelope type.") + + def compute_selfweight(self) -> float: + """Compute and returns the total selfweight of the structure based on the appropriate method.""" + + raise NotImplementedError("Implement compute_selfweight for specific envelope type.") + + # ============================================================================= + # TNA-related operations (accept formdiagram) + # ============================================================================= + + def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply selfweight to the nodes of a form diagram based on the appropriate method.""" + + raise NotImplementedError("Implement apply_selfweight_to_formdiagram for specific envelope type.") + + def apply_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply envelope bounds to a form diagram based on the appropriate method.""" + + raise NotImplementedError("Implement apply_envelope_to_formdiagram for specific envelope type.") + + def apply_target_heights_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply target heights to a form diagram based on the appropriate method.""" + + raise NotImplementedError("Implement apply_target_heights_to_formdiagram for specific envelope type.") + + def apply_reaction_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply reaction bounds to a form diagram based on the appropriate method.""" + + raise NotImplementedError("Implement apply_reaction_bounds_to_formdiagram for specific envelope type.") + + # ============================================================================= + # Numerical methods to be implemented at sub classes + # ============================================================================= + + def compute_middle(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: + """Compute the middle of the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement compute_middle for specific envelope type.") + + def compute_bounds(self, x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Compute the upper and lower bounds of the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement compute_bounds for specific envelope type.") + + def compute_bounds_derivatives(self, x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Compute the upper and lower bounds derivativesof the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement compute_bounds_derivatives for specific envelope type.") + + def compute_bound_react(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: + """Compute the reaction bounds of the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement compute_bound_react for specific envelope type.") + + def compute_bound_react_derivatives(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: + """Compute the bound_react_derivatives of the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement compute_bound_react_derivatives for specific envelope type.") diff --git a/src/compas_tna/envelope/meshenvelope.py b/src/compas_tna/envelope/meshenvelope.py new file mode 100644 index 00000000..da3898de --- /dev/null +++ b/src/compas_tna/envelope/meshenvelope.py @@ -0,0 +1,657 @@ +import math +from typing import Optional + +from numpy import asarray +from scipy.interpolate import griddata + +from compas.datastructures import Mesh +from compas_tna.diagrams import FormDiagram +from compas_tna.envelope import Envelope + + +# TODO: What if intrados and extrados are surfaces? +def interpolate_middle_mesh(intrados: Mesh, extrados: Mesh) -> Mesh: + """Interpolate a middle mesh between intrados and extrados meshes. + + This function properly calculates thickness by considering the normal vector + at each point, ensuring accurate thickness measurements on curved surfaces. + + Parameters + ---------- + intrados : Mesh + The intrados surface mesh. + extrados : Mesh + The extrados surface mesh. + + Returns + ------- + Mesh + The interpolated middle mesh with proper normal-based thickness stored. + """ + # Use the intrados as base topology + middle = intrados.copy() + + # Get point clouds for interpolation + intrados_points = asarray(intrados.vertices_attributes("xyz")) + extrados_points = asarray(extrados.vertices_attributes("xyz")) + + # Get XY coordinates of middle mesh + middle_xy = asarray(middle.vertices_attributes("xy")) + + # Interpolate Z coordinates from both surfaces + zi = griddata(intrados_points[:, :2], intrados_points[:, 2], middle_xy, method="linear") + ze = griddata(extrados_points[:, :2], extrados_points[:, 2], middle_xy, method="linear") + + # First loop: set middle Z as average + for i, key in enumerate(middle.vertices()): + middle_z = (zi[i] + ze[i]) / 2.0 + middle.vertex_attribute(key, "z", middle_z) + + # Second loop: calculate and set thickness using correct normals + for i, key in enumerate(middle.vertices()): + nx, ny, nz = middle.vertex_normal(key) + z_diff = ze[i] - zi[i] + if abs(nz) > 0.1: + thickness = abs(z_diff) * abs(nz) + else: + thickness = abs(z_diff) + middle.vertex_attribute(key, "thickness", thickness) + + return middle + + +# TODO: What if middle is a surface and not a mesh? +def offset_from_middle(middle: Mesh, fixed_xy: bool = True) -> tuple[Mesh, Mesh]: + """ + Offset a middle surface mesh to obtain extrados and intrados meshes using thickness attributes. + + This function properly handles curved surfaces by using the normal vector + and the thickness measured perpendicular to the surface. + + Parameters + ---------- + middle : Mesh + The middle surface mesh with thickness attributes per vertex. + fixed_xy : bool, optional + If True, extrados/intrados will have the same XY as the middle mesh, + and only Z will be offset (with normal correction). + If False, full 3D normal offset is used. + + Returns + ------- + tuple[Mesh, Mesh] + (intrados, extrados) offset meshes. + """ + extrados = middle.copy() + intrados = middle.copy() + + for key in middle.vertices(): + x, y, z = middle.vertex_coordinates(key) + nx, ny, nz = middle.vertex_normal(key) + + # Get thickness for this specific vertex (should be normal-based) + thickness = middle.vertex_attribute(key, "thickness") + if thickness is None: + thickness = 0.5 + half_thick = 0.5 * thickness + + if fixed_xy: + # Prevent division by zero for horizontal normals + if abs(nz) < 1e-8: + raise ValueError(f"Normal at vertex {key} is (almost) horizontal: {nx, ny, nz}") + dz = half_thick / nz + extrados_z = z + dz + intrados_z = z - dz + extrados.vertex_attribute(key, "z", extrados_z) + intrados.vertex_attribute(key, "z", intrados_z) + else: + # Full 3D normal offset - this is the most accurate for curved surfaces + extrados.vertex_attributes( + key, + "xyz", + [x + half_thick * nx, y + half_thick * ny, z + half_thick * nz], + ) + intrados.vertex_attributes( + key, + "xyz", + [x - half_thick * nx, y - half_thick * ny, z - half_thick * nz], + ) + + return intrados, extrados + + +# TODO: What if the target is a surface and not a mesh? +def project_mesh_to_target_vertical(mesh: Mesh, target: Mesh) -> None: + """Project a mesh vertically (in Z direction) onto a target mesh. + + Parameters + ---------- + mesh : Mesh + The mesh to be projected. + target : Mesh + The target mesh to project onto. + + Returns + ------- + None + The mesh is modified in place. + """ + # Get target mesh vertices for simple vertical projection + target_vertices = list(target.vertices()) + target_points = [target.vertex_point(v) for v in target_vertices] + + for vertex in mesh.vertices(): + point = mesh.vertex_point(vertex) + + # Find the closest target vertex in XY plane + min_distance = float("inf") + closest_z = point.z + + for target_point in target_points: + # Calculate XY distance (ignore Z) + xy_distance = ((point.x - target_point.x) ** 2 + (point.y - target_point.y) ** 2) ** 0.5 + + if xy_distance < min_distance: + min_distance = xy_distance + closest_z = target_point.z + + # Update vertex to closest Z value + new_point = point.copy() + new_point.z = closest_z + mesh.vertex_attributes(vertex, "xyz", new_point) + + +def pattern_inverse_height_thickness(pattern: Mesh, tmin: Optional[float] = None, tmax: Optional[float] = None) -> None: + """Set variable thickness based on inverse height. + + Parameters + ---------- + pattern : Mesh + The mesh to apply thickness to. + tmin : float, optional + Minimum thickness. If None, will be calculated as 3/1000 of the diagonal of the xy bounding box. + tmax : float, optional + Maximum thickness. If None, will be calculated as 50/1000 of the diagonal of the xy bounding box. + """ + x: list[float] = pattern.vertices_attribute(name="x") + xmin = min(x) + xmax = max(x) + dx = xmax - xmin + + y: list[float] = pattern.vertices_attribute(name="y") + ymin = min(y) + ymax = max(y) + dy = ymax - ymin + + d = (dx**2 + dy**2) ** 0.5 + + tmin = tmin or 3 * d / 1000 + tmax = tmax or 50 * d / 1000 + + pattern.update_default_vertex_attributes(thickness=0) + zvalues: list[float] = pattern.vertices_attribute(name="z") + zmin = min(zvalues) + zmax = max(zvalues) + + for vertex in pattern.vertices(): + point = pattern.vertex_point(vertex) + z = (point.z - zmin) / (zmax - zmin) + thickness = (1 - z) * (tmax - tmin) + tmin + pattern.vertex_attribute(vertex, name="thickness", value=thickness) + + +class MeshEnvelope(Envelope): + """An Envelope defined by meshes at intrados and extrados.""" + + def __init__(self, intrados: Mesh = None, extrados: Mesh = None, middle: Mesh = None, fill: Mesh = None, thickness: float = 0.5, **kwargs): + super().__init__(**kwargs) + + self.intrados = intrados + self.extrados = extrados + self.middle = middle + self.fill = fill + + # Thickness property + self._thickness = thickness + + @property + def __data__(self): + data = super().__data__ + data["intrados"] = self.intrados + data["extrados"] = self.extrados + data["middle"] = self.middle + data["fill"] = self.fill + data["thickness"] = self._thickness + return data + + def __str__(self): + return f"Envelope(name={self.name})" + + # ============================================================================= + # Factory methods + # ============================================================================= + + @classmethod + def from_meshes(cls, intrados: Mesh, extrados: Mesh, middle: Optional[Mesh] = None) -> "Envelope": + """Construct an envelope from intrados and extrados meshes. + + Parameters + ---------- + intrados : Mesh + The intrados surface mesh of the envelope. + extrados : Mesh + The extrados surface mesh of the envelope. + middle : Mesh, optional + The middle surface mesh of the envelope. + + Returns + ------- + :class:`Envelope` + + """ + envelope = cls() + envelope.intrados = intrados + envelope.extrados = extrados + if middle is not None: + envelope.middle = middle + else: + envelope.middle = interpolate_middle_mesh(intrados, extrados) + + return envelope + + @classmethod + def from_middle_mesh(cls, mesh: Mesh, thickness: Optional[float] = None) -> "Envelope": + """Construct an envelope from a mesh with specified thickness. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to create the envelope from. + thickness : float, optional + The thickness of the envelope. If None, uses thickness values stored in formdiagram vertices. + + Returns + ------- + :class:`Envelope` + + """ + envelope = cls() + + envelope.middle = mesh.copy(cls=Mesh) + + if thickness is not None: + envelope.thickness = thickness + + # Create intrados and extrados using thickness from middle mesh + intrados, extrados = offset_from_middle(envelope.middle) + envelope.intrados = intrados + envelope.extrados = extrados + + return envelope + + @classmethod + def from_formdiagram(cls, formdiagram: FormDiagram, thickness: Optional[float] = None) -> "Envelope": + """Construct an envelope from a FormDiagram with specified thickness. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to create the envelope from. + thickness : float, optional + The thickness of the envelope. If None, uses thickness values stored in formdiagram vertices. + + Returns + ------- + :class:`Envelope` + + """ + return cls.from_middle_mesh(formdiagram, thickness) + + # ============================================================================= + # Properties + # ============================================================================= + + @property + def thickness(self) -> float: + """Get the average thickness of the envelope. + + Returns + ------- + float + The average thickness of the envelope. + """ + if self.middle is not None: + # Return average thickness from middle mesh vertices + thicknesses = [] + for key in self.middle.vertices(): + thickness = self.middle.vertex_attribute(key, "thickness") + if thickness is not None: + thicknesses.append(thickness) + + if thicknesses: + return sum(thicknesses) / len(thicknesses) + + return self._thickness + + @thickness.setter + def thickness(self, value: float) -> None: + """Set a uniform thickness for all vertices of the envelope. + + Parameters + ------- + value : float + The thickness value to set for all vertices. + """ + self._thickness = value + + # Update middle mesh if it exists + if self.middle is not None: + for key in self.middle.vertices(): + self.middle.vertex_attribute(key, "thickness", value) + + @property + def is_parametric(self) -> bool: + """Check if the envelope is parametric.""" + return False + + # ============================================================================= + # Geometric operations + # ============================================================================= + + def set_variable_thickness(self, tmin: Optional[float] = None, tmax: Optional[float] = None) -> None: + """Set variable thickness based on inverse height using the pattern_inverse_height_thickness function. + + This method applies thickness variation based on the height of vertices in the middle mesh, + where higher vertices get thinner thickness and lower vertices get thicker thickness. + + Parameters + ------- + tmin : float, optional + Minimum thickness. If None, will be calculated as 3/1000 of the diagonal of the xy bounding box. + tmax : float, optional + Maximum thickness. If None, will be calculated as 50/1000 of the diagonal of the xy bounding box. + """ + if self.middle is None: + raise ValueError("Middle mesh is not available. Cannot set variable thickness.") + + # Apply the pattern_inverse_height_thickness function to the middle mesh + pattern_inverse_height_thickness(self.middle, tmin=tmin, tmax=tmax) + + # Update intrados and extrados meshes + self.intrados, self.extrados = offset_from_middle(self.middle) + + def compute_volume(self) -> float: + """Compute and returns the volume of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total volume of the structure. + + """ + if self.middle is None: + raise ValueError("Middle mesh is not available. Cannot compute volume.") + + middle = self.middle + total_volume = 0.0 + + # Use variable thickness from middle mesh vertices + for vertex in middle.vertices(): + thickness = middle.vertex_attribute(vertex, "thickness") + if thickness is None: + thickness = self._thickness + vertex_area = middle.vertex_area(vertex) # should be projected area + vertex_volume = thickness * vertex_area + total_volume += vertex_volume + + return total_volume + + def compute_selfweight(self) -> float: + """Compute and returns the total selfweight of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total selfweight of the structure. + + """ + if self.middle is None: + if self.intrados is not None and self.extrados is not None: + self.middle = interpolate_middle_mesh(self.intrados, self.extrados) + else: + raise ValueError("Middle mesh is not available and cannot be interpolated.") + + middle = self.middle + rho = self.rho + total_selfweight = 0.0 + + # Use variable thickness from middle mesh vertices + for vertex in middle.vertices(): + thickness = middle.vertex_attribute(vertex, "thickness") + if thickness is None: + thickness = self._thickness + vertex_area = middle.vertex_area(vertex) + vertex_volume = thickness * vertex_area + vertex_weight = vertex_volume * rho + total_selfweight += vertex_weight + + return total_selfweight + + def compute_area(self) -> float: + """Compute and returns the total selfweight of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total selfweight of the structure. + """ + if self.middle is None: + raise ValueError("Middle mesh is not available. Cannot compute area.") + + return self.middle.area() + + # ============================================================================= + # TNA-specific operations (accept formdiagram as parameter) + # ============================================================================= + + def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram, normalize=True) -> None: + """Apply selfweight to the nodes of a form diagram based on the middle surface and local thicknesses. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to apply selfweight to. + normalize : bool, optional + Whether or not normalize the selfweight to match the computed total selfweight, by default True + + Returns + ------- + None + The FormDiagram is modified in place + + """ + # Step 1: Check that middle mesh is present + if self.middle is None: + if self.intrados is not None and self.extrados is not None: + self.middle = interpolate_middle_mesh(self.intrados, self.extrados) + else: + raise ValueError("Middle mesh is not set. Please set the middle mesh before applying selfweight.") + + # Step 2: Compute the selfweight of the shell + total_selfweight = self.compute_selfweight() + + # Step 3: Sync thickness to the form diagram + self.sync_thickness_to_formdiagram(formdiagram) + + # Step 4: Copy the form diagram and project it onto the middle mesh vertically + form_ = formdiagram.copy() + project_mesh_to_target_vertical(form_, self.middle) + + # Step 5: Compute and lump selfweight at vertices + total_pz = 0.0 + for vertex in form_.vertices(): + # Get vertex area and thickness + vertex_area = form_.vertex_area(vertex) + thickness = form_.vertex_attribute(vertex, "thickness") + + # Compute selfweight contribution (negative for downward direction) + pz = -vertex_area * thickness * self.rho + + # Store in form diagram + formdiagram.vertex_attribute(vertex, "pz", pz) + total_pz += abs(pz) # Sum absolute values for normalization + + # Step 6: Scale to match total selfweight if normalize=True + if normalize and total_pz > 0: + scale_factor = total_selfweight / total_pz + if scale_factor != 1.0: + print(f"Scaled selfweight by factor: {scale_factor}") + + for vertex in formdiagram.vertices(): + pz = formdiagram.vertex_attribute(vertex, "pz") + formdiagram.vertex_attribute(vertex, "pz", pz * scale_factor) + + print(f"Selfweight applied to form diagram. Total load: {sum(abs(formdiagram.vertex_attribute(vertex, 'pz')) for vertex in formdiagram.vertices())}") + + def apply_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply envelope bounds to a form diagram based on the intrados and extrados surfaces. + + This method projects the form diagram onto both intrados and extrados surfaces + and assigns the heights to 'ub' (upper bound) and 'lb' (lower bound) properties. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to apply bounds to. + + Returns + ------- + None + The FormDiagram is modified in place. + """ + # Step 1: Check that intrados and extrados are present + if self.intrados is None or self.extrados is None: + raise ValueError("Intra/Extrados not set. Please set them before applying bounds.") + + # Step 2: Copy the form diagram for projection + form_ub = formdiagram.copy() # For upper bound (extrados) + form_lb = formdiagram.copy() # For lower bound (intrados) + + # Step 3: Project form diagram onto extrados (upper bound) + project_mesh_to_target_vertical(form_ub, self.extrados) + project_mesh_to_target_vertical(form_lb, self.intrados) + + # Step 4: Collect heights and assign to form diagram + for vertex in formdiagram.vertices(): + if vertex in form_ub.vertices() and vertex in form_lb.vertices(): + # Get z coordinates from projected meshes + _, _, z_ub = form_ub.vertex_coordinates(vertex) + _, _, z_lb = form_lb.vertex_coordinates(vertex) + + # Assign to form diagram + formdiagram.vertex_attribute(vertex, "ub", z_ub) + formdiagram.vertex_attribute(vertex, "lb", z_lb) + else: + print(f"Warning: Vertex {vertex} not found in projected meshes") + # Set default values if vertex not found + formdiagram.vertex_attribute(vertex, "ub", float("inf")) + formdiagram.vertex_attribute(vertex, "lb", float("-inf")) + + def apply_target_heights_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply target heights to a form diagram based on the Envelope middle surface. + + This method projects the form diagram onto the Envelope middle surface + and assigns the heights to 'target' property. This assignment can later be used to compute a bestfit optimization. + """ + if self.middle is None: + raise ValueError("Middle mesh is not set. Please set the middle mesh before applying target heights.") + + # Step 1: Copy the form diagram for projection + form_target = formdiagram.copy() # For upper bound (extrados) + + project_mesh_to_target_vertical(form_target, self.middle) + + # Step 2: Collect heights and assign to form diagram + for vertex in formdiagram.vertices(): + if vertex in form_target.vertices(): + z_target = form_target.vertex_attribute(vertex, "z") + formdiagram.vertex_attribute(vertex, "target", z_target) + + def apply_reaction_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply reaction bounds to a form diagram based on the Envelope middle surface. + + This method projects the form diagram onto the Envelope middle surface + and assigns the heights to 'target' property. This assignment can later be used to compute a bestfit optimization. + """ + + # if not self.callable_bound_react: + # raise ValueError("Callable bound reaction is not set. Please set this limit manually.") + + ## TODO: Implement this + + def sync_thickness_to_formdiagram(self, formdiagram: FormDiagram, method="linear", extrapolate=True) -> None: + """Synchronize thickness attributes from middle mesh to form diagram using continuous interpolation. + + This method creates a continuous thickness map from the middle mesh and interpolates + thickness values at the form diagram vertex locations. This ensures compatibility + even when the middle mesh and form diagram have different topologies. + + Parameters + ------- + formdiagram : FormDiagram + The form diagram to sync thickness to. + """ + if self.middle is None: + raise ValueError("Middle mesh must be set to sync thickness.") + + # Get middle mesh XY coordinates and thickness values + middle_xy = self.middle.vertices_attributes("xy") + middle_thickness = self.middle.vertices_attribute("thickness") + + # Validate data + if not middle_xy or not middle_thickness: + raise ValueError("Middle mesh must have both 'xy' and 'thickness' attributes.") + + # Convert to numpy arrays + middle_xy_array = asarray(middle_xy) + middle_thickness_array = asarray(middle_thickness) + + # Get form diagram XY coordinates + form_xy = formdiagram.vertices_attributes("xy") + if not form_xy: + raise ValueError("Form diagram must have 'xy' attributes.") + + form_xy_array = asarray(form_xy) + + # Use griddata for 2D interpolation + interpolated_thickness = griddata( + middle_xy_array, + middle_thickness_array, + form_xy_array, + method=method, + fill_value=self._thickness, # Use default thickness for points outside convex hull + # extrapolate=extrapolate + ) + + # Assign interpolated thickness values to form diagram vertices + for i, vertex in enumerate(formdiagram.vertices()): + thickness_value = float(interpolated_thickness[i]) + # Ensure thickness is positive and reasonable + if thickness_value <= 0 or math.isnan(thickness_value): + thickness_value = self._thickness + formdiagram.vertex_attribute(vertex, "thickness", thickness_value) + + def compute_middle(self, x, y): + raise NotImplementedError("Implement compute_middle for specific envelope type.") + + def compute_bounds(self, x, y, thickness): + raise NotImplementedError("Implement compute_bounds for specific envelope type.") + + def compute_bounds_derivatives(self, x, y): + raise NotImplementedError("Implement compute_bounds_derivatives for specific envelope type.") + + def compute_bound_react(self, x, y, thickness, fixed): + raise NotImplementedError("Implement compute_bound_react for specific envelope type.") + + def compute_bound_react_derivatives(self, x, y, thickness, fixed): + raise NotImplementedError("Implement compute_bound_react_derivatives for specific envelope type.") diff --git a/src/compas_tna/envelope/parametricenvelope.py b/src/compas_tna/envelope/parametricenvelope.py new file mode 100644 index 00000000..3faf93c7 --- /dev/null +++ b/src/compas_tna/envelope/parametricenvelope.py @@ -0,0 +1,216 @@ +import numpy as np + +from compas_tna.diagrams import FormDiagram +from compas_tna.envelope import Envelope + + +class ParametricEnvelope(Envelope): + """Pure geometric envelope representing masonry structure boundaries created parametrically.""" + + def __init__(self, thickness: float = 0.50, **kwargs): + super().__init__(**kwargs) + + self.is_parametric = True + self._thickness = thickness + + def __str__(self): + return f"ParametricEnvelope(name={self.name})" + + @property + def __data__(self): + data = super().__data__ + data["thickness"] = self._thickness + return data + + # ============================================================================= + # Parametric Common Prepoperties + # ============================================================================= + + @property + def thickness(self) -> float: + """Get the thickness of the envelope. + + Returns + ------- + float + The thickness of the envelope. + """ + return self._thickness + + @thickness.setter + def thickness(self, value: float) -> None: + """Set the thickness of the envelope. + + Parameters + ------- + value : float + The thickness value to set. + """ + self._thickness = value + + # ============================================================================= + # Envelope Generator + # ============================================================================= + + def update_envelope(self) -> None: + """Update the envelope based on the appropriate method.""" + + raise NotImplementedError("Implement update_envelope for specific envelope type.") + + # ============================================================================= + # Geometry operations + # ============================================================================= + + def compute_volume(self) -> float: + """Compute and returns the volume of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total volume of the structure. + + """ + if self.middle is None: + self.update_envelope() + + return self.compute_area() * self.thickness + + def compute_selfweight(self) -> float: + """Compute and returns the total selfweight of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total selfweight of the structure. + + """ + if self.middle is None: + self.update_envelope() + + return self.compute_volume() * self.rho + + def compute_area(self) -> float: + """Compute and returns the total selfweight of the structure based on the area and thickness in the data. + + Returns + ------- + float + The total selfweight of the structure. + """ + if self.middle is None: + self.update_envelope() + + return self.middle.area() + + # ============================================================================= + # TNA-specific operations (accept formdiagram as parameter) + # ============================================================================= + + def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram, normalize=True) -> None: + """Apply selfweight to the nodes of a form diagram based on the middle surface and local thicknesses. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to apply selfweight to. + normalize : bool, optional + Whether or not normalize the selfweight to match the computed total selfweight, by default True + + Returns + ------- + None + The FormDiagram is modified in place + + """ + # Step 2: Compute the selfweight of the shell + total_selfweight = self.compute_selfweight() + + # Step 3: Copy the form diagram and project it onto the middle mesh vertically + form_: FormDiagram = formdiagram.copy() + + xy = np.array(form_.vertices_attributes("xy")) + zt = list(self.compute_middle(xy[:, 0], xy[:, 1]).flatten().tolist()) + for i, key in enumerate(form_.vertices()): + form_.vertex_attribute(key, "z", zt[i]) + + # Step 4: Compute and lump selfweight at vertices + total_pz = 0.0 + for vertex in form_.vertices(): + vertex_area = form_.vertex_area(vertex) + thickness = self.thickness + pz = -vertex_area * thickness * self.rho + formdiagram.vertex_attribute(vertex, "pz", pz) + total_pz += abs(pz) + + # Step 5: Scale to match total selfweight if normalize=True + if normalize and total_pz > 0: + scale_factor = total_selfweight / total_pz + if scale_factor != 1.0: + print(f"Scaled selfweight by factor: {scale_factor}") + + for vertex in formdiagram.vertices(): + pz = formdiagram.vertex_attribute(vertex, "pz") + formdiagram.vertex_attribute(vertex, "pz", pz * scale_factor) + + print(f"Selfweight applied to form diagram. Total load: {sum(abs(formdiagram.vertex_attribute(vertex, 'pz')) for vertex in formdiagram.vertices())}") + + def apply_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply envelope bounds to a form diagram based on the intrados and extrados surfaces. + + This method projects the form diagram onto both intrados and extrados surfaces + and assigns the heights to 'ub' (upper bound) and 'lb' (lower bound) properties. + + Parameters + ---------- + formdiagram : FormDiagram + The form diagram to apply bounds to. + + Returns + ------- + None + The FormDiagram is modified in place. + """ + + xy = np.array(formdiagram.vertices_attributes("xy")) + zub, zlb = self.compute_bounds(xy[:, 0], xy[:, 1]) + for i, key in enumerate(formdiagram.vertices()): + formdiagram.vertex_attribute(key, "ub", float(zub[i])) + formdiagram.vertex_attribute(key, "lb", float(zlb[i])) + # TODO: Future Cached properties could be added here + + def apply_target_heights_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply target heights to a form diagram based on the Envelope middle surface. + + This method projects the form diagram onto the Envelope middle surface + and assigns the heights to 'target' property. This assignment can later be used to compute a bestfit optimization. + """ + + xy = np.array(formdiagram.vertices_attributes("xy")) + zt = list(self.compute_middle(xy[:, 0], xy[:, 1]).flatten().tolist()) + for i, key in enumerate(formdiagram.vertices()): + formdiagram.vertex_attribute(key, "z", zt[i]) + # TODO: Future Cached properties could be added here + + def apply_reaction_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None: + """Apply reaction bounds to a form diagram based on the Envelope middle surface. + + This method projects the form diagram onto the Envelope middle surface + and assigns the heights to 'target' property. This assignment can later be used to compute a bestfit optimization. + """ + raise NotImplementedError("Implement apply_reaction_bounds_to_formdiagram for specific envelope type.") + ## TODO: Implement this + + def compute_middle(self, x, y): + raise NotImplementedError("Implement compute_middle for specific envelope type.") + + def compute_bounds(self, x, y, thickness): + raise NotImplementedError("Implement compute_bounds for specific envelope type.") + + def compute_bounds_derivatives(self, x, y): + raise NotImplementedError("Implement compute_bounds_derivatives for specific envelope type.") + + def compute_bound_react(self, x, y, thickness, fixed): + raise NotImplementedError("Implement compute_bound_react for specific envelope type.") + + def compute_bound_react_derivatives(self, x, y, thickness, fixed): + raise NotImplementedError("Implement compute_bound_react_derivatives for specific envelope type.") diff --git a/src/compas_tna/envelope/pavillionvault.py b/src/compas_tna/envelope/pavillionvault.py new file mode 100644 index 00000000..58ab0fbc --- /dev/null +++ b/src/compas_tna/envelope/pavillionvault.py @@ -0,0 +1,462 @@ +import math + +from numpy import array +from numpy import ones +from numpy import zeros + +from compas.datastructures import Mesh +from compas_tna.diagrams.diagram_rectangular import create_cross_mesh +from compas_tna.envelope.parametricenvelope import ParametricEnvelope + + +def pavillionvault_envelope( + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, + spr_angle: float = 0.0, +): + """Create an envelope for a pavillion vault geometry with given parameters. + + Parameters + ---------- + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + thickness : float, optional + Thickness of the vault, by default 0.50 + min_lb : float, optional + Parameter for lower bound in nodes in the boundary, by default 0.0 + n : int, optional + Number of vertices for the mesh, by default 100 + spr_angle : float, optional + Springing angle, by default 0.0 + + Returns + ------- + intrados : Mesh + Intrados mesh + extrados : Mesh + Extrados mesh + middle : Mesh + Middle mesh + """ + # Create base topology + base_topology = create_cross_mesh(x_span=x_span, y_span=y_span, n=n) + xyz0, faces_i = base_topology.to_vertices_and_faces() + xi, yi, _ = array(xyz0).transpose() + + # Create middle surface + zt = pavillionvault_middle(xi, yi, x_span=x_span, y_span=y_span, spr_angle=spr_angle, tol=1e-6) + xyzt = array([xi, yi, zt.flatten()]).transpose() + middle = Mesh.from_vertices_and_faces(xyzt, faces_i) + middle.update_default_vertex_attributes(thickness=thickness) + + # Create upper and lower bounds + zub, zlb = pavillionvault_bounds(xi, yi, thickness, min_lb, x_span=x_span, y_span=y_span, spr_angle=spr_angle, tol=1e-6) + xyzlb = array([xi, yi, zlb.flatten()]).transpose() + intrados = Mesh.from_vertices_and_faces(xyzlb, faces_i) + + x_span_extra = (x_span[0] - thickness / 2 / math.cos(math.radians(spr_angle)), x_span[1] + thickness / 2 / math.cos(math.radians(spr_angle))) + y_span_extra = (y_span[0] - thickness / 2 / math.cos(math.radians(spr_angle)), y_span[1] + thickness / 2 / math.cos(math.radians(spr_angle))) + extra_topology = create_cross_mesh(x_span=x_span_extra, y_span=y_span_extra, n=n) + xyz0, faces_i = extra_topology.to_vertices_and_faces() + xi, yi, _ = array(xyz0).transpose() + zub, zlb = pavillionvault_bounds(xi, yi, thickness, min_lb, x_span=x_span, y_span=y_span, spr_angle=spr_angle, tol=1e-6) + xyzub = array([xi, yi, zub.flatten()]).transpose() + extrados = Mesh.from_vertices_and_faces(xyzub, faces_i) + + return intrados, extrados, middle + + +def pavillionvault_middle(x, y, x_span=(0.0, 10.0), y_span=(0.0, 10.0), spr_angle=0.0, tol=1e-6): + """Compute middle of a pavillion vault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + spr_angle : float, optional + Springing angle, by default 0.0 + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + z : array + Values of the middle surface in the points + """ + + x0, x1 = x_span + y0, y1 = y_span + + if spr_angle == 0.0: + z_ = 0.0 + else: + alpha = 1 / math.cos(math.radians(spr_angle)) + z_ = (x1 - x0) / 2 * math.tan(math.radians(spr_angle)) + L = x1 * alpha + Ldiff = L - x1 + x0, x1 = -Ldiff / 2, x1 + Ldiff / 2 + y0, y1 = -Ldiff / 2, y1 + Ldiff / 2 + + rx = (x1 - x0) / 2 + ry = (y1 - y0) / 2 + + z = zeros((len(x), 1)) + + for i in range(len(x)): + xi, yi = x[i], y[i] + if (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q1 + z[i] = math.sqrt((ry) ** 2 - ((yi - y0) - ry) ** 2) - z_ + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q3 + z[i] = math.sqrt((ry) ** 2 - ((yi - y0) - ry) ** 2) - z_ + elif (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q2 + z[i] = math.sqrt((rx) ** 2 - ((xi - x0) - rx) ** 2) - z_ + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q4 + z[i] = math.sqrt((rx) ** 2 - ((xi - x0) - rx) ** 2) - z_ + else: + print("Error Q. (x,y) = ({0},{1})".format(xi, yi)) + + return z + + +def pavillionvault_bounds(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), spr_angle=0.0, tol=1e-6): + """Compute upper and lower bounds of a pavillionvault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the vault + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + spr_angle : float, optional + Springing angle, by default 0.0 + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + ub : array + Values of the upper bound in the points + lb : array + Values of the lower bound in the points + """ + + x0, x1 = x_span + y0, y1 = y_span + + if spr_angle == 0.0: + z_ = 0.0 + else: + alpha = 1 / math.cos(math.radians(spr_angle)) + z_ = (x1 - x0) / 2 * math.tan(math.radians(spr_angle)) + L = x1 * alpha + Ldiff = L - x1 + x0, x1 = -Ldiff / 2, x1 + Ldiff / 2 + y0, y1 = -Ldiff / 2, y1 + Ldiff / 2 + + y1_ub = y1 + thk / 2 + y0_ub = y0 - thk / 2 + x1_ub = x1 + thk / 2 + x0_ub = x0 - thk / 2 + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + rx_ub = (x1_ub - x0_ub) / 2 + ry_ub = (y1_ub - y0_ub) / 2 + rx_lb = (x1_lb - x0_lb) / 2 + ry_lb = (y1_lb - y0_lb) / 2 + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + + for i in range(len(x)): + xi, yi = x[i], y[i] + intrados_null = False + if yi > y1_lb or xi > x1_lb or xi < x0_lb or yi < y0_lb: + intrados_null = True + if (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q1 + ub[i] = math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2) - z_ + if not intrados_null: + lb[i] = math.sqrt((ry_lb) ** 2 - ((yi - y0_lb) - ry_lb) ** 2) - z_ + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q3 + ub[i] = math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2) - z_ + if not intrados_null: + lb[i] = math.sqrt((ry_lb) ** 2 - ((yi - y0_lb) - ry_lb) ** 2) - z_ + elif (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q2 + ub[i] = math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2) - z_ + if not intrados_null: + lb[i] = math.sqrt((rx_lb) ** 2 - ((xi - x0_lb) - rx_lb) ** 2) - z_ + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q4 + ub[i] = math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2) - z_ + if not intrados_null: + lb[i] = math.sqrt((rx_lb) ** 2 - ((xi - x0_lb) - rx_lb) ** 2) - z_ + else: + print("Error Q. (x,y) = ({0},{1})".format(xi, yi)) + + return ub, lb + + +def pavillionvault_bounds_derivatives(x, y, thk, min_lb, x_span=(0.0, 10.0), y_span=(0.0, 10.0), tol=1e-6): + """Computes the sensitivities of upper and lower bounds in the x, y coordinates and thickness specified. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the vault + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + dub : array + Values of the sensitivities for the upper bound in the points + dlb : array + Values of the sensitivities for the lower bound in the points + """ + + x0, x1 = x_span + y0, y1 = y_span + + y1_ub = y1 + thk / 2 + y0_ub = y0 - thk / 2 + x1_ub = x1 + thk / 2 + x0_ub = x0 - thk / 2 + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + rx_ub = (x1_ub - x0_ub) / 2 + ry_ub = (y1_ub - y0_ub) / 2 + rx_lb = (x1_lb - x0_lb) / 2 + ry_lb = (y1_lb - y0_lb) / 2 + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + dub = zeros((len(x), 1)) + dlb = zeros((len(x), 1)) + + for i in range(len(x)): + xi, yi = x[i], y[i] + intrados_null = False + if yi > y1_lb or xi > x1_lb or xi < x0_lb or yi < y0_lb: + intrados_null = True + if (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q1 + ub[i] = math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2) + dub[i] = 1 / 2 * ry_ub / ub[i] + if not intrados_null: + lb[i] = math.sqrt((ry_lb) ** 2 - ((yi - y0_lb) - ry_lb) ** 2) + dlb[i] = -1 / 2 * ry_lb / lb[i] + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q3 + ub[i] = math.sqrt((ry_ub) ** 2 - ((yi - y0_ub) - ry_ub) ** 2) + dub[i] = 1 / 2 * ry_ub / ub[i] + if not intrados_null: + lb[i] = math.sqrt((ry_lb) ** 2 - ((yi - y0_lb) - ry_lb) ** 2) + dlb[i] = -1 / 2 * ry_lb / lb[i] + elif (yi - y0) <= y1 / x1 * (xi - x0) + tol and (yi - y0) >= (y1 - y0) - (xi - x0) - tol: # Q2 + ub[i] = math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2) + dub[i] = 1 / 2 * rx_ub / ub[i] + if not intrados_null: + lb[i] = math.sqrt((rx_lb) ** 2 - ((xi - x0_lb) - rx_lb) ** 2) + dlb[i] = -1 / 2 * rx_lb / lb[i] + elif (yi - y0) >= y1 / x1 * (xi - x0) - tol and (yi - y0) <= (y1 - y0) - (xi - x0) + tol: # Q4 + ub[i] = math.sqrt((rx_ub) ** 2 - ((xi - x0_ub) - rx_ub) ** 2) + dub[i] = 1 / 2 * rx_ub / ub[i] + if not intrados_null: + lb[i] = math.sqrt((rx_lb) ** 2 - ((xi - x0_lb) - rx_lb) ** 2) + dlb[i] = -1 / 2 * rx_lb / lb[i] + else: + print("Error Q. (x,y) = ({0},{1})".format(xi, yi)) + + return dub, dlb # ub, lb + + +def pavillionvault_bound_react(x, y, thk, fixed, x_span=(0.0, 10.0), y_span=(0.0, 10.0)): + """Computes the reaction bounds of a pavillion vault for a given thickness + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the vault + fixed : list + The list with indexes of the fixed vertices + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + + Returns + ------- + b : array + Values of the reaction bounds + """ + + x0, x1 = x_span + y0, y1 = y_span + b = zeros((len(fixed), 2)) + + for i in range(len(fixed)): + index = fixed[i] + xi, yi = x[index], y[index] + if xi == x0: + b[[i], :] += [-thk / 2, 0] + elif xi == x1: + b[i, :] += [+thk / 2, 0] + if yi == y0: + b[i, :] += [0, -thk / 2] + elif yi == y1: + b[i, :] += [0, +thk / 2] + + return abs(b) + + +def pavillionvault_bound_react_derivatives(x, y, thk, fixed, x_span=(0.0, 10.0), y_span=(0.0, 10.0)): + """Computes the sensitivities of the reaction bounds of a pavillion vault for a given thickness + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the vault + fixed : list + The list with indexes of the fixed vertices + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + + Returns + ------- + db : array + Values of the sensitivities of the reaction bounds in the points + """ + + x0, x1 = x_span + y0, y1 = y_span + db = zeros((len(fixed), 2)) + + for i in range(len(fixed)): + index = fixed[i] + xi, yi = x[index], y[index] + if xi == x0: + db[i, :] += [-1 / 2, 0] + elif xi == x1: + db[i, :] += [+1 / 2, 0] + if yi == y0: + db[i, :] += [0, -1 / 2] + elif yi == y1: + db[i, :] += [0, +1 / 2] + + return abs(db) + + +class PavillionVaultEnvelope(ParametricEnvelope): + def __init__( + self, + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, + spr_angle: float = 0.0, + **kwargs, + ): + super().__init__(thickness=thickness, **kwargs) + self.x_span = x_span + self.y_span = y_span + self.min_lb = min_lb + self.n = n + self.spr_angle = spr_angle + + self.update_envelope() # Generate the intra/extra/middle meshes + + @property + def __data__(self): + data = super().__data__ + data["x_span"] = self.x_span + data["y_span"] = self.y_span + data["min_lb"] = self.min_lb + data["n"] = self.n + data["spr_angle"] = self.spr_angle + return data + + def __str__(self): + return f"PavillionVaultEnvelope(name={self.name})" + + def update_envelope(self): + intrados, extrados, middle = pavillionvault_envelope( + x_span=self.x_span, y_span=self.y_span, thickness=self.thickness, min_lb=self.min_lb, n=self.n, spr_angle=self.spr_angle + ) + self.intrados = intrados + self.extrados = extrados + self.middle = middle + + def compute_middle(self, x, y): + return pavillionvault_middle(x, y, x_span=self.x_span, y_span=self.y_span, spr_angle=self.spr_angle, tol=1e-6) + + def compute_bounds(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pavillionvault_bounds(x, y, thickness, self.min_lb, x_span=self.x_span, y_span=self.y_span, spr_angle=self.spr_angle, tol=1e-6) + + def compute_bounds_derivatives(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pavillionvault_bounds_derivatives(x, y, thickness, self.min_lb, x_span=self.x_span, y_span=self.y_span, tol=1e-6) + + def compute_bound_react(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pavillionvault_bound_react(x, y, thickness, fixed, x_span=self.x_span, y_span=self.y_span, tol=1e-6) + + def compute_bound_react_derivatives(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pavillionvault_bound_react_derivatives(x, y, thickness, fixed, x_span=self.x_span, y_span=self.y_span, tol=1e-6) diff --git a/src/compas_tna/envelope/pointedvault.py b/src/compas_tna/envelope/pointedvault.py new file mode 100644 index 00000000..a8568413 --- /dev/null +++ b/src/compas_tna/envelope/pointedvault.py @@ -0,0 +1,660 @@ +import math + +from numpy import array +from numpy import ones +from numpy import zeros + +from compas.datastructures import Mesh +from compas_tna.diagrams.diagram_rectangular import create_cross_mesh +from compas_tna.envelope.parametricenvelope import ParametricEnvelope + + +def pointedvault_envelope( + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, + hc: float = 8.0, + he: list = None, + hm: list = None, +): + """Create an envelope for a pointed cross vault geometry with given parameters. + + Parameters + ---------- + cls : class + The Envelope class to use for creating the envelope. + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + thickness : float, optional + Thickness of the vault, by default 0.50 + min_lb : float, optional + Parameter for lower bound in nodes in the boundary, by default 0.0 + n : int, optional + Number of vertices for the mesh, by default 100 + hc : float, optional + Height in the middle point of the vault, by default 8.0 + he : list, optional + Height of the opening mid-span for each of the quadrants, by default None + hm : list, optional + Height of each quadrant center (spadrel), by default None + rho : float, optional + Density of the material in kN/m³, by default 25.0 + + Returns + ------- + envelope : Envelope + The created envelope with intrados, extrados, and middle meshes. + """ + # Create base topology + base_topology = create_cross_mesh(x_span=x_span, y_span=y_span, n=n) + xyz0, faces_i = base_topology.to_vertices_and_faces() + xi, yi, _ = array(xyz0).transpose() + + # Create middle surface + zt = pointedvault_middle(xi, yi, min_lb, x_span=x_span, y_span=y_span, hc=hc, he=he, hm=hm, tol=1e-6) + xyzt = array([xi, yi, zt.flatten()]).transpose() + middle = Mesh.from_vertices_and_faces(xyzt, faces_i) + middle.update_default_vertex_attributes(thickness=thickness) + + # Create upper and lower bounds + zub, zlb = pointedvault_bounds(xi, yi, thickness, min_lb, x_span=x_span, y_span=y_span, hc=hc, he=he, hm=hm, tol=1e-6) + xyzub = array([xi, yi, zub.flatten()]).transpose() + xyzlb = array([xi, yi, zlb.flatten()]).transpose() + + extrados = Mesh.from_vertices_and_faces(xyzub, faces_i) + intrados = Mesh.from_vertices_and_faces(xyzlb, faces_i) + + return intrados, extrados, middle + + +def pointedvault_middle( + x, + y, + min_lb, + x_span=(0.0, 10.0), + y_span=(0.0, 10.0), + hc=8.0, + he=None, + hm=None, + tol=1e-6, +): + """Compute middle of a pointed vault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + hc : float, optional + Height in the middle point of the vault, by default 8.0 + he : [float, float, float, float], optional + Height of the opening mid-span for each of the quadrants, by default None + hm : [float, float, float, float], optional + Height of each quadrant center (spadrel), by default None + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + middle : array + Values of the middle surface in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + lx = x1 - x0 + ly = y1 - y0 + + if he and hm is None: + h1, k1, r1 = _circle_3points_xy([x0, he[1]], [(x1 + x0) / 2, hc], [x1, he[0]]) + h2, k2, r2 = h1, k1, r1 + h3, k3, r3 = _circle_3points_xy([y0, he[3]], [(y1 + y0) / 2, hc], [y1, he[2]]) + h4, k4, r4 = h3, k3, r3 + elif hm and he: + h1, k1, r1 = _circle_3points_xy([(x1 + x0) / 2, hc], [3 * (x1 + x0) / 4, hm[0]], [x1, he[0]]) + h2, k2, r2 = _circle_3points_xy([(x1 + x0) / 2, hc], [1 * (x1 + x0) / 4, hm[1]], [x0, he[1]]) + h3, k3, r3 = _circle_3points_xy([(y1 + y0) / 2, hc], [3 * (y1 + y0) / 4, hm[2]], [y1, he[2]]) + h4, k4, r4 = _circle_3points_xy([(y1 + y0) / 2, hc], [1 * (y1 + y0) / 4, hm[3]], [y0, he[3]]) + + middle = zeros((len(x), 1)) + + for i in range(len(x)): + xi, yi = x[i], y[i] + + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + # Equation (xi - hx) ** 2 + (hi - kx) ** 2 = rx **2 to find the height of the pointed part (middle of quadrant) with that height one find the equivalent radius + if he: + hi = k1 + math.sqrt(r1**2 - (xi - h1) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, ly) # This in the equation ri ** 2 = (xi - xc_) ** 2 + (zi - zc_) ** 2 -> zc = 0.0 and xc_ = (x0 + x1)/2 + if yi <= (y1 + y0) / 2: + zi = _sqrt((ri) ** 2 - (yi - (y0 + ri)) ** 2) + else: + zi = _sqrt((ri) ** 2 - (yi - (y1 - ri)) ** 2) + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + # Equation (xi - hy) ** 2 + (hi - ky) ** 2 = ry **2 to find the height of the pointed part (middle of quadrant) with that height one find the equivalent radius + if he: + hi = k3 + math.sqrt(r3**2 - (yi - h3) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) # This in the equation ri ** 2 = (xi - xc_) ** 2 + (zi - zc_) ** 2 -> zc = 0.0 and xc_ = (x0 + x1)/2 + if xi <= (x0 + x1) / 2: + zi = _sqrt((ri) ** 2 - (xi - (x0 + ri)) ** 2) + else: + zi = _sqrt((ri) ** 2 - (xi - (x1 - ri)) ** 2) + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + if he: + hi = k2 + math.sqrt(r2**2 - (xi - h2) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, ly) + if yi <= (y1 + y0) / 2: + zi = _sqrt((ri) ** 2 - (yi - (y0 + ri)) ** 2) + else: + zi = _sqrt((ri) ** 2 - (yi - (y1 - ri)) ** 2) + + elif yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + if he: + hi = k4 + math.sqrt(r4**2 - (yi - h4) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) + if xi <= (x0 + x1) / 2: + zi = _sqrt((ri) ** 2 - (xi - (x0 + ri)) ** 2) + else: + zi = _sqrt((ri) ** 2 - (xi - (x1 - ri)) ** 2) + + else: + print("Vertex did not belong to any Q. (x,y) = ({0},{1})".format(xi, yi)) + + middle[i] = zi + + return middle + + +def pointedvault_bounds( + x, + y, + thk, + min_lb, + x_span=(0.0, 10.0), + y_span=(0.0, 10.0), + hc=8.0, + he=None, + hm=None, + tol=1e-6, +): + """Compute upper and lower bounds of a pointed vault based on the parameters. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the arch + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + hc : float, optional + Height in the middle point of the vault, by default 8.0 + he : [float, float, float, float], optional + Height of the opening mid-span for each of the quadrants, by default None + hm : [float, float, float, float], optional + Height of each quadrant center (spadrel), by default None + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + ub : array + Values of the upper bound in the points + lb : array + Values of the lower bound in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + lx = x1 - x0 + ly = y1 - y0 + + if he: + he_ub = he.copy() + he_lb = he.copy() + for i in range(len(he)): + he_ub[i] += thk / 2 + he_lb[i] -= thk / 2 + if hm: + raise NotImplementedError() + + if he and hm is None: + h1, k1, r1 = _circle_3points_xy([x0, he[1]], [(x1 + x0) / 2, hc], [x1, he[0]]) + h2, k2, r2 = h1, k1, r1 + h3, k3, r3 = _circle_3points_xy([y0, he[3]], [(y1 + y0) / 2, hc], [y1, he[2]]) + h4, k4, r4 = h3, k3, r3 + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + + for i in range(len(x)): + xi, yi = x[i], y[i] + + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + if he: + hi = k1 + math.sqrt(r1**2 - (xi - h1) ** 2) + else: + hi = hc + + ri = _find_r_given_h_l(hi, ly) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if yi <= (y1 + y0) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y1 - ri)) ** 2) + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + if he: + hi = k3 + math.sqrt(r3**2 - (yi - h3) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if xi <= (x0 + x1) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x1 - ri)) ** 2) + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + if he: + hi = k2 + math.sqrt(r2**2 - (xi - h2) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, ly) + ri_lb = ri - thk / 2 + ri_ub = ri + thk / 2 + if yi <= (y1 + y0) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y1 - ri)) ** 2) + + elif yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + if he: + hi = k4 + math.sqrt(r4**2 - (yi - h4) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if xi <= (x0 + x1) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x1 - ri)) ** 2) + else: + print("Vertex did not belong to any Q. (x,y) = ({0},{1})".format(xi, yi)) + + if ((yi) > y1_lb and ((xi) > x1_lb or (xi) < x0_lb)) or ((yi) < y0_lb and ((xi) > x1_lb or (xi) < x0_lb)): + lb[i] = -1 * min_lb + + return ub, lb + + +def pointedvault_bounds_derivatives( + x, + y, + thk, + min_lb, + x_span=(0.0, 10.0), + y_span=(0.0, 10.0), + hc=8.0, + he=None, + hm=None, + tol=1e-6, +): + """Computes the sensitivities of upper and lower bounds in the x, y coordinates and thickness specified. + + Parameters + ---------- + x : list + x-coordinates of the points + y : list + y-coordinates of the points + thk : float + Thickness of the arch + min_lb : float + Parameter for lower bound in nodes in the boundary + x_span : tuple, optional + Span of the vault in x direction, by default (0.0, 10.0) + y_span : tuple, optional + Span of the vault in y direction, by default (0.0, 10.0) + hc : float, optional + Height in the middle point of the vault, by default 8.0 + he : [float, float, float, float], optional + Height of the opening mid-span for each of the quadrants, by default None + hm : [float, float, float, float], optional + Height of each quadrant center (spadrel), by default None + tol : float, optional + Tolerance, by default 1e-6 + + Returns + ------- + dub : array + Values of the sensitivities for the upper bound in the points + dlb : array + Values of the sensitivities for the lower bound in the points + """ + + y1 = y_span[1] + y0 = y_span[0] + x1 = x_span[1] + x0 = x_span[0] + + y1_lb = y1 - thk / 2 + y0_lb = y0 + thk / 2 + x1_lb = x1 - thk / 2 + x0_lb = x0 + thk / 2 + + lx = x1 - x0 + ly = y1 - y0 + + if he: + he_ub = he.copy() + he_lb = he.copy() + for i in range(len(he)): + he_ub[i] += thk / 2 + he_lb[i] -= thk / 2 + if hm: + raise NotImplementedError() + hm_ub = hm.copy() + hm_lb = hm.copy() + for i in range(len(hm)): + hm_ub[i] += thk / 2 + hm_lb[i] -= thk / 2 + + if he and hm is None: + h1, k1, r1 = _circle_3points_xy([x0, he[1]], [(x1 + x0) / 2, hc], [x1, he[0]]) + h2, k2, r2 = h1, k1, r1 + h3, k3, r3 = _circle_3points_xy([y0, he[3]], [(y1 + y0) / 2, hc], [y1, he[2]]) + h4, k4, r4 = h3, k3, r3 + + ub = ones((len(x), 1)) + lb = ones((len(x), 1)) * -min_lb + dub = zeros((len(x), 1)) + dlb = zeros((len(x), 1)) + + for i in range(len(x)): + xi, yi = x[i], y[i] + + if yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q1 + if he: + hi = k1 + math.sqrt(r1**2 - (xi - h1) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, ly) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if yi <= (y1 + y0) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y1 - ri)) ** 2) + dub[i] = 1 / 2 * ri_ub / ub[i] + dlb[i] = -1 / 2 * ri_lb / lb[i] + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi >= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) - tol: # Q3 + if he: + hi = k3 + math.sqrt(r3**2 - (yi - h3) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if xi <= (x0 + x1) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x1 - ri)) ** 2) + dub[i] = 1 / 2 * ri_ub / ub[i] + dlb[i] = -1 / 2 * ri_lb / lb[i] + + elif yi >= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) - tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q2 + if he: + hi = k2 + math.sqrt(r2**2 - (xi - h2) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, ly) + ri_lb = ri - thk / 2 + ri_ub = ri + thk / 2 + if yi <= (y1 + y0) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (yi - (y1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (yi - (y1 - ri)) ** 2) + dub[i] = 1 / 2 * ri_ub / ub[i] + dlb[i] = -1 / 2 * ri_lb / lb[i] + + elif yi <= y0 + (y1 - y0) / (x1 - x0) * (xi - x0) + tol and yi <= y1 - (y1 - y0) / (x1 - x0) * (xi - x0) + tol: # Q4 + if he: + hi = k4 + math.sqrt(r4**2 - (yi - h4) ** 2) + else: + hi = hc + ri = _find_r_given_h_l(hi, lx) + ri_ub = ri + thk / 2 + ri_lb = ri - thk / 2 + if xi <= (x0 + x1) / 2: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x0 + ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x0 + ri)) ** 2) + else: + ub[i] = _sqrt((ri_ub) ** 2 - (xi - (x1 - ri)) ** 2) + lb[i] = _sqrt((ri_lb) ** 2 - (xi - (x1 - ri)) ** 2) + dub[i] = 1 / 2 * ri_ub / ub[i] + dlb[i] = -1 / 2 * ri_lb / lb[i] + else: + print("Vertex did not belong to any Q. (x,y) = ({0},{1})".format(xi, yi)) + + if ((yi) > y1_lb and ((xi) > x1_lb or (xi) < x0_lb)) or ((yi) < y0_lb and ((xi) > x1_lb or (xi) < x0_lb)): + lb[i] = -1 * min_lb + dlb[i] = 0.0 + + return dub, dlb # ub, lb + + +def pointedvault_bound_react( + x, + y, + thk, + min_lb, + x_span=(0.0, 10.0), + y_span=(0.0, 10.0), + hc=8.0, + he=None, + hm=None, + tol=1e-6, +): + """Compute the bounds on the reaction vector of the pointed cross vault.""" + pass + + +def pointedvault_bound_react_derivatives( + x, + y, + thk, + min_lb, + x_span=(0.0, 10.0), + y_span=(0.0, 10.0), + hc=8.0, + he=None, + hm=None, + tol=1e-6, +): + """Compute the sensitivities of the bounds on the reaction vector of the pointed cross vault.""" + pass + + +def _find_r_given_h_l(h, length): + r = h**2 / length + length / 4 + + return r + + +def _circle_3points_xy(p1, p2, p3): + x1 = p1[0] + z1 = p1[1] + x2 = p2[0] + z2 = p2[1] + x3 = p3[0] + z3 = p3[1] + + x12 = x1 - x2 + x13 = x1 - x3 + z12 = z1 - z2 + z13 = z1 - z3 + z31 = z3 - z1 + z21 = z2 - z1 + x31 = x3 - x1 + x21 = x2 - x1 + + sx13 = x1**2 - x3**2 + sz13 = z1**2 - z3**2 + sx21 = x2**2 - x1**2 + sz21 = z2**2 - z1**2 + + f = ((sx13) * (x12) + (sz13) * (x12) + (sx21) * (x13) + (sz21) * (x13)) / (2 * ((z31) * (x12) - (z21) * (x13))) + g = ((sx13) * (z12) + (sz13) * (z12) + (sx21) * (z13) + (sz21) * (z13)) / (2 * ((x31) * (z12) - (x21) * (z13))) + c = -(x1**2) - z1**2 - 2 * g * x1 - 2 * f * z1 + h = -g + k = -f + r2 = h * h + k * k - c + r = math.sqrt(r2) + + return h, k, r + + +def _sqrt(x): + try: + sqrt_x = math.sqrt(x) + except BaseException: + if x > -10e4: + sqrt_x = math.sqrt(abs(x)) + else: + sqrt_x = 0.0 + return sqrt_x + + +class PointedVaultEnvelope(ParametricEnvelope): + def __init__( + self, + x_span: tuple = (0.0, 10.0), + y_span: tuple = (0.0, 10.0), + thickness: float = 0.50, + min_lb: float = 0.0, + n: int = 100, + hc: float = 5.0, + he: list = None, + hm: list = None, + **kwargs, + ): + super().__init__(thickness=thickness, **kwargs) + self.x_span = x_span + self.y_span = y_span + self.min_lb = min_lb + self.n = n + self.hc = hc + self.he = he + self.hm = hm + + self.update_envelope() # Generate the intra/extra/middle meshes + + @property + def __data__(self): + data = super().__data__ + data["x_span"] = self.x_span + data["y_span"] = self.y_span + data["min_lb"] = self.min_lb + data["n"] = self.n + data["hc"] = self.hc + data["he"] = self.he + data["hm"] = self.hm + return data + + def __str__(self): + return f"PointedVaultEnvelope(name={self.name})" + + def update_envelope(self): + intrados, extrados, middle = pointedvault_envelope( + x_span=self.x_span, y_span=self.y_span, thickness=self.thickness, min_lb=self.min_lb, n=self.n, hc=self.hc, he=self.he, hm=self.hm + ) + self.intrados = intrados + self.extrados = extrados + self.middle = middle + + def compute_middle(self, x, y): + return pointedvault_middle(x, y, self.min_lb, self.x_span, self.y_span, self.hc, self.he, self.hm) + + def compute_bounds(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pointedvault_bounds(x, y, thickness, self.min_lb, self.x_span, self.y_span, self.hc, self.he, self.hm) + + def compute_bounds_derivatives(self, x, y, thickness=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pointedvault_bounds_derivatives(x, y, thickness, self.min_lb, self.x_span, self.y_span, self.hc, self.he, self.hm) + + def compute_bound_react(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pointedvault_bound_react(x, y, thickness, self.min_lb, self.x_span, self.y_span, self.hc, self.he, self.hm) + + def compute_bound_react_derivatives(self, x, y, thickness=None, fixed=None): + if thickness is None: + thickness = self.thickness + else: + self.thickness = thickness + return pointedvault_bound_react_derivatives(x, y, thickness, self.min_lb, self.x_span, self.y_span, self.hc, self.he, self.hm)