From ddee9b34268f62f6293ec4f6ce386230f0a7b9c6 Mon Sep 17 00:00:00 2001
From: robbievanleeuwen <robbie.vanleeuwen@gmail.com>
Date: Tue, 7 Nov 2023 23:43:33 +1100
Subject: [PATCH 1/3] Add beam validation tests

---
 src/planestress/pre/__init__.py             |   1 +
 src/planestress/pre/geometry.py             |   4 +-
 tests/validation/test_beam.py               | 176 ++++++++++++++++++++
 tests/validation/test_validation_strand7.py |   5 +
 todo.md                                     |   2 +-
 5 files changed, 185 insertions(+), 3 deletions(-)
 create mode 100644 tests/validation/test_beam.py

diff --git a/src/planestress/pre/__init__.py b/src/planestress/pre/__init__.py
index 4f18c65..a826fc2 100644
--- a/src/planestress/pre/__init__.py
+++ b/src/planestress/pre/__init__.py
@@ -11,3 +11,4 @@
 )
 from planestress.pre.geometry import Geometry
 from planestress.pre.material import Material
+from planestress.pre.mesh import BoxField, DistanceField
diff --git a/src/planestress/pre/geometry.py b/src/planestress/pre/geometry.py
index b4cbe81..1a01739 100644
--- a/src/planestress/pre/geometry.py
+++ b/src/planestress/pre/geometry.py
@@ -1003,14 +1003,14 @@ def plot_geometry(
             # plot point tags
             if "points" in tags:
                 for pt in self.points:
-                    ax.annotate(str(pt.idx + 1), xy=(pt.x, pt.y), color="r")
+                    ax.annotate(str(pt.idx), xy=(pt.x, pt.y), color="r")
 
             # plot facet tags
             if "facets" in tags:
                 for fct in self.facets:
                     xy = (fct.pt1.x + fct.pt2.x) / 2, (fct.pt1.y + fct.pt2.y) / 2
 
-                    ax.annotate(str(fct.idx + 1), xy=xy, color="b")
+                    ax.annotate(str(fct.idx), xy=xy, color="b")
 
             # plot the analysis case
             if analysis_case is not None:
diff --git a/tests/validation/test_beam.py b/tests/validation/test_beam.py
new file mode 100644
index 0000000..d716b50
--- /dev/null
+++ b/tests/validation/test_beam.py
@@ -0,0 +1,176 @@
+"""Validation tests for simple beams in flexure."""
+
+import pytest
+import pytest_check as check
+
+from planestress.analysis import PlaneStress
+from planestress.pre import AnalysisCase, LineLoad, LineSupport, NodeLoad, NodeSupport
+from planestress.pre.library import rectangle
+
+
+@pytest.mark.parametrize("el_type", ["Tri6", "Quad8", "Quad9"])
+def test_simple_beam_point_load(el_type):
+    """Tests a simply supported beam with a point load."""
+    # set parameters
+    l = 1
+    d = 0.05
+    p = 1.0
+
+    # calculate theoretical deflection and stress
+    ixx = 1 * d**3 / 12
+    u = p * l**3 / (48 * ixx)
+    mx = p * l / 4
+    sig = mx * 0.5 * d / ixx
+
+    # create geometry
+    geom = rectangle(d=d, b=l, n_x=2)
+
+    # supports and load
+    lhs_support = NodeSupport(point=(0, 0), direction="xy")
+    rhs_support = NodeSupport(point=(l, 0), direction="y")
+    load = NodeLoad(point=(0.5 * l, d), direction="y", value=-p)
+    case = AnalysisCase([lhs_support, rhs_support, load])
+
+    # create mesh
+    ms = d / 2  # two elements thick
+
+    if el_type == "Quad9":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2)
+    elif el_type == "Quad8":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2, serendipity=True)
+    elif el_type == "Tri6":
+        geom.create_mesh(mesh_sizes=ms, mesh_order=2)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    check.almost_equal(min(res.u), -u, rel=1e-2)
+    check.almost_equal(max(res.get_nodal_stresses()[:, 0]), sig, rel=1.5e-2)
+
+
+@pytest.mark.parametrize("el_type", ["Tri6", "Quad8", "Quad9"])
+def test_simple_udl(el_type):
+    """Tests a simply supported beam with a uniformly distributed load."""
+    # set parameters
+    l = 1
+    d = 0.05
+    w = 1.0
+
+    # calculate theoretical deflection and stress
+    ixx = 1 * d**3 / 12
+    u = 5 * w * l**4 / (384 * ixx)
+    mx = w * l**2 / 8
+    sig = mx * 0.5 * d / ixx
+
+    # create geometry
+    geom = rectangle(d=d, b=l)
+
+    # supports and load
+    lhs_support = NodeSupport(point=(0, 0), direction="xy")
+    rhs_support = NodeSupport(point=(l, 0), direction="y")
+    load = LineLoad(point1=(0, d), point2=(l, d), direction="y", value=-w)
+    case = AnalysisCase([lhs_support, rhs_support, load])
+
+    # create mesh
+    ms = d / 2  # two elements thick
+
+    if el_type == "Quad9":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2)
+    elif el_type == "Quad8":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2, serendipity=True)
+    elif el_type == "Tri6":
+        geom.create_mesh(mesh_sizes=ms, mesh_order=2)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    check.almost_equal(min(res.u), -u, rel=1e-2)
+    check.almost_equal(max(res.get_nodal_stresses()[:, 0]), sig, rel=1.5e-2)
+
+
+@pytest.mark.parametrize("el_type", ["Tri6", "Quad8", "Quad9"])
+def test_cantilever_point_load(el_type):
+    """Tests a cantilever beam with a point load."""
+    # set parameters
+    l = 1
+    d = 0.05
+    p = 1.0
+
+    # calculate theoretical deflection and stress
+    ixx = 1 * d**3 / 12
+    u = p * l**3 / (3 * ixx)
+    mx = p * l
+    sig = mx * 0.5 * d / ixx
+
+    # create geometry
+    geom = rectangle(d=d, b=l)
+
+    # supports and load
+    x_support = LineSupport(point1=(0, 0), point2=(0, d), direction="x")
+    y_support = NodeSupport(point=(0, 0), direction="y")
+    load = NodeLoad(point=(l, d), direction="y", value=-p)
+    case = AnalysisCase([x_support, y_support, load])
+
+    # create mesh
+    ms = d / 2  # two elements thick
+
+    if el_type == "Quad9":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2)
+    elif el_type == "Quad8":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2, serendipity=True)
+    elif el_type == "Tri6":
+        geom.create_mesh(mesh_sizes=ms, mesh_order=2)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    check.almost_equal(min(res.u), -u, rel=1e-2)
+    check.almost_equal(max(res.get_nodal_stresses()[:, 0]), sig, rel=1.5e-2)
+
+
+@pytest.mark.parametrize("el_type", ["Tri6", "Quad8", "Quad9"])
+def test_cantilever_udl(el_type):
+    """Tests a cantilever beam with a uniformly distributed load."""
+    # set parameters
+    l = 1
+    d = 0.05
+    w = 1.0
+
+    # calculate theoretical deflection and stress
+    ixx = 1 * d**3 / 12
+    u = w * l**4 / (8 * ixx)
+    mx = w * l**2 / 2
+    sig = mx * 0.5 * d / ixx
+
+    # create geometry
+    geom = rectangle(d=d, b=l)
+
+    # supports and load
+    x_support = LineSupport(point1=(0, 0), point2=(0, d), direction="x")
+    y_support = NodeSupport(point=(0, 0), direction="y")
+    load = LineLoad(point1=(0, d), point2=(l, d), direction="y", value=-w)
+    case = AnalysisCase([x_support, y_support, load])
+
+    # create mesh
+    ms = d / 2  # two elements thick
+
+    if el_type == "Quad9":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2)
+    elif el_type == "Quad8":
+        geom.create_mesh(mesh_sizes=ms, quad_mesh=True, mesh_order=2, serendipity=True)
+    elif el_type == "Tri6":
+        geom.create_mesh(mesh_sizes=ms, mesh_order=2)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    check.almost_equal(min(res.u), -u, rel=1e-2)
+    check.almost_equal(max(res.get_nodal_stresses()[:, 0]), sig, rel=1.5e-2)
diff --git a/tests/validation/test_validation_strand7.py b/tests/validation/test_validation_strand7.py
index 0fd794a..fe2030d 100644
--- a/tests/validation/test_validation_strand7.py
+++ b/tests/validation/test_validation_strand7.py
@@ -121,3 +121,8 @@ def test_vls9(el_type):
     target_stress = -53.2
 
     check.almost_equal(target_stress, sig_yy, rel=rel)
+
+
+@pytest.mark.parametrize("el_type", ["Quad4", "Tri6", "Quad8", "Quad9"])
+def test_vls11(el_type):
+    """VLS11: xxx."""
diff --git a/todo.md b/todo.md
index 2251db9..cc1385b 100644
--- a/todo.md
+++ b/todo.md
@@ -6,6 +6,7 @@
 
 - Testing overlapping facets, other testing? (see deltares pandamesh)
 - Persistent analysis case
+- Line load normal to curves
 
 ### Analysis
 
@@ -27,7 +28,6 @@
 
 - Meshing options
   - Get mesh statistics
-- Node(?) & line load normal to curves
 - LineBC over multiple line tags?
 - Linear loads
 

From e1be7c0e8f950faf03f1d806e426996d52b9cee4 Mon Sep 17 00:00:00 2001
From: robbievanleeuwen <robbie.vanleeuwen@gmail.com>
Date: Wed, 8 Nov 2023 23:07:57 +1100
Subject: [PATCH 2/3] Fix D matrix, add normal line lines, plot embedded
 geometry, fix mesh orientation check, add st7 validation tests

---
 noxfile.py                                    |   4 +-
 .../analysis/finite_elements/lines.py         |  37 ++-
 src/planestress/post/plotting.py              |  27 +-
 src/planestress/pre/boundary_condition.py     |  31 +-
 src/planestress/pre/geometry.py               |  14 +
 src/planestress/pre/material.py               |  39 +--
 src/planestress/pre/mesh.py                   |  29 +-
 tests/validation/test_validation_strand7.py   | 277 +++++++++++++++---
 todo.md                                       |   2 +-
 9 files changed, 360 insertions(+), 100 deletions(-)

diff --git a/noxfile.py b/noxfile.py
index 8c823cf..40d361e 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -161,7 +161,9 @@ def tests(session: Session) -> None:
     Args:
         session: Nox session
     """
-    session.install(".")
+    session.run_always(
+        "poetry", "install", "--only", "main", "--extras", "pardiso", external=True
+    )
 
     # linux CI needs the no X Windows versions of gmsh
     if platform.system().lower() == "linux":
diff --git a/src/planestress/analysis/finite_elements/lines.py b/src/planestress/analysis/finite_elements/lines.py
index 2639ca3..17952d2 100644
--- a/src/planestress/analysis/finite_elements/lines.py
+++ b/src/planestress/analysis/finite_elements/lines.py
@@ -53,6 +53,30 @@ def __repr__(self) -> str:
             f"tag: {self.line_tag}."
         )
 
+    def get_normal_vector(self) -> tuple[float, float]:
+        """Gets a unit vector in the direction normal to the line.
+
+        Returns:
+            Normal unit vector.
+        """
+        dx = self.coords[0, 1] - self.coords[0, 0]
+        dy = self.coords[1, 1] - self.coords[1, 0]
+        length = (dx**2 + dy**2) ** 0.5
+
+        return -dy / length, dx / length
+
+    def get_tangent_vector(self) -> tuple[float, float]:
+        """Gets a unit vector in the direction tangent to the line.
+
+        Returns:
+            Tangent unit vector.
+        """
+        dx = self.coords[0, 1] - self.coords[0, 0]
+        dy = self.coords[1, 1] - self.coords[1, 0]
+        length = (dx**2 + dy**2) ** 0.5
+
+        return dx / length, dy / length
+
     @staticmethod
     @cache
     @njit(cache=True, nogil=True)  # type: ignore
@@ -79,7 +103,8 @@ def element_load_vector(
         """Assembles the load vector for the line element.
 
         Args:
-            direction: Direction of the line load, ``"x"``, ``"y"`` or ``"xy"``.
+            direction: Direction of the line load, ``"x"``, ``"y"``, ``"xy"``, ``"n"``
+                or ``"t"``.
             value: Value of the line load.
 
         Returns:
@@ -97,8 +122,16 @@ def element_load_vector(
             b = np.array([1.0, 0.0])
         elif direction == "y":
             b = np.array([0.0, 1.0])
-        else:
+        elif direction == "xy":
             b = np.array([1.0, 1.0])
+        elif direction == "n":
+            b = np.array(self.get_normal_vector())
+        elif direction == "t":
+            b = np.array(self.get_tangent_vector())
+        else:
+            raise ValueError(
+                f"'direction' must be 'x', 'y', 'xy', 'n' or 't', not {direction}."
+            )
 
         b *= value
 
diff --git a/src/planestress/post/plotting.py b/src/planestress/post/plotting.py
index 495ac72..d003b00 100644
--- a/src/planestress/post/plotting.py
+++ b/src/planestress/post/plotting.py
@@ -338,8 +338,31 @@ def plot_line_loads(
         y1 = line_load.point1[1]
         x2 = line_load.point2[0]
         y2 = line_load.point2[1]
-        dx = arrow_length if line_load.direction in ["x", "xy"] else 0.0
-        dy = arrow_length if line_load.direction in ["y", "xy"] else 0.0
+
+        if line_load.direction == "x":
+            dx = arrow_length
+            dy = 0.0
+        elif line_load.direction == "y":
+            dx = 0.0
+            dy = arrow_length
+        elif line_load.direction == "xy":
+            dx = arrow_length
+            dy = arrow_length
+        elif line_load.direction == "n":
+            d_x = x2 - x1
+            d_y = y2 - y1
+            length = (d_x**2 + d_y**2) ** 0.5
+            dx = -d_y / length * arrow_length
+            dy = d_x / length * arrow_length
+        elif line_load.direction == "t":
+            d_x = x2 - x1
+            d_y = y2 - y1
+            length = (d_x**2 + d_y**2) ** 0.5
+            dx = d_x / length * arrow_length
+            dy = d_y / length * arrow_length
+        else:
+            dx = 0.0
+            dy = 0.0
 
         # check to see if arrow tip is in geometry or on boundary
         pt = shapely.Point(0.5 * (x1 + x2) + dx, 0.5 * (y1 + y2) + dy)
diff --git a/src/planestress/pre/boundary_condition.py b/src/planestress/pre/boundary_condition.py
index ff68fb6..ed33846 100644
--- a/src/planestress/pre/boundary_condition.py
+++ b/src/planestress/pre/boundary_condition.py
@@ -507,8 +507,9 @@ def __init__(
         Args:
             point1: Point location (``x``, ``y``) of the start of the line load.
             point2: Point location (``x``, ``y``) of the end of the line load.
-            direction: Direction of the line load, ``"x"``, ``"y"`` or ``"xy"`` (two
-                line loads in both ``x`` and ``y`` directions).
+            direction: Direction of the line load, ``"x"``, ``"y"``, ``"xy"`` (two
+                line loads in both ``x`` and ``y`` directions), ``"n"`` (normal to the
+                line) and ``"t"`` (tangent to the line).
             value: Line load per unit length.
 
         Example:
@@ -538,12 +539,32 @@ def apply_bc(
         if self.mesh_tag is None:
             raise RuntimeError("Mesh tag is not assigned.")
 
+        # correct direction if normal or tangent
+        if self.direction in ["n", "t"]:
+            user_pt1 = self.point1
+            user_pt2 = self.point2
+            mesh_pt1 = self.mesh_tag.tagged_nodes[0]
+
+            # calculate distances to mesh_pt1
+            dist_pt1 = (
+                (user_pt1[0] - mesh_pt1.x) ** 2 + (user_pt1[1] - mesh_pt1.y) ** 2
+            ) ** 0.5
+            dist_pt2 = (
+                (user_pt2[0] - mesh_pt1.x) ** 2 + (user_pt2[1] - mesh_pt1.y) ** 2
+            ) ** 0.5
+
+            # if the points are flipped, flip direction of load
+            if dist_pt1 > dist_pt2:
+                val = -self.value
+            else:
+                val = self.value
+        else:
+            val = self.value
+
         # loop through all line elements
         for element in self.mesh_tag.elements:
             # get element load vector
-            f_el = element.element_load_vector(
-                direction=self.direction, value=self.value
-            )
+            f_el = element.element_load_vector(direction=self.direction, value=val)
 
             # get element degrees of freedom
             el_dofs = dof_map(node_idxs=tuple(element.node_idxs))
diff --git a/src/planestress/pre/geometry.py b/src/planestress/pre/geometry.py
index 1a01739..f1799b3 100644
--- a/src/planestress/pre/geometry.py
+++ b/src/planestress/pre/geometry.py
@@ -993,6 +993,20 @@ def plot_geometry(
                 )
                 label = None
 
+            # plot the embedded geometry
+            label = "Embedded Geometry"
+            for embed_geom in self.embedded_geometry:
+                if isinstance(embed_geom, Point):
+                    ax.plot(embed_geom.x, embed_geom.y, "k*", label=label)
+                elif isinstance(embed_geom, Facet):
+                    ax.plot(
+                        [embed_geom.pt1.x, embed_geom.pt2.x],
+                        [embed_geom.pt1.y, embed_geom.pt2.y],
+                        "k*-",
+                        label=label,
+                    )
+                label = None
+
             # plot the holes
             label = "Holes"
             for hl in self.holes:
diff --git a/src/planestress/pre/material.py b/src/planestress/pre/material.py
index 9310517..5288263 100644
--- a/src/planestress/pre/material.py
+++ b/src/planestress/pre/material.py
@@ -33,28 +33,6 @@ class Material:
     density: float = 1.0
     color: str = "w"
 
-    @property
-    def mu(self) -> float:
-        r"""Returns Lamé parameter mu.
-
-        Returns:
-            Lamé parameter :math:`\mu`.
-        """
-        return self.elastic_modulus / (2 * (1 + self.poissons_ratio))
-
-    @property
-    def lda(self) -> float:
-        r"""Returns Lamé parameter lambda.
-
-        Returns:
-            Lamé parameter :math:`\lambda`.
-        """
-        return (
-            self.poissons_ratio
-            * self.elastic_modulus
-            / ((1 + self.poissons_ratio) * (1 - 2 * self.poissons_ratio))
-        )
-
     @cache
     def get_d_matrix(self) -> npt.NDArray[np.float64]:
         r"""Returns the constitutive matrix for plane-stress.
@@ -67,16 +45,15 @@ def get_d_matrix(self) -> npt.NDArray[np.float64]:
         Returns:
             Constitutive matrix.
         """
-        mu = self.mu
-        lda = self.lda
-
-        d_mat = np.zeros((3, 3))  # allocate D matrix
-        d_mat[0:2, 0:2] = lda + 2 * mu
-        d_mat[2, 2] = mu
-        d_mat[0, 1] = lda
-        d_mat[1, 0] = lda
+        d_mat = np.array(
+            [
+                [1, self.poissons_ratio, 0],
+                [self.poissons_ratio, 1, 0],
+                [0, 0, (1 - self.poissons_ratio) / 2],
+            ]
+        )
 
-        return d_mat
+        return d_mat * self.elastic_modulus / (1 - self.poissons_ratio**2)
 
 
 DEFAULT_MATERIAL = Material()
diff --git a/src/planestress/pre/mesh.py b/src/planestress/pre/mesh.py
index c526ca5..f765d35 100644
--- a/src/planestress/pre/mesh.py
+++ b/src/planestress/pre/mesh.py
@@ -215,19 +215,6 @@ def create_mesh(
                     dim=1, tags=[ln], inDim=2, inTag=geo.pt1.poly_idxs[0]
                 )
 
-        # check surface orientation:
-        # list describing if surface is correctly oriented
-        surface_orientated: list[bool] = []
-
-        for _, tag in gmsh.model.get_entities(dim=2):
-            normal = gmsh.model.get_normal(tag, (0, 0))
-
-            # if surface is incorrectly oriented, re-orient
-            if normal[2] < 0:
-                surface_orientated.append(False)
-            else:
-                surface_orientated.append(True)
-
         # calculate bounding box
         self.bbox = gmsh.model.get_bounding_box(dim=-1, tag=-1)
 
@@ -253,6 +240,22 @@ def create_mesh(
         # set mesh order
         gmsh.model.mesh.set_order(order=mesh_order)
 
+        # check surface orientation:
+        # list describing if surface is correctly oriented
+        surface_orientated: list[bool] = []
+
+        for _, tag in gmsh.model.get_entities(dim=2):
+            tags, coord, param = gmsh.model.mesh.getNodes(
+                dim=2, tag=tag, includeBoundary=True
+            )
+            normal = gmsh.model.get_normal(tag=tag, parametricCoord=param[0:2])
+
+            # if surface is incorrectly oriented, re-orient
+            if normal[2] < 0:
+                surface_orientated.append(False)
+            else:
+                surface_orientated.append(True)
+
         # view model - for debugging
         # gmsh.fltk.run()
 
diff --git a/tests/validation/test_validation_strand7.py b/tests/validation/test_validation_strand7.py
index fe2030d..8785e1f 100644
--- a/tests/validation/test_validation_strand7.py
+++ b/tests/validation/test_validation_strand7.py
@@ -3,17 +3,76 @@
 Insert reference here...
 """
 
+from typing import Callable
+
 import numpy as np
 import pytest
 import pytest_check as check
+from shapely import Polygon
 
 import planestress.pre.boundary_condition as bc
-from planestress.analysis.plane_stress import PlaneStress
-from planestress.pre.analysis_case import AnalysisCase
+from planestress.analysis import PlaneStress
+from planestress.pre import AnalysisCase, Geometry, Material
 from planestress.pre.library import circle, rectangle, steel_material
 
 
-def test_vls1():
+@pytest.fixture
+def circular_membrane() -> Callable:
+    """Creates a meshed circular membrane for VLS8 & VLS9.
+
+    Returns:
+        Generator function, returning a ``Geometry`` object.
+    """
+
+    def _generate(el_type: str) -> Geometry:
+        """Generates the unit square.
+
+        Args:
+            lc: Characterisic mesh length.
+            element_type: Element type, can be ``"Tri3"``, ``"Tri6"``, ``"Quad4"``,
+                ``"Quad8"`` or ``"Quad9"``.
+
+        Returns:
+            Geometry object.
+        """
+        # define materials - use N & mm
+        steel = steel_material(thickness=1000.0)
+
+        # define geometry
+        circle_outer = circle(r=11e3, n=128, material=steel)
+        circle_inner = circle(r=10e3, n=128)
+        bbox = rectangle(d=12e3, b=12e3)  # bounding box
+        geom = (circle_outer - circle_inner) & bbox
+
+        # create mesh
+        if el_type == "Quad4":
+            geom.create_mesh(
+                mesh_sizes=1000.0, quad_mesh=True, mesh_order=1, mesh_algorithm=11
+            )
+        elif el_type == "Tri6":
+            geom.create_mesh(mesh_sizes=1000.0, mesh_order=2)
+        elif el_type == "Quad8":
+            geom.create_mesh(
+                mesh_sizes=1000.0,
+                quad_mesh=True,
+                mesh_order=2,
+                serendipity=True,
+                mesh_algorithm=11,
+            )
+        elif el_type == "Quad9":
+            geom.create_mesh(
+                mesh_sizes=1000.0, quad_mesh=True, mesh_order=2, mesh_algorithm=11
+            )
+        else:
+            raise ValueError(f"{el_type} element not supported for this test.")
+
+        return geom
+
+    return _generate
+
+
+@pytest.mark.parametrize("el_type", ["Quad4", "Tri6", "Quad8", "Quad9"])
+def test_vls1(el_type):
     """VLS1: Elliptic Membrane.
 
     An elliptical plate with an elliptical hole is analysed.
@@ -27,13 +86,81 @@ def test_vls1():
     Materials: Steel, E=200e3, v=0.3
 
     Target value - tangential stress at (x=2, y=0) of 92.7 MPa.
-
-    TODO - must first implement line load normal to curve.
     """
-    pass
+    rel = 0.03  # aim for 3% error
+
+    # define materials - use N & mm
+    steel = steel_material(thickness=1)
+
+    # create geometry
+    shell = []
+    n_curve = 24
+
+    # inner curve
+    a = 2.0
+    b = 1.0
+
+    for idx in range(n_curve):
+        theta = idx / (n_curve - 1) * np.pi / 2
+        shell.append((a * np.cos(theta), b * np.sin(theta)))
+
+    # outer curve
+    a = 3.25
+    b = 2.75
+
+    for idx in range(n_curve):
+        theta = np.pi / 2 - idx / (n_curve - 1) * np.pi / 2
+        shell.append((a * np.cos(theta), b * np.sin(theta)))
+
+    geom = Geometry(polygons=Polygon(shell=shell), materials=steel)
+
+    # create loads
+    loads = []
+    sig = 10
+
+    for idx in range(n_curve - 1):
+        theta1 = idx / (n_curve - 1) * np.pi / 2
+        theta2 = (idx + 1) / (n_curve - 1) * np.pi / 2
+        pt1 = (a * np.cos(theta1), b * np.sin(theta1))
+        pt2 = (a * np.cos(theta2), b * np.sin(theta2))
+        loads.append(bc.LineLoad(pt1, pt2, "n", -sig))
+
+    # create supports
+    lhs_support = bc.LineSupport(point1=(0.0, 1), point2=(0.0, 2.75), direction="x")
+    rhs_support = bc.LineSupport(point1=(2, 0.0), point2=(3.25, 0.0), direction="y")
+    bcs = [lhs_support, rhs_support]
+    bcs.extend(loads)
+    case = AnalysisCase(bcs)
+
+    # create mesh
+    if el_type == "Quad4":
+        geom.create_mesh(quad_mesh=True, mesh_order=1, mesh_algorithm=11)
+    elif el_type == "Tri6":
+        geom.create_mesh(mesh_sizes=0.15, mesh_order=2)
+    elif el_type == "Quad8":
+        geom.create_mesh(
+            quad_mesh=True, mesh_order=2, serendipity=True, mesh_algorithm=8
+        )
+    elif el_type == "Quad9":
+        geom.create_mesh(quad_mesh=True, mesh_order=2, mesh_algorithm=8)
+    else:
+        raise ValueError(f"{el_type} element not supported for this test.")
 
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
 
-def test_vls8():
+    # get stress at (x=2, y=0)
+    node_idx = ps.mesh.get_node_idx_by_coordinates(2.0, 0.0)
+    sig_yy = res.get_nodal_stresses()[node_idx][1]
+    target_stress = 92.7
+
+    check.almost_equal(target_stress, sig_yy, rel=rel)
+
+
+@pytest.mark.parametrize("el_type", ["Quad4", "Tri6", "Quad8", "Quad9"])
+def test_vls8(el_type, circular_membrane):
     """VLS8: Circular Membrane - Edge Pressure.
 
     A ring under uniform external pressure of 100 MPa is analysed.
@@ -46,14 +173,46 @@ def test_vls8():
     Materials: Steel, E=200e3, v=0.3
 
     Target value - tangential stress at (x=10, y=0) of -1150 MPa.
-
-    TODO - must first implement line load normal to curve.
     """
-    pass
+    rel = 0.03  # aim for 3% error
+
+    # create meshed geometry
+    geom = circular_membrane(el_type)
+
+    # create supports
+    lhs_support = bc.LineSupport(point1=(0.0, 10e3), point2=(0.0, 11e3), direction="x")
+    rhs_support = bc.LineSupport(point1=(10e3, 0.0), point2=(11e3, 0.0), direction="y")
+
+    # create loads
+    loads = []
+    sig = 100  # MPa
+    p = sig * 1000  # N/mm
+    for idx in range(32):
+        theta1 = np.pi / 64 * idx
+        theta2 = np.pi / 64 * (idx + 1)
+        pt1 = (11e3 * np.cos(theta1), 11e3 * np.sin(theta1))
+        pt2 = (11e3 * np.cos(theta2), 11e3 * np.sin(theta2))
+        loads.append(bc.LineLoad(pt1, pt2, "n", p))
+
+    bcs = [lhs_support, rhs_support]
+    bcs.extend(loads)
+    case = AnalysisCase(bcs)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    # get stress at (x=10, y=0)
+    node_idx = ps.mesh.get_node_idx_by_coordinates(10e3, 0.0)
+    sig_yy = res.get_nodal_stresses()[node_idx][1]
+    target_stress = -1150.0
+
+    check.almost_equal(target_stress, sig_yy, rel=rel)
 
 
 @pytest.mark.parametrize("el_type", ["Quad4", "Tri6", "Quad8", "Quad9"])
-def test_vls9(el_type):
+def test_vls9(el_type, circular_membrane):
     """VLS9: Circular Membrane - Point Load.
 
     A ring under concentrated forces is analysed. (10000 kN at 4 x 45 deg).
@@ -69,14 +228,8 @@ def test_vls9(el_type):
     """
     rel = 0.03  # aim for 3% error
 
-    # define materials - use N & mm
-    steel = steel_material(thickness=1000.0)
-
-    # define geometry
-    circle_outer = circle(r=11e3, n=128, material=steel)
-    circle_inner = circle(r=10e3, n=128)
-    bbox = rectangle(d=12e3, b=12e3)  # bounding box
-    geom = (circle_outer - circle_inner) & bbox
+    # create meshed geometry
+    geom = circular_membrane(el_type)
 
     # create supports
     lhs_support = bc.LineSupport(point1=(0.0, 10e3), point2=(0.0, 11e3), direction="x")
@@ -88,31 +241,9 @@ def test_vls9(el_type):
     load = bc.NodeLoad(point=pt, direction="xy", value=-force)
     case = AnalysisCase([lhs_support, rhs_support, load])
 
-    # create mesh
-    if el_type == "Quad4":
-        geom.create_mesh(
-            mesh_sizes=1000.0, quad_mesh=True, mesh_order=1, mesh_algorithm=11
-        )
-    elif el_type == "Tri6":
-        geom.create_mesh(mesh_sizes=1000.0, mesh_order=2)
-    elif el_type == "Quad8":
-        geom.create_mesh(
-            mesh_sizes=1000.0,
-            quad_mesh=True,
-            mesh_order=2,
-            serendipity=True,
-            mesh_algorithm=11,
-        )
-    elif el_type == "Quad9":
-        geom.create_mesh(
-            mesh_sizes=1000.0, quad_mesh=True, mesh_order=2, mesh_algorithm=11
-        )
-    else:
-        raise ValueError(f"{el_type} element not supported for this test.")
-
     # solve
     ps = PlaneStress(geom, [case])
-    results_list = ps.solve()
+    results_list = ps.solve(solver_type="pardiso")
     res = results_list[0]
 
     # get stress at (x=10, y=0)
@@ -123,6 +254,62 @@ def test_vls9(el_type):
     check.almost_equal(target_stress, sig_yy, rel=rel)
 
 
-@pytest.mark.parametrize("el_type", ["Quad4", "Tri6", "Quad8", "Quad9"])
-def test_vls11(el_type):
-    """VLS11: xxx."""
+def test_vls11():
+    """VLS11: Plate Patch Test.
+
+    Enforced displacements are applied to acheive a uniform strain of 1e-3.
+
+    Material properties: E = 1e6, nu = 0.25, t = 0.001.
+
+    sig_xx = sig_yy = 1333
+    sig_xy = 400
+    """
+    # create material
+    material = Material(elastic_modulus=1e6, poissons_ratio=0.25, thickness=0.001)
+
+    # create geometry
+    geom = rectangle(d=0.12, b=0.24, material=material)
+
+    # add nodes
+    geom.embed_point(x=0.04, y=0.02)
+    geom.embed_point(x=0.08, y=0.08)
+    geom.embed_point(x=0.16, y=0.08)
+    geom.embed_point(x=0.18, y=0.03)
+
+    # create mesh
+    geom.create_mesh(mesh_sizes=0.24)
+
+    # apply loads
+    loads = []
+    pts = [
+        (0, 0),
+        (0, 0.12),
+        (0.24, 0.12),
+        (0.24, 0),
+        (0.04, 0.02),
+        (0.08, 0.08),
+        (0.16, 0.08),
+        (0.18, 0.03),
+    ]
+
+    for pt in pts:
+        dx = 1e-3 * (pt[0] + 0.5 * pt[1])
+        dy = 1e-3 * (pt[1] + 0.5 * pt[0])
+        loads.append(bc.NodeSupport(point=pt, direction="x", value=dx))
+        loads.append(bc.NodeSupport(point=pt, direction="y", value=dy))
+
+    case = AnalysisCase(loads)
+
+    # solve
+    ps = PlaneStress(geom, [case])
+    results_list = ps.solve(solver_type="pardiso")
+    res = results_list[0]
+
+    # check stresses
+    sig_xx = res.get_nodal_stresses()[0][0]
+    sig_yy = res.get_nodal_stresses()[0][1]
+    sig_xy = res.get_nodal_stresses()[0][2]
+
+    check.almost_equal(4000 / 3, sig_xx)
+    check.almost_equal(4000 / 3, sig_yy)
+    check.almost_equal(400, sig_xy)
diff --git a/todo.md b/todo.md
index cc1385b..bf07f2a 100644
--- a/todo.md
+++ b/todo.md
@@ -6,7 +6,7 @@
 
 - Testing overlapping facets, other testing? (see deltares pandamesh)
 - Persistent analysis case
-- Line load normal to curves
+- Rethink boundary conditions (rather than assign to nodes, assign to tags!)
 
 ### Analysis
 

From 79795978631c62f788939afa028b24c5ef814375 Mon Sep 17 00:00:00 2001
From: robbievanleeuwen <robbie.vanleeuwen@gmail.com>
Date: Wed, 8 Nov 2023 23:21:38 +1100
Subject: [PATCH 3/3] Recent mkl only has mac distros

---
 poetry.lock    | 12 ++++++++----
 pyproject.toml |  3 ++-
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/poetry.lock b/poetry.lock
index 4b491c4..5716d1f 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1879,12 +1879,16 @@ files = [
 
 [[package]]
 name = "mkl"
-version = "2023.2.1"
+version = "2023.2.0"
 description = "Intel® oneAPI Math Kernel Library"
 optional = true
 python-versions = "*"
 files = [
-    {file = "mkl-2023.2.1-py2.py3-none-macosx_10_15_x86_64.macosx_11_0_x86_64.whl", hash = "sha256:953a55a778d87cea9fd44f777c3b0a7a5804e9e810dc5f22852b2c4a42a79bcb"},
+    {file = "mkl-2023.2.0-py2.py3-none-macosx_10_15_x86_64.macosx_11_0_x86_64.whl", hash = "sha256:a8e448a49cb636afc589c952b5b3de6c94ea19fe3c50e6c4e7e915f47adf3070"},
+    {file = "mkl-2023.2.0-py2.py3-none-manylinux1_i686.whl", hash = "sha256:60e301ee5f836c00cc48bbcef803c57032e15ebdeffdf78190e828c9adf41206"},
+    {file = "mkl-2023.2.0-py2.py3-none-manylinux1_x86_64.whl", hash = "sha256:577ac0fd075f87351d989d141f46efcb462bd2de9851ffe1c0cae3fa315b0ff7"},
+    {file = "mkl-2023.2.0-py2.py3-none-win32.whl", hash = "sha256:2cf51d9e2104f874f48a5e06e4dd3f07e756773f118ad555bb489edb91adef23"},
+    {file = "mkl-2023.2.0-py2.py3-none-win_amd64.whl", hash = "sha256:7b93b9a20c51e8c60f891097924f3f2be786797d4b587dc848725efacd559a0b"},
 ]
 
 [package.dependencies]
@@ -3888,9 +3892,9 @@ docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.link
 testing = ["big-O", "jaraco.functools", "jaraco.itertools", "more-itertools", "pytest (>=6)", "pytest-black (>=0.3.7)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-ignore-flaky", "pytest-mypy (>=0.9.1)", "pytest-ruff"]
 
 [extras]
-pardiso = ["intel-openmp", "pypardiso"]
+pardiso = ["intel-openmp", "mkl", "pypardiso"]
 
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.9.0,<3.12"
-content-hash = "b44641fb1014dfc4d2c61caefe3dd88143ff6db03bb8fdce7d784af085d37c79"
+content-hash = "f151da1cdd1df9add34d588759d4386ccbc0af490e1718fb28b31bf487dc23dd"
diff --git a/pyproject.toml b/pyproject.toml
index 23c175d..38944ff 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -51,6 +51,7 @@ rich = "^13.5.0"
 click = "^8.1.7"
 pypardiso = { version = "^0.4.3", optional = true }
 intel-openmp = { version = "==2023.2.0", optional = true }
+mkl = { version = "==2023.2.0", optional = true }
 
 [tool.poetry.group.dev.dependencies]
 black = "^23.10.1"
@@ -85,7 +86,7 @@ sphinx-copybutton = "^0.5.2"
 sphinxext-opengraph = "^0.8.2"
 
 [tool.poetry.extras]
-pardiso = ["pypardiso", "intel-openmp"]
+pardiso = ["pypardiso", "intel-openmp", "mkl"]
 
 [tool.poetry.scripts]
 planestress = "planestress.__main__:main"