diff --git a/pyhope/mesh/transform/mesh_transform.py b/pyhope/mesh/transform/mesh_transform.py index 488974c6..fd363093 100644 --- a/pyhope/mesh/transform/mesh_transform.py +++ b/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': diff --git a/pyhope/mesh/transform/templates/cylinder.py b/pyhope/mesh/transform/templates/cylinder.py index dec11691..5e92c621 100644 --- a/pyhope/mesh/transform/templates/cylinder.py +++ b/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 diff --git a/pyhope/mesh/transform/templates/sphere.py b/pyhope/mesh/transform/templates/sphere.py index 56ec23e4..87014719 100644 --- a/pyhope/mesh/transform/templates/sphere.py +++ b/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: diff --git a/pyhope/mesh/transform/templates/torus.py b/pyhope/mesh/transform/templates/torus.py new file mode 100644 index 00000000..d801294c --- /dev/null +++ b/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 . + +# ================================================================================================================================== +# 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 diff --git a/pyhope/readintools/readintools.py b/pyhope/readintools/readintools.py index 55bb80d9..4f7ec18c 100644 --- a/pyhope/readintools/readintools.py +++ b/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: diff --git a/tutorials/1-06-curved_cylinder/analyze.toml b/tutorials/1-06-curved_cylinder/analyze.toml new file mode 100644 index 00000000..322919ca --- /dev/null +++ b/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 diff --git a/tutorials/1-06-curved_cylinder/parameter.ini b/tutorials/1-06-curved_cylinder/parameter.ini new file mode 100644 index 00000000..c2ab2af3 --- /dev/null +++ b/tutorials/1-06-curved_cylinder/parameter.ini @@ -0,0 +1,70 @@ +DEFVAR=(INT): i0 = 002 ! no. elems in inner square i0xi0 +DEFVAR=(INT): ir = 002 ! no. elems in r +DEFVAR=(INT): iz = 004 ! no. elems in z +DEFVAR=(REAL): ri = 0.5 ! inner square dim +DEFVAR=(REAL): rm = 1. ! middle square dim +DEFVAR=(REAL): r0 = 2. ! outer square dim +DEFVAR=(REAL): lz = 20. ! length of domain in z + +!================================================================================================================================= ! +! OUTPUT +!================================================================================================================================= ! +ProjectName = 1-06-curved_cylinder ! Name of output files +Debugvisu = F ! Launch the GMSH GUI to visualize the mesh +OutputFormat = HDF5 ! Mesh output format (HDF5 VTK) + +!================================================================================================================================= ! +! MESH +!================================================================================================================================= ! +Mode = Internal ! Mode for Cartesian boxes +NGeo = 4 ! Polynomial order of the mapping +nZones = 5 ! Number of boxes + +! center +Corner = (/-ri,-ri,0. ,,ri,-ri,0. ,,ri,ri,0. ,, -ri,ri,0.,, -ri,-ri,lz ,,ri,-ri,lz ,,ri,ri,lz ,, -ri,ri,lz /) +nElems = (/i0,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,0,0,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! left +Corner = (/-r0,-r0,0. ,,-ri,-ri,0. ,,-ri,ri,0. ,, -r0,r0,0.,, -r0,-r0,lz ,,-ri,-ri,lz ,,-ri,ri,lz ,, -r0,r0,lz /) +nElems = (/ir,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,0,5,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! right +Corner = (/ri,-ri,0. ,,r0,-r0,0. ,,r0,r0,0. ,, ri,ri,0.,, ri,-ri,lz ,,r0,-r0,lz ,,r0,r0,lz ,, ri,ri,lz /) +nElems = (/ir,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,3,0,0,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! upper +Corner = (/-ri,ri,0. ,,ri,ri,0. ,,r0,r0,0. ,, -r0,r0,0.,, -ri,ri,lz ,,ri,ri,lz ,,r0,r0,lz ,, -r0,r0,lz /) +nElems = (/i0,ir,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,4,0,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! lower +Corner = (/-r0,-r0,0. ,,r0,-r0,0. ,,ri,-ri,0. ,, -ri,-ri,0.,, -r0,-r0,lz ,,r0,-r0,lz ,,ri,-ri,lz ,, -ri,-ri,lz /) +nElems = (/i0,ir,iz/) ! Number of elements in each direction +BCIndex = (/1,2,0,0,0,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) + +!================================================================================================================================= ! +! BOUNDARY CONDITIONS +!================================================================================================================================= ! +BoundaryName = BC_zminus ! BC index 1 (from position in parameterfile) +BoundaryType = (/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /) +BoundaryName = BC_yminus ! BC index 2 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_xplus ! BC index 3 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_yplus ! BC index 4 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_xminus ! BC index 5 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_zplus ! BC index 6 +BoundaryType = (/1,0,0,-1/) +vv = (/0.,0.,lz/) ! Vector for periodic BC + +!================================================================================================================================= ! +! MESH POST DEFORM +!================================================================================================================================= ! +MeshPostDeform = Cylinder ! Deforms [-1,1]^2 to a cylinder with radius Postdeform_R0 +PostDeform_R0 = 1. ! Domain [-2,2]^2 is mapped to a cylinder with radius 0.5*2 = 1 diff --git a/tutorials/1-07-curved_torus/analyze.toml b/tutorials/1-07-curved_torus/analyze.toml new file mode 100644 index 00000000..624a2d21 --- /dev/null +++ b/tutorials/1-07-curved_torus/analyze.toml @@ -0,0 +1,23 @@ +[ElemInfo] +min = 0.0 +max = 12288.0 +mean = 2280.527777777778 +stddev = 3426.4279625525073 + +[GlobalNodeIDs] +min = 1.0 +max = 5875.0 +mean = 2934.2105305989585 +stddev = 1662.6456197571297 + +[NodeCoords] +min = -6.0 +max = 6.0 +mean = 0.031206087055271378 +stddev = 2.908074017320844 + +[SideInfo] +min = -638.0 +max = 640.0 +mean = 30.497635135135134 +stddev = 169.8852533880788 diff --git a/tutorials/1-07-curved_torus/parameter.ini b/tutorials/1-07-curved_torus/parameter.ini new file mode 100644 index 00000000..98c0d1fa --- /dev/null +++ b/tutorials/1-07-curved_torus/parameter.ini @@ -0,0 +1,71 @@ +DEFVAR=(INT): i0 = 002 ! no. elems in inner square i0xi0 +DEFVAR=(INT): ir = 002 ! no. elems in r +DEFVAR=(INT): iz = 008 ! no. elems in toroidal dir +DEFVAR=(REAL): ri = 0.5 ! inner square dim +DEFVAR=(REAL): rm = 1. ! middle square dim +DEFVAR=(REAL): r0 = 2. ! outer square dim +DEFVAR=(REAL): s0 = 0.5 ! = 1/r0 +DEFVAR=(REAL): rz = 5. ! torus radius + +!================================================================================================================================= ! +! OUTPUT +!================================================================================================================================= ! +ProjectName = 1-07-torus ! Name of output files +Debugvisu = F ! Launch the GMSH GUI to visualize the mesh +OutputFormat = HDF5 ! Mesh output format (HDF5 VTK) + +!================================================================================================================================= ! +! MESH +!================================================================================================================================= ! +Mode = Internal ! Mode for Cartesian boxes +NGeo = 3 ! Polynomial order of the mapping +nZones = 5 ! Number of boxes + +Corner = (/-ri,-ri,0. ,,ri,-ri,0. ,,ri,ri,0. ,, -ri,ri,0.,, -ri,-ri,1. ,,ri,-ri,1. ,,ri,ri,1. ,, -ri,ri,1. /) +nElems = (/i0,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,0,0,3/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! left +Corner = (/-r0,-r0,0. ,,-ri,-ri,0. ,,-ri,ri,0. ,, -r0,r0,0.,, -r0,-r0,1. ,,-ri,-ri,1. ,,-ri,ri,1. ,, -r0,r0,1. /) +nElems = (/ir,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,0,2,3/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! right +Corner = (/ri,-ri,0. ,,r0,-r0,0. ,,r0,r0,0. ,, ri,ri,0.,, ri,-ri,1. ,,r0,-r0,1. ,,r0,r0,1. ,, ri,ri,1. /) +nElems = (/ir,i0,iz/) ! Number of elements in each direction +BCIndex = (/1,0,2,0,0,3/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! upper +Corner = (/-ri,ri,0. ,,ri,ri,0. ,,r0,r0,0. ,, -r0,r0,0.,, -ri,ri,1. ,,ri,ri,1. ,,r0,r0,1. ,, -r0,r0,1. /) +nElems = (/4,ir,iz/) ! Number of elements in each direction +BCIndex = (/1,0,0,2,0,3/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +! lower +Corner = (/-r0,-r0,0. ,,r0,-r0,0. ,,ri,-ri,0. ,, -ri,-ri,0.,, -r0,-r0,1. ,,r0,-r0,1. ,,ri,-ri,1. ,, -ri,-ri,1. /) +nElems = (/i0,ir,iz/) ! Number of elements in each direction +BCIndex = (/1,2,0,0,0,3/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) + +!================================================================================================================================= ! +! BOUNDARY CONDITIONS +!================================================================================================================================= ! +BoundaryName = BC_zminus ! BC index 1 (from position in parameterfile) +BoundaryType = (/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /) +BoundaryName = BC_outer ! BC index 2 +BoundaryType = (/2,1,0,0/) +BoundaryName = BC_zplus ! BC index 3 +BoundaryType = (/1,0,0,-1/) +vv = (/0.,0.,1./) ! Vector for periodic BC + +!================================================================================================================================= ! +! MESH POST DEFORM +!================================================================================================================================= ! +MeshPostDeform = Torus ! Deforms [-1,1]^2 to a cylinder with radius Postdeform_R0 +PostDeform_R0 = s0 ! Domain [-2,2]^2 is mapped to a cylinder with radius 0.5*2 = 1 +PostDeform_Rtorus = rz ! z must be inside [0,1] and periodic + +!================================================================================================================================= ! +! Mesh Checks +!================================================================================================================================= ! +CheckWatertightness = F ! Check if the mesh is watertight + ! > Disabled for torus since trigonometric functions cannot be fully watertight diff --git a/tutorials/1-08-curved_sphere/analyze.toml b/tutorials/1-08-curved_sphere/analyze.toml new file mode 100644 index 00000000..816acdbb --- /dev/null +++ b/tutorials/1-08-curved_sphere/analyze.toml @@ -0,0 +1,23 @@ +[ElemInfo] +min = 0.0 +max = 116000.0 +mean = 20297.330459770114 +stddev = 32963.348084757534 + +[GlobalNodeIDs] +min = 1.0 +max = 60281.0 +mean = 30134.372862068965 +stddev = 17293.308334402507 + +[NodeCoords] +min = -2.0 +max = 2.0 +mean = -0.002304199398700667 +stddev = 0.5604248254465676 + +[SideInfo] +min = -2830.0 +max = 2832.0 +mean = 104.39727011494253 +stddev = 761.4726810886073 diff --git a/tutorials/1-08-curved_sphere/parameter.ini b/tutorials/1-08-curved_sphere/parameter.ini new file mode 100644 index 00000000..3c95b62f --- /dev/null +++ b/tutorials/1-08-curved_sphere/parameter.ini @@ -0,0 +1,129 @@ +DEFVAR=(INT): i0 = 004 ! no. elems in inner cube i0xi0xi0 +DEFVAR=(INT): ir = 003 ! no. elems in r inside +DEFVAR=(INT): jr = 006 ! no. elems in r outside +DEFVAR=(REAL): ri = 0.5 ! inner square dim +DEFVAR=(REAL): rm = 1. ! middle square dim +DEFVAR=(REAL): r0 = 4. ! outer square dim +DEFVAR=(REAL): f1 = 1.2 ! strech factor first ring +DEFVAR=(REAL): f2 = 1.2 ! strech factor second ring + +!================================================================================================================================= ! +! OUTPUT +!================================================================================================================================= ! +ProjectName = 1-08-sphere ! Name of output files +Debugvisu = F ! Launch the GMSH GUI to visualize the mesh +OutputFormat = HDF5 ! Mesh output format (HDF5 VTK) + +!================================================================================================================================= ! +! MESH +!================================================================================================================================= ! +Mode = Internal ! Mode for Cartesian boxes +NGeo = 4 ! Polynomial order of the mapping +nZones = 13 ! Number of boxes + +! center +Corner = (/-ri,-ri,-ri ,,ri,-ri,-ri ,,ri,ri,-ri ,, -ri,ri,-ri,, -ri,-ri,ri ,,ri,-ri,ri ,,ri,ri,ri ,, -ri,ri,ri /) +nElems = (/i0,i0,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,1.,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! left (-x) +Corner = (/-rm,-rm,-rm ,,-ri,-ri,-ri ,,-ri,ri,-ri ,, -rm,rm,-rm,, -rm,-rm,rm ,,-ri,-ri,ri ,,-ri,ri,ri ,, -rm,rm,rm /) +nElems = (/ir,i0,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/-f1,1.,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! right (+x) +Corner = (/ri,-ri,-ri ,,rm,-rm,-rm ,,rm,rm,-rm ,, ri,ri,-ri,, ri,-ri,ri ,,rm,-rm,rm ,,rm,rm,rm ,, ri,ri,ri /) +nElems = (/ir,i0,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/f1,1.,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! lower (-y) +Corner = (/-rm,-rm,-rm ,,rm,-rm,-rm ,,ri,-ri,-ri ,, -ri,-ri,-ri,, -rm,-rm,rm ,,rm,-rm,rm ,,ri,-ri,ri ,, -ri,-ri,ri /) +nElems = (/i0,ir,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,-f1,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! upper (+y) +Corner = (/-ri,ri,-ri ,,ri,ri,-ri ,,rm,rm,-rm ,, -rm,rm,-rm,, -ri,ri,ri ,,ri,ri,ri ,,rm,rm,rm ,, -rm,rm,rm /) +nElems = (/i0,ir,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,f1,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! back (-z) +Corner = (/-rm,-rm,-rm ,,rm,-rm,-rm ,,rm,rm,-rm ,, -rm,rm,-rm,, -ri,-ri,-ri ,,ri,-ri,-ri ,,ri,ri,-ri ,, -ri,ri,-ri/) +nElems = (/i0,i0,ir/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,1.,-f1/) ! Element stretching, either with a constant growth factor (+/- changes direction) +! front (+z) +Corner = (/ -ri,-ri,ri ,,ri,-ri,ri ,,ri,ri,ri ,, -ri,ri,ri,, -rm,-rm,rm ,,rm,-rm,rm ,,rm,rm,rm ,, -rm,rm,rm/) +nElems = (/i0,i0,ir/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,1.,f1/) ! Element stretching, either with a constant growth factor (+/- changes direction) +! left outer (- x) +Corner = (/-r0,-r0,-r0 ,,-rm,-rm,-rm ,,-rm,rm,-rm ,, -r0,r0,-r0,, -r0,-r0,r0 ,,-rm,-rm,rm ,,-rm,rm,rm ,, -r0,r0,r0 /) +nElems = (/jr,i0,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,5,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/-f2,1.,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! right outer ( +x) +Corner = (/rm,-rm,-rm ,,r0,-r0,-r0 ,,r0,r0,-r0 ,, rm,rm,-rm,, rm,-rm,rm ,,r0,-r0,r0 ,,r0,r0,r0 ,, rm,rm,rm /) +nElems = (/jr,i0,i0/) ! Number of elements in each direction +BCIndex = (/0,0,3,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/f2,1.,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! lower outer ( -y) +Corner = (/-r0,-r0,-r0 ,,r0,-r0,-r0 ,,rm,-rm,-rm ,, -rm,-rm,-rm,, -r0,-r0,r0 ,,r0,-r0,r0 ,,rm,-rm,rm ,, -rm,-rm,rm /) +nElems = (/i0,jr,i0/) ! Number of elements in each direction +BCIndex = (/0,2,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,-f2,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! upper outer ( +y) +Corner = (/-rm,rm,-rm ,,rm,rm,-rm ,,r0,r0,-r0 ,, -r0,r0,-r0,, -rm,rm,rm ,,rm,rm,rm ,,r0,r0,r0 ,, -r0,r0,r0 /) +nElems = (/i0,jr,i0/) ! Number of elements in each direction +BCIndex = (/0,0,0,4,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,f2,1./) ! Element stretching, either with a constant growth factor (+/- changes direction) +! back outer (- z) +Corner = (/-r0,-r0,-r0 ,,r0,-r0,-r0 ,,r0,r0,-r0 ,, -r0,r0,-r0,, -rm,-rm,-rm ,,rm,-rm,-rm ,,rm,rm,-rm ,, -rm,rm,-rm/) +nElems = (/i0,i0,jr/) ! Number of elements in each direction +BCIndex = (/1,0,0,0,0,0/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,1.,-f2/) ! Element stretching, either with a constant growth factor (+/- changes direction) +! front outer ( +z) +Corner = (/ -rm,-rm,rm ,,rm,-rm,rm ,,rm,rm,rm ,, -rm,rm,rm,, -r0,-r0,r0 ,,r0,-r0,r0 ,,r0,r0,r0 ,, -r0,r0,r0/) +nElems = (/i0,i0,jr/) ! Number of elements in each direction +BCIndex = (/0,0,0,0,0,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) +Factor = (/1.,1.,f2/) ! Element stretching, either with a constant growth factor (+/- changes direction) + +!================================================================================================================================= ! +! BOUNDARY CONDITIONS +!================================================================================================================================= ! +BoundaryName = BC_zminus ! BC index 1 (from position in parameterfile) +BoundaryType = (/2,0,0,0/) ! (/ Type, curveIndex, State, alpha /) +BoundaryName = BC_yminus ! BC index 2 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_xplus ! BC index 3 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_yplus ! BC index 4 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_xminus ! BC index 5 +BoundaryType = (/2,0,0,0/) +BoundaryName = BC_zplus ! BC index 6 +BoundaryType = (/2,0,0,0/) + +!================================================================================================================================= ! +! MESH POST DEFORM +!================================================================================================================================= ! +MeshPostDeform = sphere +PostDeform_R0 = 0.5 + +!================================================================================================================================= ! +! Mesh Checks +!================================================================================================================================= ! +CheckWatertightness = F ! Check if the mesh is watertight + ! > Disabled for sphere since trigonometric functions cannot be fully watertight diff --git a/tutorials/1-13-curved_periodic_hill/analyze.toml b/tutorials/1-13-curved_periodic_hill/analyze.toml new file mode 100644 index 00000000..e9f4be19 --- /dev/null +++ b/tutorials/1-13-curved_periodic_hill/analyze.toml @@ -0,0 +1,23 @@ +[ElemInfo] +min = 0.0 +max = 131072.0 +mean = 23928.166666666668 +stddev = 36784.532756838016 + +[GlobalNodeIDs] +min = 1.0 +max = 60025.0 +mean = 30014.710418701172 +stddev = 17001.52278169127 + +[NodeCoords] +min = 0.0 +max = 9.0 +mean = 2.7948472471761776 +stddev = 2.208208550209938 + +[SideInfo] +min = -6272.0 +max = 6272.0 +mean = 221.87063802083333 +stddev = 1687.11359440716 diff --git a/tutorials/1-13-curved_periodic_hill/parameter.ini b/tutorials/1-13-curved_periodic_hill/parameter.ini new file mode 100644 index 00000000..bdf2f32b --- /dev/null +++ b/tutorials/1-13-curved_periodic_hill/parameter.ini @@ -0,0 +1,51 @@ + +! Parameters for the Ercoftac periodic hill testcase +! http://www.kbwiki.ercoftac.org/w/index.php/Abstr:2D_Periodic_Hill_Flow +! For Re=700, a resolution of 16x16x8 elems with bell stretching and DXmaxToDXmin = (/1.5,10.,1./) +! yields approximately y+=0.8 and x+=z+=12 at the hill cusp and about half of that everywhere else. + +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = 1-13-curved_periodic_hill ! Name of output files +Debugvisu = F ! Launch the GMSH GUI to visualize the mesh +OutputFormat = HDF5 ! Mesh output format (HDF5 VTK) + +!=============================================================================== ! +! MESH +!=============================================================================== ! +Mode = Internal ! Mode for Cartesian boxes +NGeo = 3 ! Polynomial order of the mapping +nZones = 1 ! Number of boxes + +X0 = (/0.,0. ,0. /) ! Origin of a zone. Equivalent to a corner node. +DX = (/9.,3.035,4.5/) ! Extension of the zone in each spatial direction starting from the origin X0 corner node +nElems = (/16,16 ,8 /) ! Number of elements in each direction +BCIndex = (/1,2,3,4,5,6/) ! Indices of boundary conditions for six boundary faces (z-,y-,x+,y+,x-,z+) +ElemType = 108 ! Element type (108: Hexahedral) + +StretchType = (/3 ,3 ,0 /) ! Stretching type for individual zone per spatial direction. +DXmaxToDXmin = (/1.5,10.,1./) ! Ratio between the smallest and largest element per spatial direction + +!================================================================================================================================= ! +! MESH POST DEFORM +!================================================================================================================================= ! +MeshPostDeform = PHill ! Periodic hill Ercoftac test case + +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! +BoundaryName = BC_periodicz- ! BC index 1 (from position in parameterfile) +BoundaryType = (/1,0,0,2/) ! (/ Type, curveIndex, State, alpha /) ! (/bc-type,curved index,bc state,displacement vector/) +BoundaryName = BC_wall_lower ! BC index 2 +BoundaryType = (/4,0,1,0/) +BoundaryName = BC_periodicx+ ! BC index 3 +BoundaryType = (/1,0,0,-1/) +BoundaryName = BC_wall_upper ! BC index 4 +BoundaryType = (/4,0,1,0/) +BoundaryName = BC_periodicx- ! BC index 5 +BoundaryType = (/1,0,0,1/) +BoundaryName = BC_periodicz+ ! BC index 6 +BoundaryType = (/1,0,0,-2/) ! (/bc-type,curved index,bc state,displacement vector/) +vv = (/9.,0.,0. /) ! Vector for periodic BC +vv = (/0.,0.,4.5/) ! Vector for periodic BC