Improve Mesh

README.md

@@ -58,13 +58,16 @@ All notable changes to this project will be documented in this file. The format
 - Add `RegionBoundary` along with template regions for `Quad` and `Hexahedron` and `GaussLegendreBoundary`.
 - Add optional normal vector argument for function and gradient methods of `AreaChange`.
 - Add a new Mesh-tool `triangulate()`, applicable on Quad and Hexahedron meshes.
+- Add a new Mesh-method `Mesh.as_meshio()`.
 
 ### Changed
 - Enforce consistent arguments for functions inside `mesh` (`points, cells, cell_data` or `Mesh`).
 - Rename Numba-`parallel` assembly to `jit`.
 - Move single element shape functions and their derivatives from `region.h` to `region.element.h` and `region.dhdr` to `region.element.dhdr`.
 - [Repeat](https://numpy.org/doc/stable/reference/generated/numpy.tile.html) element shape functions and their derivatives for each cell (as preparation for an upcoming `RegionBoundary`).
 - Improve `mesh.convert()` by using the function decorator `@mesh_or_data`.
+- Allow an array to be passed as the expansion arguments of `mesh.expand()` and `mesh.revolve()`.
+- Allow optional keyword args to be passed to `Mesh.save(**kwargs)`, acts as a wrapper for `Mesh.as_meshio(**kwargs).write()`.
 
 ### Fixed
 - Fix area normal vectors of `RegionBoundary`.

felupe/mesh/_mesh.py

@@ -107,8 +107,8 @@ def disconnect(self):
 
         return Mesh(points, cells, cell_type=self.cell_type)
 
-    def save(self, filename="mesh.vtk"):
-        """Export the mesh as VTK file. For XDMF-export please ensure to have ``h5py`` (as an optional dependancy of ``meshio``) installed.
+    def as_meshio(self, **kwargs):
+        """Export the mesh as ``meshio.Mesh``.
 
         Parameters
         ----------
@@ -123,7 +123,20 @@ def save(self, filename="mesh.vtk"):
             import meshio
 
         cells = {self.cell_type: self.cells}
-        meshio.Mesh(self.points, cells).write(filename)
+        return meshio.Mesh(self.points, cells, **kwargs)
+
+    def save(self, filename="mesh.vtk", **kwargs):
+        """Export the mesh as VTK file. For XDMF-export please ensure to have
+        ``h5py`` (as an optional dependancy of ``meshio``) installed.
+
+        Parameters
+        ----------
+        filename : str, optional
+            The filename of the mesh (default is ``mesh.vtk``).
+
+        """
+
+        self.as_meshio(**kwargs).write(filename)
 
     def copy(self):
         """Return a deepcopy of the mesh.

felupe/mesh/_tools.py

@@ -46,9 +46,10 @@ def expand(points, cells, cell_type, n=11, z=1):
     n : int, optional
         Number of n-point repetitions or (n-1)-cell repetitions,
         default is 11.
-    z : int, optional
-        Total expand dimension (edge length in expand direction is z / n),
-        default is 1.
+    z : float or ndarray, optional
+        Total expand dimension as float (edge length in expand direction is z / n),
+        default is 1. Optionally, if an array is passed these entries are
+        taken as expansion and `n` is ignored.
 
     Returns
     -------
@@ -79,7 +80,13 @@ def expand(points, cells, cell_type, n=11, z=1):
     p = np.pad(points, ((0, 0), (0, 1)))
 
     # generate new points array for every thickness expansion ``h``
-    points_new = np.vstack([p + np.array([*zeros, h]) for h in np.linspace(0, z, n)])
+    if isinstance(z, int) or isinstance(z, float):
+        points_z = np.linspace(0, z, n)
+    else:
+        points_z = z
+        n = len(z)
+
+    points_new = np.vstack([p + np.array([*zeros, h]) for h in points_z])
 
     # generate new cells array
     c = [cells + len(p) * a for a in np.arange(n)]
@@ -148,7 +155,7 @@ def revolve(points, cells, cell_type, n=11, phi=180, axis=0):
     n : int, optional
         Number of n-point revolutions (or (n-1) cell revolutions),
         default is 11.
-    phi : int, optional
+    phi : float or ndarray, optional
         Revolution angle in degree (default is 180).
     axis : int, optional
         Revolution axis (default is 0).
@@ -174,19 +181,23 @@ def revolve(points, cells, cell_type, n=11, phi=180, axis=0):
         "quad": ("hexahedron", slice(None, None, None)),
     }[cell_type]
 
-    if abs(phi) > 360:
+    if isinstance(phi, int) or isinstance(phi, float):
+        points_phi = np.linspace(0, phi, n)
+    else:
+        points_phi = phi
+        n = len(points_phi)
+
+    if abs(points_phi[-1]) > 360:
         raise ValueError("phi must be within |phi| <= 360 degree.")
 
     p = np.pad(points, ((0, 0), (0, 1)))
     R = rotation_matrix
 
-    points_new = np.vstack(
-        [(R(angle, dim + 1) @ p.T).T for angle in np.linspace(0, phi, n)]
-    )
+    points_new = np.vstack([(R(angle, dim + 1) @ p.T).T for angle in points_phi])
 
     c = [cells + len(p) * a for a in np.arange(n)]
 
-    if phi == 360:
+    if points_phi[-1] == 360:
         c[-1] = c[0]
         points_new = points_new[: len(points_new) - len(points)]
 

tests/test_mesh.py

@@ -81,10 +81,14 @@ def test_meshes():
     with pytest.raises(NotImplementedError):
         fe.mesh.convert(m, order=1)
 
-    fe.mesh.revolve(m, n=11, phi=180, axis=0)
+    mr1 = fe.mesh.revolve(m, n=11, phi=180, axis=0)
     fe.mesh.revolve(m.points, m.cells, m.cell_type, n=11, phi=180, axis=0)
     fe.mesh.revolve(m.points, m.cells, m.cell_type, n=11, phi=360, axis=0)
 
+    mr2 = fe.mesh.revolve(m, phi=np.linspace(0, 180, 11), axis=0)
+
+    assert np.allclose(mr1.points, mr2.points)
+
     with pytest.raises(ValueError):
         fe.mesh.revolve(m.points, m.cells, m.cell_type, n=11, phi=361, axis=0)
 
@@ -94,6 +98,13 @@ def test_meshes():
     fe.mesh.expand(points=m.points, cells=m.cells, cell_type=m.cell_type)
     fe.mesh.expand(m)
 
+    me1 = fe.mesh.expand(m, n=3, z=1)
+    me2 = fe.mesh.expand(m, n=3, z=1.0)
+    me3 = fe.mesh.expand(m, z=np.linspace(0, 1, 3))
+
+    assert np.allclose(me1.points, me2.points)
+    assert np.allclose(me2.points, me3.points)
+
     m = fe.Cube(a=(-1, -2, -0.5), b=(2, 3.1, 1), n=(4, 9, 5))
     assert m.points.shape == (4 * 9 * 5, 3)
     assert m.cells.shape == (3 * 8 * 4, 8)
@@ -127,6 +138,7 @@ def test_meshes():
     fe.mesh.sweep(m)
     fe.mesh.sweep(m.points, m.cells, m.cell_type, decimals=4)
 
+    m.as_meshio(point_data={"data": m.points}, cell_data={"cell_data": [m.cells[:, 0]]})
     m.save()
 
     m.cell_type = None