Merge pull request #43 from ricardoavelino/fill

Add Fill Mesh commands

Author: tomvanmele

CHANGELOG.md

@@ -35,6 +35,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added comprehensive mesh creation functions in `diagram_rectangular.py`:
 * Added circular mesh creation functions in `diagram_circular.py` and corrected oculus and diagonal properties.
 * Added arch mesh creation functions in `diagram_arch.py`:
+* Added option to add `fill` to envelope. TODO: properly deal with 'pz' summation in future
 
 ### Changed
 

src/compas_tna/envelope/dome.py

@@ -225,18 +225,17 @@ def dome_bound_react(x, y, thk, fixed, center=(5.0, 5.0), radius=5.0):
 
     Returns
     -------
-    b : array
+    b : array (nb x 2)
         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))
+        x_ = thk / 2 * math.cos(theta)
+        y_ = thk / 2 * math.sin(theta)
         b[i, 0] = x_
         b[i, 1] = y_
 

src/compas_tna/envelope/envelope.py

@@ -10,10 +10,12 @@
 class Envelope(Data):
     """Pure geometric envelope representing masonry structure boundaries."""
 
-    def __init__(self, rho: Optional[float] = 20.0, is_parametric: bool = False, **kwargs):
+    def __init__(self, rho: Optional[float] = 20.0, rho_fill: Optional[float] = 14.0, is_parametric: bool = False, **kwargs):
         super().__init__(**kwargs)
 
         self.rho = rho
+        self.fill = None
+        self.rho_fill = rho_fill
         self._is_parametric = is_parametric
 
         # Computed properties (cached)
@@ -25,6 +27,7 @@ def __init__(self, rho: Optional[float] = 20.0, is_parametric: bool = False, **k
     def __data__(self):
         data = {}
         data["rho"] = self.rho
+        data["rho_fill"] = self.rho_fill
         data["is_parametric"] = self.is_parametric
         return data
 
@@ -89,6 +92,11 @@ def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram) -> None:
 
         raise NotImplementedError("Implement apply_selfweight_to_formdiagram for specific envelope type.")
 
+    def apply_fill_weight_to_formdiagram(self, formdiagram: FormDiagram) -> None:
+        """Apply selfweight to the nodes of a form diagram based on the appropriate method."""
+
+        raise NotImplementedError("Implement apply_fill_weight_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."""
 

src/compas_tna/envelope/meshenvelope.py

@@ -9,6 +9,29 @@
 from compas_tna.envelope import Envelope
 
 
+def griddata_project(xy: list[list[float]], xyz_target: list[list[float]], method="linear"):
+    """Project a point cloud onto a target point cloud using griddata interpolation.
+
+    Parameters
+    ----------
+    xy : list[list[float]]
+        The XY coordinates of the points to project.
+    xyz_target : list[list[float]]
+        The XYZ coordinates of the target points.
+    method : str, optional
+        The method to use for interpolation. Default is "linear".
+
+    Returns
+    -------
+    list[float]
+        The projected Z coordinates.
+    """
+    xy = asarray(xy)
+    xy_target = asarray(xyz_target)[:, :2]
+    z_target = asarray(xyz_target)[:, 2]
+    return griddata(xy_target, z_target, xy, method=method).tolist()
+
+
 def interpolate_middle_mesh(intrados: Mesh, extrados: Mesh) -> Mesh:
     """Interpolate a middle mesh between intrados and extrados meshes.
 
@@ -30,25 +53,18 @@ def interpolate_middle_mesh(intrados: Mesh, extrados: Mesh) -> Mesh:
     # 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 Z coordinates from both surfaces based on the same XY coordinates
+    zi = middle.vertices_attribute("z")
+    ze = griddata_project(middle.vertices_attributes("xy"), extrados.vertices_attributes("xyz"), method="linear")
 
-    # 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
+    # First loop: set middle Z as average of intrados and extrados
     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)
+        _, _, nz = middle.vertex_normal(key)
         z_diff = ze[i] - zi[i]
         if abs(nz) > 0.1:
             thickness = abs(z_diff) * abs(nz)
@@ -158,7 +174,6 @@ def project_mesh_to_target_vertica_nearest(mesh: Mesh, target: Mesh) -> None:
         mesh.vertex_attributes(vertex, "xyz", new_point)
 
 
-# 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.
 
@@ -174,14 +189,7 @@ def project_mesh_to_target_vertical(mesh: Mesh, target: Mesh) -> None:
     None
         The mesh is modified in place.
     """
-    # Get point clouds for interpolation
-    target_points = asarray(target.vertices_attributes("xyz"))
-
-    # Get XY coordinates of middle mesh
-    mesh_xy = asarray(mesh.vertices_attributes("xy"))
-
-    # Interpolate Z coordinates from both surfaces
-    z_target = griddata(target_points[:, :2], target_points[:, 2], mesh_xy, method="linear").tolist()
+    z_target = griddata_project(mesh.vertices_attributes("xy"), target.vertices_attributes("xyz"), method="linear")
 
     for key, i in enumerate(mesh.vertices()):
         mesh.vertex_attribute(key, "z", z_target[i])
@@ -477,7 +485,7 @@ def compute_area(self) -> float:
         return self.middle.area()
 
     # =============================================================================
-    # TNA-specific operations (accept formdiagram as parameter)
+    # Loads operations
     # =============================================================================
 
     def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram, normalize=True) -> None:
@@ -539,6 +547,51 @@ def apply_selfweight_to_formdiagram(self, formdiagram: FormDiagram, normalize=Tr
 
         print(f"Selfweight applied to form diagram. Total load: {sum(abs(formdiagram.vertex_attribute(vertex, 'pz')) for vertex in formdiagram.vertices())}")
 
+    def apply_fill_weight_to_formdiagram(self, formdiagram: FormDiagram) -> None:
+        """Apply fill weight to the nodes of a form diagram based on the fill surface and local thicknesses."""
+        if self.fill is None or self.extrados is None:
+            raise ValueError("Fill mesh is not set. Please set the fill mesh and extrados before applying fill weight.")
+
+        # Step 2: Copy the form diagram for projection
+        form_fill = formdiagram.copy()  # For upper bound (extrados)
+        form_ub = formdiagram.copy()  # For lower bound (intrados)
+        form_zero = formdiagram.copy()  # For zero bound (extrados)
+        form_zero.vertices_attribute("z", 0.0)
+
+        # Step 3: Project form diagram onto extrados (upper bound)
+        project_mesh_to_target_vertical(form_fill, self.fill)
+        project_mesh_to_target_vertical(form_ub, self.extrados)
+
+        fill_weight = 0.0
+
+        # Step 4: Collect heights and assign to form diagram
+        for vertex in formdiagram.vertices():
+            if vertex in form_ub.vertices() and vertex in form_fill.vertices():
+                # Get z coordinates from projected meshes
+                _, _, z_fill = form_fill.vertex_coordinates(vertex)
+                _, _, z_ub = form_ub.vertex_coordinates(vertex)
+                a0 = form_zero.vertex_area(vertex)
+
+                if z_fill < z_ub:
+                    z_fill = z_ub
+
+                zdiff = z_fill - z_ub
+                pz_fill = -a0 * zdiff * self.rho_fill
+                pz0 = formdiagram.vertex_attribute(vertex, "pz")
+                formdiagram.vertex_attribute(vertex, "pz", pz_fill + pz0)
+                fill_weight += abs(pz_fill)
+
+                # Store z_fill
+                formdiagram.vertex_attribute(vertex, "zfill", z_fill)
+            else:
+                print(f"Warning: Vertex {vertex} not found in projected meshes")
+
+        print(f"Fill weight applied to form diagram. Total load: {fill_weight}")
+
+    # =============================================================================
+    # Envelope and target projection operations
+    # =============================================================================
+
     def apply_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None:
         """Apply envelope bounds to a form diagram based on the intrados and extrados surfaces.
 

src/compas_tna/envelope/parametricenvelope.py

@@ -208,7 +208,7 @@ def apply_reaction_bounds_to_formdiagram(self, formdiagram: FormDiagram) -> None
         None
             The FormDiagram is modified in place.
         """
-        fixed: list[int] = formdiagram.vertices_where({"is_support": True})
+        fixed: list[int] = list(formdiagram.vertices_where({"is_support": True}))
         xy = np.array(formdiagram.vertices_attributes("xy"))
         bound_react = self.compute_bound_react(xy[:, 0], xy[:, 1], self.thickness, fixed)
         for i, key in enumerate(fixed):