Merge pull request #524 from adtzlr/add-mesh-circle

Add `Circle(radius, centerpoint, n)` for creating quad-meshed circles

CHANGELOG.md

@@ -7,6 +7,7 @@ All notable changes to this project will be documented in this file. The format
 - Add optional normalization of mesh runouts (which are then indents) by `mesh.runouts(normalize=False)`.
 - Add `LinearElasticLargeStrain(E=None, nu=None, parallel=False)`, suitable for large-rotation analyses. This is based on `NeoHooke()` with converted Lamé-constants.
 - Add a simple boundary-based quad- or hex mesher: A mesh tool for filling the face or volume between two line- or quad-meshes `mesh.fill_between(mesh, other_mesh, n=11)`.
+- Add `Circle(radius, centerpoint, n)` for the creation of a quad-meshed circle.
 
 ### Changed
 - Update the mesh also with a new points array: this changes the Mesh-update function `mesh.update(cells, cell_type=None)` to `mesh.update(points=None, cells=None, cell_type=None)`. Note that this could break old scripts which use `mesh.update(new_cells)` instead of `mesh.update(cells=new_cells)`.

docs/howto/meshgen.rst

@@ -56,8 +56,8 @@ Alternatively, these mesh-related tools are also provided as methods of a :class
     :width: 400px
 
 
-Lines, rectangles and cubes
-***************************
+Lines, rectangles, cubes and circles
+************************************
 
 Of course lines, rectangles, cubes and cylinders do not have to be constructed manually each time. Instead, some easier to use classes are povided by FElupe like :class:`felupe.mesh.Line`, :class:`felupe.Rectangle` or :class:`felupe.Cube`. For non equi-distant points per axis use :class:`felupe.Grid`.
 
@@ -68,6 +68,14 @@ Of course lines, rectangles, cubes and cylinders do not have to be constructed m
 ..  image:: images/cube.png
     :width: 400px
 
+There is also :class:`felupe.Circle` for the creation of a quad-mesh for a circle.
+
+..  code-block:: python
+
+    circle = fem.Circle(radius=1.5, centerpoint=[1, 2], n=6, sections=[0, 90, 180, 270])
+
+..  image:: images/circle.png
+    :width: 400px
 
 Cylinders
 *********
@@ -86,6 +94,7 @@ Cylinders are created by a revolution of a rectangle.
 ..  image:: images/cylinder.png
     :width: 400px
 
+
 Fill between boundaries
 ***********************
 

src/felupe/__init__.py

@@ -82,7 +82,7 @@
     StateNearlyIncompressible,
     Step,
 )
-from .mesh import Cube, Grid, Mesh, MeshContainer, Rectangle
+from .mesh import Circle, Cube, Grid, Mesh, MeshContainer, Rectangle
 from .quadrature import GaussLegendre, GaussLegendreBoundary
 from .quadrature import Tetrahedron as TetrahedronQuadrature
 from .quadrature import Triangle as TriangleQuadrature
@@ -142,6 +142,7 @@
     "FieldsMixed",
     "AreaChange",
     "LinearElastic",
+    "LinearElasticLargeStrain",
     "LinearElasticPlaneStrain",
     "LinearElasticPlaneStress",
     "LinearElasticPlasticIsotropicHardening",
@@ -184,6 +185,7 @@
     "Triangle",
     "TriangleMINI",
     "TriQuadraticHexahedron",
+    "Circle",
     "Cube",
     "Grid",
     "Mesh",

src/felupe/mesh/__init__.py

@@ -10,6 +10,7 @@
 )
 from ._dual import dual
 from ._geometry import (
+    Circle,
     Cube,
     CubeArbitraryOrderHexahedron,
     Grid,
@@ -50,6 +51,7 @@
     "collect_volumes",
     "convert",
     "dual",
+    "Circle",
     "Cube",
     "CubeArbitraryOrderHexahedron",
     "Grid",
@@ -61,6 +63,7 @@
     "concatenate",
     "expand",
     "flip",
+    "fill_between",
     "mirror",
     "revolve",
     "rotate",

src/felupe/mesh/_geometry.py

@@ -17,11 +17,10 @@
 """
 
 import numpy as np
-from scipy.interpolate import griddata
 
-from ..math import transpose
 from ._line_rectangle_cube import cube_hexa, line_line, rectangle_quad
 from ._mesh import Mesh
+from ._tools import concatenate
 
 
 class Line(Mesh):
@@ -216,3 +215,40 @@ def search_face(p, a, x, vertices, edges):
         cells = np.hstack((vertices, edges, faces, volume)).reshape(1, -1)
 
         super().__init__(points, cells, cell_type="VTK_LAGRANGE_HEXAHEDRON")
+
+
+class Circle(Mesh):
+    def __init__(
+        self,
+        radius=1,
+        centerpoint=[0, 0],
+        n=6,
+        sections=[0, 90, 180, 270],
+        value=0.15,
+        exponent=2,
+        decimals=10,
+    ):
+        rect = Rectangle(b=(0.5, 0.5), n=(n, n))
+        right = rect.points[:, 0] == 0.5
+
+        rect = rect.add_runouts(values=[-value], axis=0, exponent=exponent)
+        rect.points[:, 0] = rect.points[:, 1].reshape(n, n, order="F").ravel()
+
+        line = Line(n=n)
+        phi = np.linspace(np.pi / 4, 0, n)
+
+        bottom = line.copy(points=rect.points[right][::-1])
+        top = line.copy(points=np.vstack([np.cos(phi), np.sin(phi)]).T)
+        face = bottom.fill_between(top, n=n)
+
+        quarter = concatenate([face, face.mirror(normal=[-1, 1]), rect])
+        circle = concatenate([quarter.rotate(alpha, 2) for alpha in sections]).sweep(
+            decimals=decimals
+        )
+
+        circle.points *= radius
+        circle.points += np.array(centerpoint)
+
+        super().__init__(
+            points=circle.points, cells=circle.cells, cell_type=circle.cell_type
+        )

src/felupe/mesh/_mesh.py

@@ -16,7 +16,6 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-from copy import deepcopy
 from functools import wraps
 
 import numpy as np

tests/test_mesh.py

@@ -410,6 +410,10 @@ def test_mesh_fill_between():
 
     assert np.allclose(face1.points, face2.points)
 
+    region = fe.RegionQuad(face1)
+
+    assert np.all(region.dV > 0)
+
     # 3d
     face = fe.Rectangle(n=(6, 7))
 
@@ -425,6 +429,40 @@ def test_mesh_fill_between():
 
     assert np.allclose(solid1.points, solid2.points)
 
+    region = fe.RegionHexahedron(solid1)
+
+    assert np.all(region.dV > 0)
+
+
+def test_circle():
+    centerpoint = [0, 0]
+    radius = 1.5
+    x0, y0 = centerpoint
+
+    mesh = fe.Circle(
+        radius=radius,
+        centerpoint=centerpoint,
+        n=6,
+        sections=[0, 90, 180],
+        exponent=2,
+        value=0.15,
+        decimals=12,
+    )
+
+    x, y = mesh.points.T
+
+    r = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
+
+    assert np.all(r - 1e-10 <= 1.5)
+
+    region = fe.RegionQuad(mesh)
+
+    assert np.all(region.dV > 0)
+
+    boundary = fe.RegionQuadBoundary(mesh)
+
+    assert boundary.mesh.ncells == 50
+
 
 if __name__ == "__main__":
     test_meshes()
@@ -439,3 +477,4 @@ def test_mesh_fill_between():
     test_mesh_methods()
     test_read_nocells(filename="mesh_no-cells.bdf")
     test_mesh_fill_between()
+    test_circle()