More validation tests, add normal line loads

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":

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"

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
 

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)

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

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))

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:
@@ -1003,14 +1017,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:

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()

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()