Feature: Tutorials for PostDeform

pyhope/mesh/transform/mesh_transform.py

@@ -158,6 +158,7 @@ def TransformMesh() -> None:
 
     # Read in the mesh post-deformation flag
     meshPostDeform = GetStr('MeshPostDeform') if CountOption('MeshPostDeform') != 0 else 'none'
+    meshPostDeform = meshPostDeform.lower()
 
     # Leave if no transformation is required
     if all(x == 0 for x in [nMeshScale, nMeshTrans, nMeshRot]) and meshPostDeform == 'none':

pyhope/mesh/transform/templates/cylinder.py

@@ -43,72 +43,90 @@ def PostDeform(points: np.ndarray) -> np.ndarray:
     This function applies a deformation transformation to the input points based on the given Fortran logic.
     The transformation maps a 2D square region to a cylindrical or toroidal coordinate system.
     """
+    # Local imports ----------------------------------------
+    import pyhope.output.output as hopout
+    from pyhope.readintools.readintools import CreateReal, GetReal
+    # ------------------------------------------------------
 
-    # TODO: Readin parameters from a configuration file
-    PostDeform_R0 = 1.0       # Define appropriate values
-    PostDeform_Rtorus = -1.0  # Adjust as needed
-    PostDeform_Lz = 1.0       # Cylinder height scale
-    PostDeform_sq = 0.0       # Spiral rotation parameter
-    MeshPostDeform = 1        # Deformation mode
+    hopout.sep()
+    CreateReal('PostDeform_R0',     default= 1.0, help='Scaling factor for the cylinder/torus radius')                         # noqa: E251
+    CreateReal('PostDeform_sq',     default= 0.0, help='Spiral factor along z. 0: no spiral, 1: one rotation for z in [0,1]')  # noqa: E251
+    CreateReal('PostDeform_Rtorus', default=-1.0, help='>0 for torus major radius (around z); <0 for cylinder mode')           # noqa: E251
+    CreateReal('PostDeform_Lz',     default= 1.0, help='Axial scaling for cylinder mode')                                      # noqa: E251
+
+    PostDeform_R0     = GetReal('PostDeform_R0')
+    PostDeform_sq     = GetReal('PostDeform_sq')
+    PostDeform_Rtorus = GetReal('PostDeform_Rtorus')
+    PostDeform_Lz     = GetReal('PostDeform_Lz')
 
     nTotal = points.shape[0]
-    X_out = np.zeros_like(points)
+    X_out  = np.zeros_like(points)
+    Pi     = np.pi
 
     for i in range(nTotal):
-        x = points[i, :].copy()
+        x  = points[i, :].copy()
         rr = max(abs(x[0]), abs(x[1]))
 
         if rr < 0.5:
-            dx1 = np.array([
-                0.5 * np.sqrt(2) * np.cos(0.25 * np.pi * x[1] / 0.5) - 0.5,
-                0.5 * np.sqrt(2) * np.sin(0.25 * np.pi * x[1] / 0.5) - x[1]
-            ])
-            dx2 = np.array([
-                0.5 * np.sqrt(2) * np.sin(0.25 * np.pi * x[0] / 0.5) - x[0],
-                0.5 * np.sqrt(2) * np.cos(0.25 * np.pi * x[0] / 0.5) - 0.5
-            ])
+            # inside [-0.5,0.5]^2
+            # right side at x=0.5
+            dx1_0 = 0.5 * np.sqrt(2.0) * np.cos(0.25 * Pi * x[1] / 0.5) - 0.5
+            dx1_1 = 0.5 * np.sqrt(2.0) * np.sin(0.25 * Pi * x[1] / 0.5) - x[1]
+            dx1 = np.array([dx1_0, dx1_1])
+
+            # upper side at y=0.5
+            dx2_0 = 0.5 * np.sqrt(2.0) * np.sin(0.25 * Pi * x[0] / 0.5) - x[0]
+            dx2_1 = 0.5 * np.sqrt(2.0) * np.cos(0.25 * Pi * x[0] / 0.5) - 0.5
+            dx2 = np.array([dx2_0, dx2_1])
+
             alpha = 0.35
-            dx = alpha * (dx1 * np.array([2 * x[0], 1.]) + dx2 * np.array([1., 2 * x[1]]))
+            # coons mapping of edges, dx=0 at the corners
+            dx = alpha * (dx1 * np.array([2.0 * x[0], 1.0]) + dx2 * np.array([1.0, 2.0 * x[1]]))
         else:
-            if abs(x[1]) < abs(x[0]):
-                dx = np.array([
-                    x[0] * np.sqrt(2) * np.cos(0.25 * np.pi * x[1] / x[0]) - x[0],
-                    x[0] * np.sqrt(2) * np.sin(0.25 * np.pi * x[1] / x[0]) - x[1]
-                ])
-            else:
-                dx = np.array([
-                    x[1] * np.sqrt(2) * np.sin(0.25 * np.pi * x[0] / x[1]) - x[0],
-                    x[1] * np.sqrt(2) * np.cos(0.25 * np.pi * x[0] / x[1]) - x[1]
-                ])
-            alpha = min(1., 2. * rr - 1.)
-            alpha = np.sin(0.5 * np.pi * alpha)
-            alpha = 1.0 * alpha + 0.35 * (1. - alpha)
+            # outside [-0.5,0.5]^2
+            if abs(x[1]) < abs(x[0]):  # left and right quarter
+                dx0 = x[0] * np.sqrt(2.0) * np.cos(0.25 * Pi * x[1] / x[0]) - x[0]
+                dx1 = x[0] * np.sqrt(2.0) * np.sin(0.25 * Pi * x[1] / x[0]) - x[1]
+                dx = np.array([dx0, dx1])
+            else:  # upper and lower quarter (ABS(x(2)).GE.ABS(x(1)))
+                dx0 = x[1] * np.sqrt(2.0) * np.sin(0.25 * Pi * x[0] / x[1]) - x[0]
+                dx1 = x[1] * np.sqrt(2.0) * np.cos(0.25 * Pi * x[0] / x[1]) - x[1]
+                dx = np.array([dx0, dx1])
+
+            # maps [0.5,1] --> [0,1] and alpha=1 outside [-1,1]^2
+            alpha = min(1.0, 2.0 * rr - 1.0)
+            # smooth transition at the outer boundary max(|x|,|y|)=1
+            alpha = np.sin(0.5 * Pi * alpha)
+            # alpha=1 at max(|x|,|y|)=1, and alpha=0.35 at max(|x|,|y|)=0.5
+            alpha = 1.0 * alpha + 0.35 * (1.0 - alpha)
             dx *= alpha
 
-        xout = PostDeform_R0 * np.sqrt(0.5) * (x[:2] + dx)
+        xout     = np.zeros(3, dtype=points.dtype)
+        xout[:2] = PostDeform_R0 * np.sqrt(0.5) * (x[:2] + dx)
 
-        if MeshPostDeform == 1:
-            arg = 2. * np.pi * x[2] * PostDeform_sq
-        elif MeshPostDeform == 11:
-            arg = 2. * np.pi * x[2] * PostDeform_sq * np.sum(xout**2)
-        elif MeshPostDeform == 12:
-            arg = 2. * np.pi * x[2] * PostDeform_sq * np.sum(xout**2) * (1 + 0.5 * xout[0])
-        else:
-            arg = 0
+        # Spiral rotation along z sq=0. no spiral, sq=1: 1 rotation along z [0,1]
+        arg = 2.0 * Pi * x[2] * PostDeform_sq
 
-        rotmat = np.array([
-            [np.cos(arg), -np.sin(arg)],
-            [np.sin(arg),  np.cos(arg)]
-        ])
-        xout = np.matmul(rotmat, xout)
+        # Rotation matrix
+        c = np.cos(arg)
+        s = np.sin(arg)
+        rotmat   = np.array([[c, -s],
+                             [s,  c]], dtype=xout.dtype)
+        xout[:2] = rotmat @ xout[:2]
 
-        if PostDeform_Rtorus < 0:
-            xout = np.append(xout, x[2] * PostDeform_Lz)
+        if PostDeform_Rtorus < 0.0:
+            # cylinder
+            xout[2] = x[2] * PostDeform_Lz
         else:
-            temp_z = xout[1]
-            xout[1] = -(xout[0] + PostDeform_Rtorus) * np.sin(2 * np.pi * x[2])
-            xout[0] = (xout[0] + PostDeform_Rtorus) * np.cos(2 * np.pi * x[2])
-            xout = np.append(xout, temp_z)
+            # torus, z_in must be [0,1] and periodic
+            # torus around z axis ,x =R*cos(phi), y=-R*sin(phi)
+            # store Z
+            z_val = xout[1]
+            R_minor = xout[0] + PostDeform_Rtorus
+            phi = 2.0 * Pi * x[2]
+            xout[2] = z_val
+            xout[1] = -R_minor * np.sin(phi)
+            xout[0] =  R_minor * np.cos(phi)
 
         X_out[i, :] = xout
 

pyhope/mesh/transform/templates/sphere.py

@@ -44,17 +44,21 @@ def PostDeform(points: np.ndarray) -> np.ndarray:
     3D box, x,y in [-1,1]^3, to Sphere with radius PostDeform_R0
     all points outside [-1,1]^4 will be mapped directly to a sphere
     """
+    # Local imports ----------------------------------------
+    import pyhope.output.output as hopout
+    from pyhope.readintools.readintools import CreateReal, GetReal
+    # ------------------------------------------------------
 
-    # TODO: Readin parameters from a configuration file
-    PostDeform_R0 = 1.0  # Define based on the expected scaling factor
-
+    hopout.sep()
+    CreateReal('PostDeform_R0', default=1.0, help='Scaling factor for the radius of the sphere')
+    PostDeform_R0 = GetReal('PostDeform_R0')
 
     nTotal = points.shape[0]
-    X_out = np.zeros_like(points)
-    Pi = np.pi
+    X_out  = np.zeros_like(points)
+    Pi     = np.pi
 
     for i in range(nTotal):
-        x = points[i, :].copy()
+        x  = points[i, :].copy()
         rr = max(abs(x[0]), abs(x[1]), abs(x[2]))
 
         if rr <= 0.5:

pyhope/mesh/transform/templates/torus.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+# This file is part of PyHOPE
+#
+# Copyright (c) 2024 Numerics Research Group, University of Stuttgart, Prof. Andrea Beck
+#
+# PyHOPE is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# PyHOPE is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# PyHOPE. If not, see <http://www.gnu.org/licenses/>.
+
+# ==================================================================================================================================
+# Mesh generation library
+# ==================================================================================================================================
+# ----------------------------------------------------------------------------------------------------------------------------------
+# Standard libraries
+# ----------------------------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------------------------
+# Third-party libraries
+# ----------------------------------------------------------------------------------------------------------------------------------
+import numpy as np
+# ----------------------------------------------------------------------------------------------------------------------------------
+# Local imports
+# ----------------------------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------------------------
+# Local definitions
+# ----------------------------------------------------------------------------------------------------------------------------------
+# ==================================================================================================================================
+
+
+def PostDeform(points: np.ndarray) -> np.ndarray:
+    """
+    Apply post-deformation transformation to input points.
+
+    2D box, x,y in [-1,1]^2, to cylinder with radius PostDeform_R0 (with PostDeform_Rtorus>0 to a torus, with zperiodic [0,1])
+    all points outside [-1,1]^2 will be mapped directly to a circle (p.e. 2,2 => sqrt(0.5)*PostDeform_R0*(2,2) )
+    inside [-1,1]^2 and outside [-0.5,0.5]^2 there will be a blending from a circle to a square
+    the inner square [-0.5,0.5]^2 will be a linear blending of the bounding curves
+    """
+    # Local imports ----------------------------------------
+    import pyhope.output.output as hopout
+    from pyhope.readintools.readintools import CreateReal, GetReal
+    # ------------------------------------------------------
+
+    hopout.sep()
+    CreateReal('PostDeform_R0',     default= 1.0, help='Scaling factor for the cylinder/torus radius')                         # noqa: E251
+    CreateReal('PostDeform_sq',     default= 0.0, help='Spiral factor along z. 0: no spiral, 1: one rotation for z in [0,1]')  # noqa: E251
+    CreateReal('PostDeform_Rtorus', default=-1.0, help='>0 for torus major radius (around z); <0 for cylinder mode')           # noqa: E251
+    CreateReal('PostDeform_Lz',     default= 1.0, help='Axial scaling for cylinder mode')                                      # noqa: E251
+
+    PostDeform_R0     = GetReal('PostDeform_R0')
+    PostDeform_sq     = GetReal('PostDeform_sq')
+    PostDeform_Rtorus = GetReal('PostDeform_Rtorus')
+    PostDeform_Lz     = GetReal('PostDeform_Lz')
+
+    nTotal = points.shape[0]
+    X_out  = np.zeros_like(points)
+    Pi     = np.pi
+
+    for i in range(nTotal):
+        x  = points[i, :].copy()
+        rr = max(abs(x[0]), abs(x[1]))
+
+        if rr < 0.5:
+            # inside [-0.5,0.5]^2
+            # right side at x=0.5
+            dx1_0 = 0.5 * np.sqrt(2.0) * np.cos(0.25 * Pi * x[1] / 0.5) - 0.5
+            dx1_1 = 0.5 * np.sqrt(2.0) * np.sin(0.25 * Pi * x[1] / 0.5) - x[1]
+            dx1 = np.array([dx1_0, dx1_1])
+
+            # upper side at y=0.5
+            dx2_0 = 0.5 * np.sqrt(2.0) * np.sin(0.25 * Pi * x[0] / 0.5) - x[0]
+            dx2_1 = 0.5 * np.sqrt(2.0) * np.cos(0.25 * Pi * x[0] / 0.5) - 0.5
+            dx2 = np.array([dx2_0, dx2_1])
+
+            alpha = 0.35
+            # coons mapping of edges, dx=0 at the corners
+            dx = alpha * (dx1 * np.array([2.0 * x[0], 1.0]) + dx2 * np.array([1.0, 2.0 * x[1]]))
+        else:
+            # outside [-0.5,0.5]^2
+            if abs(x[1]) < abs(x[0]):  # left and right quarter
+                dx0 = x[0] * np.sqrt(2.0) * np.cos(0.25 * Pi * x[1] / x[0]) - x[0]
+                dx1 = x[0] * np.sqrt(2.0) * np.sin(0.25 * Pi * x[1] / x[0]) - x[1]
+                dx = np.array([dx0, dx1])
+            else:  # upper and lower quarter (ABS(x(2)).GE.ABS(x(1)))
+                dx0 = x[1] * np.sqrt(2.0) * np.sin(0.25 * Pi * x[0] / x[1]) - x[0]
+                dx1 = x[1] * np.sqrt(2.0) * np.cos(0.25 * Pi * x[0] / x[1]) - x[1]
+                dx = np.array([dx0, dx1])
+
+            # maps [0.5,1] --> [0,1] and alpha=1 outside [-1,1]^2
+            alpha = min(1.0, 2.0 * rr - 1.0)
+            # smooth transition at the outer boundary max(|x|,|y|)=1
+            alpha = np.sin(0.5 * Pi * alpha)
+            # alpha=1 at max(|x|,|y|)=1, and alpha=0.35 at max(|x|,|y|)=0.5
+            alpha = 1.0 * alpha + 0.35 * (1.0 - alpha)
+            dx *= alpha
+
+        xout     = np.zeros(3, dtype=points.dtype)
+        xout[:2] = PostDeform_R0 * np.sqrt(0.5) * (x[:2] + dx)
+
+        # Spiral rotation along z sq=0. no spiral, sq=1: 1 rotation along z [0,1]
+        arg = 2.0 * Pi * x[2] * PostDeform_sq
+
+        # Rotation matrix
+        c = np.cos(arg)
+        s = np.sin(arg)
+        rotmat   = np.array([[c, -s],
+                             [s,  c]], dtype=xout.dtype)
+        xout[:2] = rotmat @ xout[:2]
+
+        if PostDeform_Rtorus < 0.0:
+            # cylinder
+            xout[2] = x[2] * PostDeform_Lz
+        else:
+            # torus, z_in must be [0,1] and periodic
+            # torus around z axis ,x =R*cos(phi), y=-R*sin(phi)
+            # store Z
+            z_val = xout[1]
+            R_minor = xout[0] + PostDeform_Rtorus
+            phi = 2.0 * Pi * x[2]
+            xout[2] = z_val
+            xout[1] = -R_minor * np.sin(phi)
+            xout[0] =  R_minor * np.cos(phi)
+
+        X_out[i, :] = xout
+
+    return X_out

pyhope/readintools/readintools.py

@@ -380,7 +380,7 @@ def GetStr(name: str, default: Optional[str] = None, number: Optional[int] = Non
 
 def GetReal(name: str, default: Optional[str] = None, number: Optional[int] = None) -> float:
     value = GetParam(name=name, default=default, number=number, calltype='real')
-    return strToFloatOrPi(value)
+    return strToFloatOrPi(str(value))
 
 
 def GetInt(name: str, default: Optional[str] = None, number: Optional[int] = None) -> int:

tutorials/1-06-curved_cylinder/analyze.toml

@@ -0,0 +1,23 @@
+[ElemInfo]
+min = 0.0
+max = 10000.0
+mean = 1781.8333333333333
+stddev = 2823.094691803463
+
+[GlobalNodeIDs]
+min = 1.0
+max = 5729.0
+mean = 2861.0498
+stddev = 1623.2432874710926
+
+[NodeCoords]
+min = -2.0000000000000004
+max = 20.0
+mean = 3.333333333330626
+stddev = 5.8444859976977375
+
+[SideInfo]
+min = -255.0
+max = 256.0
+mean = 17.120833333333334
+stddev = 69.01824323543417