Feature/radial mesh controls (#21)

* Started on adding radial mesh controls

* Ran black

* Progress on mesh constraints

* Changed element type to quad

* Added element type enum

* Quad updates

* Updated github action ubuntu version:

Co-authored-by: Ben Cannell <[email protected]>

Author: thisistheplace

.github/workflows/test.yml

@@ -4,7 +4,7 @@ on: [push]
 
 jobs:
   test:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-22.04
     if: github.ref == 'refs/heads/main'
 
     strategy:

.vscode/settings.json

@@ -1,6 +1,6 @@
 {
     "python.testing.pytestArgs": [
-        "."
+        "src"
     ],
     "python.testing.unittestEnabled": false,
     "python.testing.pytestEnabled": true

Dockerfile-debug

@@ -15,6 +15,6 @@ RUN pip install pytest
 
 # Install Python dependencies
 ADD requirements.txt .
-RUN pip install --no-cache-dir -r requirements.txt && pip install --no-cache-dir gunicorn
+RUN pip install --no-cache-dir -r requirements.txt
 
 WORKDIR /workspaces/joint_model

docker-compose.yml

@@ -4,7 +4,7 @@ services:
   fastapi:
     container_name: jointrest
     environment:
-      - RESTAPI_URL=http://127.0.0.1:8000
+      - RESTAPI_URL=http://127.0.0.1:8020
       - VIEWER_URL=http://127.0.0.1:8050
     restart: always
     build:
@@ -13,15 +13,15 @@ services:
     expose:
       - "8000"
     ports:
-      - "8000:8000"
+      - "8020:8000"
     command: uvicorn app.app_rest:app --host 0.0.0.0 --port 8000 --proxy-headers --workers 1
 
   viewer:
     depends_on:
       - fastapi
     container_name: jointviewer
     environment:
-      - RESTAPI_URL=http://host.docker.internal:8000
+      - RESTAPI_URL=http://host.docker.internal:8020
       - VIEWER_URL=http://host.docker.internal:8050
     restart: always
     build:

src/app/converters/mesh.py

@@ -1,34 +1,35 @@
 from collections import defaultdict
+from enum import Enum
 import gmsh
 
 from ..interfaces import DashVtkMesh
 
+class ElementType(str, Enum):
+    TRIANGLE = "Triangle"
+    QUADRANGLE = "Quadrangle"
 
-def process_mesh(mesh: gmsh.model.mesh):
-    # Like elements, mesh edges and faces are described by (an ordered list of)
-    # their nodes. Let us retrieve the edges and the (triangular) faces of all the
-    # first order tetrahedra in the mesh:
-    elementType = mesh.getElementType("Triangle", 1)
-    # edgeNodes = mesh.getElementEdgeNodes(elementType)
-    faceNodes = mesh.getElementFaceNodes(elementType, 3)
+NUM_NODES = {
+    ElementType.TRIANGLE: 3,
+    ElementType.QUADRANGLE: 4
+}
 
-    # Edges and faces are returned for each element as a list of nodes corresponding
-    # to the canonical orientation of the edges and faces for a given element type.
-    # Edge and face tags can then be retrieved by providing their nodes:
-    # edgeTags, edgeOrientations = mesh.getEdges(edgeNodes)
-    faceTags, _ = mesh.getFaces(3, faceNodes)
+def process_mesh(mesh: gmsh.model.mesh, eltype: ElementType):
+    elementType = mesh.getElementType(eltype.value, 1)
+    faceNodes = mesh.getElementFaceNodes(elementType, NUM_NODES[eltype])
+
+    faceTags, _ = mesh.getFaces(NUM_NODES[eltype], faceNodes)
     elementTags, _ = mesh.getElementsByType(elementType)
 
     faces2Elements = defaultdict(list)
 
-    for i in range(len(faceTags)):  # 4 faces per tetrahedron
-        faces2Elements[faceTags[i]].append(elementTags[i // 4])
+    for i in range(len(faceTags)):
+        faces2Elements[faceTags[i]].append(elementTags[i])
 
     return faces2Elements
 
 
-def mesh_to_dash_vtk(mesh: gmsh.model.mesh) -> DashVtkMesh:
-    face2el = process_mesh(mesh)
+def mesh_to_dash_vtk(mesh: gmsh.model.mesh, eltype: ElementType) -> DashVtkMesh:
+    face2el = process_mesh(mesh, eltype)
 
     outerfaces = [k for k in face2el.keys()]
 
@@ -41,7 +42,7 @@ def mesh_to_dash_vtk(mesh: gmsh.model.mesh) -> DashVtkMesh:
         _, node_ids = mesh.getElement(element)
         # if eltype != 2 or len(node_ids) != 3:
         #     continue
-        line = [len(node_ids)]
+        line = [len(node_ids) + 1]
         poly = [len(node_ids)]
         for nid in node_ids:
             if nid not in nid2pointidx:
@@ -50,6 +51,7 @@ def mesh_to_dash_vtk(mesh: gmsh.model.mesh) -> DashVtkMesh:
                 nid2pointidx[nid] = int(len(points) / 3 - 1)
             line.append(nid2pointidx[nid])
             poly.append(nid2pointidx[nid])
+        line.append(line[1])
 
         lines += line
         polys += poly

src/app/converters/model.py

@@ -1,10 +1,10 @@
-from .mesh import mesh_to_dash_vtk
+from .mesh import mesh_to_dash_vtk, ElementType
 
 from ..interfaces import *
 from ..modelling.mesher.mesh import mesh_model
 
 
 def convert_model_to_dash_vtk(model: Model) -> DashVtkModel:
     with mesh_model(model, MeshSpecs(size=0.1)) as mesh:
-        mesh = mesh_to_dash_vtk(mesh)
+        mesh = mesh_to_dash_vtk(mesh, ElementType.QUADRANGLE)
         return DashVtkModel(name=model.name, mesh=mesh)

src/app/interfaces/examples/joints.py

@@ -7,17 +7,18 @@
         joint=Joint(
             name="TJoint",
             master=Tubular(
-                    name="Vertical",
-                    axis=Axis3D(
-                        point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
-                    ),
-                    diameter=0.5,
+                name="Vertical",
+                axis=Axis3D(
+                    point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
                 ),
+                diameter=0.5,
+            ),
             slaves=[
                 Tubular(
                     name="Horizontal",
                     axis=Axis3D(
-                        point=Point3D(x=0., y=0.25, z=0), vector=Vector3D(x=0, y=3, z=0)
+                        point=Point3D(x=0.0, y=0.25, z=0),
+                        vector=Vector3D(x=0, y=3, z=0),
                     ),
                     diameter=0.25,
                 ),
@@ -29,17 +30,18 @@
         joint=Joint(
             name="TAngle",
             master=Tubular(
-                    name="Vertical",
-                    axis=Axis3D(
-                        point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
-                    ),
-                    diameter=0.5,
+                name="Vertical",
+                axis=Axis3D(
+                    point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
                 ),
+                diameter=0.5,
+            ),
             slaves=[
                 Tubular(
                     name="Horizontal",
                     axis=Axis3D(
-                        point=Point3D(x=0., y=0.25, z=0), vector=Vector3D(x=0, y=3, z=3)
+                        point=Point3D(x=0.0, y=0.25, z=0),
+                        vector=Vector3D(x=0, y=3, z=3),
                     ),
                     diameter=0.25,
                 ),
@@ -51,12 +53,12 @@
         joint=Joint(
             name="TOffset",
             master=Tubular(
-                    name="Vertical",
-                    axis=Axis3D(
-                        point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
-                    ),
-                    diameter=0.5,
+                name="Vertical",
+                axis=Axis3D(
+                    point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
                 ),
+                diameter=0.5,
+            ),
             slaves=[
                 Tubular(
                     name="Horizontal",
@@ -73,12 +75,12 @@
         joint=Joint(
             name="KJoint",
             master=Tubular(
-                    name="Chord",
-                    axis=Axis3D(
-                        point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
-                    ),
-                    diameter=0.6,
+                name="Chord",
+                axis=Axis3D(
+                    point=Point3D(x=0, y=0, z=-5), vector=Vector3D(x=0, y=0, z=10)
                 ),
+                diameter=0.6,
+            ),
             slaves=[
                 Tubular(
                     name="Brace1",
@@ -90,21 +92,24 @@
                 Tubular(
                     name="Brace2",
                     axis=Axis3D(
-                        point=Point3D(x=-0.3, y=0.3, z=1), vector=Vector3D(x=-2, y=2, z=2)
+                        point=Point3D(x=-0.3, y=0.3, z=1),
+                        vector=Vector3D(x=-2, y=2, z=2),
                     ),
                     diameter=0.25,
                 ),
                 Tubular(
                     name="Brace3",
                     axis=Axis3D(
-                        point=Point3D(x=0.3, y=0.3, z=-1), vector=Vector3D(x=2, y=2, z=-2)
+                        point=Point3D(x=0.3, y=0.3, z=-1),
+                        vector=Vector3D(x=2, y=2, z=-2),
                     ),
                     diameter=0.4,
                 ),
                 Tubular(
                     name="Brace4",
                     axis=Axis3D(
-                        point=Point3D(x=-0.3, y=0.3, z=-1), vector=Vector3D(x=-2, y=2, z=-2)
+                        point=Point3D(x=-0.3, y=0.3, z=-1),
+                        vector=Vector3D(x=-2, y=2, z=-2),
                     ),
                     diameter=0.4,
                 ),

src/app/modelling/geometry/intersections.py

@@ -75,6 +75,7 @@ def intersection(master: NpTubular, slave: NpTubular) -> np.ndarray:
 
     return intersect3D
 
+
 def circle_intersect(
     center: np.ndarray, diameter: float, point: np.ndarray, vector: np.ndarray
 ) -> np.ndarray:
@@ -143,31 +144,37 @@ def flat_tube_intersection(master: NpTubular, slave: NpTubular) -> np.ndarray:
     # Get intersections on master surface
     point = intersection(master, slave)
     master_circle = sympy.Circle(master.axis.point.array[:2], master.diameter / 2.0)
-    angle = arc_angle_signed(master_circle, point) * -1.0 # -1 since rotation is X to Y
+    angle = arc_angle_signed(master_circle, point) * -1.0  # -1 since rotation is X to Y
     if angle > math.pi:
         angle = 2 * math.pi - angle
     elif angle < -1 * math.pi:
         angle = -2 * math.pi - angle
-    circ_length = (master_circle.circumference / 2.) * angle / math.pi
+    circ_length = (master_circle.circumference / 2.0) * angle / math.pi
     point[0] = circ_length
+    # Y coordinate is at tube radius
+    point[1] = master.diameter / 2.0
     return point
 
+
 def arc_angle_signed(circle: sympy.Circle, point: np.ndarray) -> float:
     """Seam (0 rads) is aligned with Y axis
-    
+
     Positive rotation from X to Y axis (clockwise!)
 
     Circle is in 2D X/Y plane
     """
-    seam = sympy.Point2D(0., circle.radius)
+    seam = sympy.Point2D(0.0, circle.radius)
     seg_vector = np.empty((3,))
     seg_vector[:2] = point[:2] - np.array(seam.coordinates)
     seg_vector[2] = 0.0
     seg_length = np.linalg.norm(seg_vector)
     sub_angle = math.acos(seg_length / 2.0 / circle.radius)
     return np.sign(circle.center[0] - point[0]) * (math.pi - 2 * sub_angle)
 
-def plane_intersect(axis: np.ndarray, point: np.ndarray, plane: sympy.Plane) -> np.ndarray:
+
+def plane_intersect(
+    axis: np.ndarray, point: np.ndarray, plane: sympy.Plane
+) -> np.ndarray:
     """Calculate 3D point where slave intersects plane
 
     Plane is in X/Z plane.
@@ -185,15 +192,17 @@ def plane_intersect(axis: np.ndarray, point: np.ndarray, plane: sympy.Plane) ->
     plane_point = np.array(plane.p1, dtype=float)
     plane_normal = np.array(plane.normal_vector, dtype=float)
     epsilon = 1e-6
-    
+
     try:
         ndotu = plane_normal.dot(axis)
         if abs(ndotu) < epsilon:
-            raise IntersectionError(f"Cannot compute intersect of vector {axis} which lies in plane {plane}")
+            raise IntersectionError(
+                f"Cannot compute intersect of vector {axis} which lies in plane {plane}"
+            )
         w = point - plane_point
         si = -plane_normal.dot(w) / ndotu
         return w + si * axis + plane_point
     except Exception as e:
         # Handle some specific exception here where an intersection is not found
         msg = f"Could not find intersection point of vector {axis} with plane {plane}.\nEncountered error:\n{e}"
-        raise IntersectionError(msg)
+        raise IntersectionError(msg)

src/app/modelling/geometry/vectors.py

@@ -30,4 +30,4 @@ def rotate(point: np.ndarray, axis: np.ndarray, angle: float) -> np.ndarray:
     rotation = transforms.rotation_matrix(angle, axis, point)
     rpoint = np.zeros((4,))
     rpoint[:3] = point[:3]
-    return np.dot(rotation, rpoint)[:3]
+    return np.dot(rotation, rpoint)[:3]

src/app/modelling/geometry/weld.py

@@ -4,20 +4,26 @@
 from typing import Generator
 
 from .vectors import rotate, unit_vector
-from .intersections import flat_tube_intersection, plane_intersect, intersection, arc_angle_signed
+from .intersections import (
+    flat_tube_intersection,
+    plane_intersect,
+    intersection,
+    arc_angle_signed,
+)
 from ...interfaces import *
 
 # TODO: look into https://mathcurve.com/courbes2d.gb/alain/alain.shtml
 # TODO:https://math.stackexchange.com/questions/141593/formula-for-cylinder?newreg=52c6b3d9cb1d47c1aaa3ff0706aa2298
 # Plane is XZ so X = Z, Y = offset
 
+
 def get_weld_intersect_points(
     master: NpTubular, slave: NpTubular, angle_inc: int = 10
 ) -> Generator[np.ndarray, np.ndarray, None]:
     """X/Z plane"""
     intersect = intersection(master, slave)
     # Radius point is on Y axis at radius
-    radius_point = np.array([0., master.diameter / 2., intersect[2]])
+    radius_point = np.array([0.0, master.diameter / 2.0, intersect[2]])
     # Create X/Z plane at radius point
     p2 = deepcopy(radius_point)
     p2[0] += 1.0
@@ -31,24 +37,31 @@ def get_weld_intersect_points(
     radius_arc_angle = arc_angle_signed(circle, radius_point)
     intersect_arc_angle = arc_angle_signed(circle, intersect)
     arc_angle = radius_arc_angle + intersect_arc_angle
-    slave_vector = rotate(unit_vector(slave.axis.vector.array), np.array([0.0, 0.0, 1.0]), arc_angle * -1.)
-    perp = unit_vector(np.cross(master.axis.vector.array, slave_vector)) * slave.diameter / 2.
+    slave_vector = rotate(
+        unit_vector(slave.axis.vector.array),
+        np.array([0.0, 0.0, 1.0]),
+        arc_angle * -1.0,
+    )
+    perp = (
+        unit_vector(np.cross(master.axis.vector.array, slave_vector))
+        * slave.diameter
+        / 2.0
+    )
 
     # Get intersect on flattened tube plane
     flat_intersect = flat_tube_intersection(master, slave)
-    flat_intersect[1] = radius_point[1]
 
-    def calc_angle(angle):
+    def pnt_intersect_at_angle(angle):
         # calculate new vector and point
         rotated_point = flat_intersect + rotate(perp, slave_vector, angle)
         pnt_intersect = plane_intersect(slave_vector, rotated_point, plane)
         pnt_intersect = np.array(pnt_intersect, dtype=float)
         return pnt_intersect
 
     for degrees in range(0, 361, angle_inc):
-        angle = degrees * math.pi / 180.
+        angle = degrees * math.pi / 180.0
         if abs(degrees - 360) < 1e-8 or degrees > 360:
             continue
-        yield calc_angle(angle)
+        yield pnt_intersect_at_angle(angle)
 
-    return calc_angle(0.0)
+    return pnt_intersect_at_angle(0.0)