diff --git a/CHANGELOG.md b/CHANGELOG.md
index ee756524..26f165b8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,8 @@ All notable changes to this project will be documented in this file. The format
### Changed
- Change `einsumt` from an optional to a required dependency.
- Vectorize implementations of `MultiPointConstraint` and `MultiPointContact` and re-implement both as `scipy.sparse.lil_matrix()`.
+- Rename old `Mesh` to `DiscreteGeometry` and rebase new `Mesh` on `DiscreteGeometry`.
+- Simplify the usage of explicit mesh-related tools by adding them as methods to `Mesh`, i.e. `mesh.tools.expand(Rectangle())` is equivalent to `Rectangle().expand()`.
### Fixed
- Fix `tools.project()` for higher-order quad- and hexahedron elements.
diff --git a/src/felupe/mesh/_convert.py b/src/felupe/mesh/_convert.py
index 7c52c96e..77530790 100644
--- a/src/felupe/mesh/_convert.py
+++ b/src/felupe/mesh/_convert.py
@@ -27,8 +27,8 @@
import numpy as np
+from ._discrete_geometry import DiscreteGeometry
from ._helpers import mesh_or_data
-from ._mesh import Mesh
@mesh_or_data
diff --git a/src/felupe/mesh/_discrete_geometry.py b/src/felupe/mesh/_discrete_geometry.py
new file mode 100644
index 00000000..26511465
--- /dev/null
+++ b/src/felupe/mesh/_discrete_geometry.py
@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+"""
+ _______ _______ ___ __ __ _______ _______
+| || || | | | | || || |
+| ___|| ___|| | | | | || _ || ___|
+| |___ | |___ | | | |_| || |_| || |___
+| ___|| ___|| |___ | || ___|| ___|
+| | | |___ | || || | | |___
+|___| |_______||_______||_______||___| |_______|
+
+This file is part of felupe.
+
+Felupe 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.
+
+Felupe 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 Felupe. If not, see .
+
+"""
+
+import numpy as np
+
+
+class DiscreteGeometry:
+ """A discrete geometry with points, cells and optional a specified cell type.
+
+ Parameters
+ ----------
+ points : ndarray
+ Point coordinates.
+ cells : ndarray
+ Point-connectivity of cells.
+ cell_type : str or None, optional
+ An optional string in VTK-convention that specifies the cell type (default is None). Necessary when a discrete geometry is saved to a file.
+
+ Attributes
+ ----------
+ points : ndarray
+ Point coordinates.
+ cells : ndarray
+ Point-connectivity of cells.
+ cell_type : str or None
+ A string in VTK-convention that specifies the cell type.
+ npoints : int
+ Amount of points.
+ dim : int
+ Dimension of mesh point coordinates.
+ ndof : int
+ Amount of degrees of freedom.
+ ncells : int
+ Amount of cells.
+ points_with_cells : array
+ Array with points connected to cells.
+ points_without_cells : array
+ Array with points not connected to cells.
+ cells_per_point : array
+ Array which counts connected cells per point. Used for averaging results.
+
+ """
+
+ def __init__(self, points, cells, cell_type=None):
+ self.points = np.array(points)
+ self.cells = np.array(cells)
+ self.cell_type = cell_type
+
+ self.update(self.cells)
+
+ def update(self, cells, cell_type=None):
+ "Update the cell and dimension attributes with a given cell array."
+ self.cells = cells
+
+ if cell_type is not None:
+ self.cell_type = cell_type
+
+ # obtain dimensions
+ self.npoints, self.dim = self.points.shape
+ self.ndof = self.points.size
+ self.ncells = self.cells.shape[0]
+
+ # get number of cells per point
+ points_in_cell, self.cells_per_point = np.unique(cells, return_counts=True)
+
+ # check if there are points without cells
+ if self.npoints != len(self.cells_per_point):
+ self.point_has_cell = np.isin(np.arange(self.npoints), points_in_cell)
+ # update "cells_per_point" ... cells per point
+ cells_per_point = -np.ones(self.npoints, dtype=int)
+ cells_per_point[points_in_cell] = self.cells_per_point
+ self.cells_per_point = cells_per_point
+
+ self.points_without_cells = np.arange(self.npoints)[~self.point_has_cell]
+ self.points_with_cells = np.arange(self.npoints)[self.point_has_cell]
+ else:
+ self.points_without_cells = np.array([], dtype=int)
+ self.points_with_cells = np.arange(self.npoints)
diff --git a/src/felupe/mesh/_helpers.py b/src/felupe/mesh/_helpers.py
index 1934bd65..22aee006 100644
--- a/src/felupe/mesh/_helpers.py
+++ b/src/felupe/mesh/_helpers.py
@@ -27,12 +27,10 @@
from functools import wraps
-from ._mesh import Mesh
-
def mesh_or_data(meshfun):
- """If a ``Mesh`` is passed to a mesh function, extract ``points`` and
- ``cells`` arrays along with the ``cell_type`` and return a ``Mesh``
+ """If a ``DiscreteGeometry`` is passed to a mesh function, extract ``points`` and
+ ``cells`` arrays along with the ``cell_type`` and return a ``DiscreteGeometry``
as a result."""
@wraps(meshfun)
@@ -44,8 +42,8 @@ def mesh_or_points_cells_type(*args, **kwargs):
# check if unnamed args are passed
if len(args) > 0:
- # meshfun(Mesh)
- if isinstance(args[0], Mesh):
+ # meshfun(DiscreteGeometry)
+ if hasattr(args[0], "__mesh__"):
# set mesh flag
is_mesh = True
@@ -55,6 +53,9 @@ def mesh_or_points_cells_type(*args, **kwargs):
cells = args[0].cells
cell_type = args[0].cell_type
+ # get mesh class
+ Mesh = args[0].__mesh__
+
# remove Mesh from args
args = args[1:]
@@ -104,7 +105,7 @@ def mesh_or_points_cells_type(*args, **kwargs):
# call mesh manipulation function
points, cells, cell_type = meshfun(points, cells, cell_type, *args, **kwargs)
- # return a Mesh if a Mesh was passed
+ # return a DiscreteGeometry if a DiscreteGeometry was passed
if is_mesh:
return Mesh(points=points, cells=cells, cell_type=cell_type)
diff --git a/src/felupe/mesh/_mesh.py b/src/felupe/mesh/_mesh.py
index 8d3e6b79..fe389e79 100644
--- a/src/felupe/mesh/_mesh.py
+++ b/src/felupe/mesh/_mesh.py
@@ -26,11 +26,29 @@
"""
from copy import deepcopy
+from functools import wraps
import numpy as np
+from ._convert import (
+ add_midpoints_edges,
+ add_midpoints_faces,
+ add_midpoints_volumes,
+ collect_edges,
+ collect_faces,
+ collect_volumes,
+ convert,
+)
+from ._discrete_geometry import DiscreteGeometry
+from ._tools import expand, mirror, revolve, rotate, runouts, sweep, triangulate
-class Mesh:
+
+def as_mesh(obj):
+ "Convert a ``DiscreteGeometry`` object to a ``Mesh`` object."
+ return Mesh(points=obj.points, cells=obj.cells, cell_type=obj.cell_type)
+
+
+class Mesh(DiscreteGeometry):
"""A mesh with points, cells and optional a specified cell type.
Parameters
@@ -50,20 +68,6 @@ class Mesh:
Point-connectivity of cells.
cell_type : str or None
A string in VTK-convention that specifies the cell type.
- npoints : int
- Amount of points.
- dim : int
- Dimension of mesh point coordinates.
- ndof : int
- Amount of degrees of freedom.
- ncells : int
- Amount of cells.
- points_with_cells : array
- Array with points connected to cells.
- points_without_cells : array
- Array with points not connected to cells.
- cells_per_point : array
- Array which counts connected cells per point. Used for averging results.
"""
@@ -72,36 +76,20 @@ def __init__(self, points, cells, cell_type=None):
self.cells = np.array(cells)
self.cell_type = cell_type
- self.update(self.cells)
-
- def update(self, cells, cell_type=None):
- "Update the cell and dimension attributes with a given cell array."
- self.cells = cells
+ super().__init__(points=points, cells=cells, cell_type=cell_type)
- if cell_type is not None:
- self.cell_type = cell_type
+ self.__mesh__ = Mesh
- # obtain dimensions
- self.npoints, self.dim = self.points.shape
- self.ndof = self.points.size
- self.ncells = self.cells.shape[0]
-
- # get number of cells per point
- points_in_cell, self.cells_per_point = np.unique(cells, return_counts=True)
+ def __repr__(self):
+ header = ""
+ points = f" Number of points: {len(self.points)}"
+ cells_header = " Number of cells:"
+ cells = [f" {self.cell_type}: {self.ncells}"]
- # check if there are points without cells
- if self.npoints != len(self.cells_per_point):
- self.point_has_cell = np.isin(np.arange(self.npoints), points_in_cell)
- # update "cells_per_point" ... cells per point
- cells_per_point = -np.ones(self.npoints, dtype=int)
- cells_per_point[points_in_cell] = self.cells_per_point
- self.cells_per_point = cells_per_point
+ return "\n".join([header, points, cells_header, *cells])
- self.points_without_cells = np.arange(self.npoints)[~self.point_has_cell]
- self.points_with_cells = np.arange(self.npoints)[self.point_has_cell]
- else:
- self.points_without_cells = np.array([], dtype=int)
- self.points_with_cells = np.arange(self.npoints)
+ def __str__(self):
+ return self.__repr__()
def disconnect(self, points_per_cell=None, calc_points=True):
"""Return a new instance of a Mesh with disconnected cells. Optionally, the
@@ -160,13 +148,88 @@ def copy(self):
return deepcopy(self)
- def __repr__(self):
- header = ""
- points = f" Number of points: {len(self.points)}"
- cells_header = " Number of cells:"
- cells = [f" {self.cell_type}: {self.ncells}"]
+ @wraps(expand)
+ def expand(self, n=11, z=1):
+ return as_mesh(expand(self, n=n, z=z))
+
+ @wraps(rotate)
+ def rotate(self, angle_deg, axis, center=None):
+ return as_mesh(rotate(self, angle_deg=angle_deg, axis=axis, center=center))
+
+ @wraps(revolve)
+ def revolve(self, n=11, phi=180, axis=0):
+ return as_mesh(revolve(self, n=n, phi=phi, axis=axis))
+
+ @wraps(sweep)
+ def sweep(self, decimals=None):
+ return as_mesh(sweep(self, decimals=decimals))
+
+ @wraps(mirror)
+ def mirror(self, normal=[1, 0, 0], centerpoint=[0, 0, 0], axis=None):
+ return as_mesh(mirror(self, normal=normal, centerpoint=centerpoint, axis=axis))
+
+ @wraps(triangulate)
+ def triangulate(self, mode=3):
+ return as_mesh(triangulate(self, mode=mode))
+
+ @wraps(runouts)
+ def add_runouts(
+ self,
+ values=[0.1, 0.1],
+ centerpoint=[0, 0, 0],
+ axis=0,
+ exponent=5,
+ mask=slice(None),
+ ):
+ return as_mesh(
+ runouts(
+ self,
+ values=values,
+ centerpoint=centerpoint,
+ axis=axis,
+ exponent=exponent,
+ mask=mask,
+ )
+ )
+
+ @wraps(convert)
+ def convert(
+ self,
+ order=0,
+ calc_points=False,
+ calc_midfaces=False,
+ calc_midvolumes=False,
+ ):
+ return as_mesh(
+ convert(
+ self,
+ order=order,
+ calc_points=calc_points,
+ calc_midfaces=calc_midfaces,
+ calc_midvolumes=calc_midvolumes,
+ )
+ )
- return "\n".join([header, points, cells_header, *cells])
+ @wraps(collect_edges)
+ def collect_edges(self):
+ return collect_edges
- def __str__(self):
- return self.__repr__()
+ @wraps(collect_faces)
+ def collect_faces(self):
+ return collect_faces
+
+ @wraps(collect_volumes)
+ def collect_volumes(self):
+ return collect_volumes
+
+ @wraps(add_midpoints_edges)
+ def add_midpoints_edges(self):
+ return add_midpoints_edges
+
+ @wraps(add_midpoints_faces)
+ def add_midpoints_faces(self):
+ return add_midpoints_faces
+
+ @wraps(add_midpoints_volumes)
+ def add_midpoints_volumes(self):
+ return add_midpoints_volumes
diff --git a/src/felupe/mesh/_tools.py b/src/felupe/mesh/_tools.py
index f1561aaf..11440fbf 100644
--- a/src/felupe/mesh/_tools.py
+++ b/src/felupe/mesh/_tools.py
@@ -29,7 +29,6 @@
from ..math import rotation_matrix
from ._helpers import mesh_or_data
-from ._mesh import Mesh
@mesh_or_data
@@ -332,6 +331,8 @@ def mirror(
def concatenate(meshes):
"Join a sequence of meshes with identical cell types."
+ Mesh = meshes[0].__mesh__
+
points = np.vstack([mesh.points for mesh in meshes])
offsets = np.cumsum(np.insert([mesh.npoints for mesh in meshes][:-1], 0, 0))
cells = np.vstack([offset + mesh.cells for offset, mesh in zip(offsets, meshes)])
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index ce35604d..73abf195 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -52,6 +52,7 @@ def test_meshes():
fe.mesh.convert(m, order=0, calc_points=True)
fe.mesh.convert(m, order=2)
fe.mesh.convert(m, order=2, calc_midfaces=True)
+ m.convert(order=2)
m = fe.Mesh(
points=np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]),
@@ -69,6 +70,7 @@ def test_meshes():
assert m.cells.shape == (4, 2)
mr = fe.mesh.revolve(m, n=11, phi=180, axis=2)
+ mr = m.revolve(n=11, phi=180, axis=2)
assert mr.ncells == 4 * 10
m = fe.Rectangle(a=(-1.2, -2), b=(2, 3.1), n=(4, 9))
@@ -105,6 +107,7 @@ def test_meshes():
fe.mesh.expand(m.points, cells=m.cells, cell_type=m.cell_type)
fe.mesh.expand(points=m.points, cells=m.cells, cell_type=m.cell_type)
fe.mesh.expand(m)
+ m.expand()
me1 = fe.mesh.expand(m, n=3, z=1)
me2 = fe.mesh.expand(m, n=3, z=1.0)
@@ -130,6 +133,7 @@ def test_meshes():
fe.mesh.rotate(m, angle_deg=10, axis=0, center=None)
fe.mesh.rotate(m.points, m.cells, m.cell_type, angle_deg=10, axis=0, center=None)
fe.mesh.rotate(m, angle_deg=10, axis=1, center=[0, 0, 0])
+ m.rotate(angle_deg=10, axis=0, center=None)
fe.mesh.CubeArbitraryOrderHexahedron()
fe.mesh.RectangleArbitraryOrderQuad()
@@ -150,6 +154,7 @@ def test_meshes():
fe.mesh.sweep(m)
fe.mesh.sweep(m.points, m.cells, m.cell_type, decimals=4)
+ m.sweep()
m.as_meshio(point_data={"data": m.points}, cell_data={"cell_data": [m.cells[:, 0]]})
m.save()
@@ -176,6 +181,7 @@ def test_mirror():
m = fe.mesh.Line()
r = fe.Region(m, fe.Line(), fe.GaussLegendre(1, 1))
n = fe.mesh.mirror(m, **kwargs)
+ n = m.mirror(**kwargs)
s = fe.Region(n, fe.Line(), fe.GaussLegendre(1, 1))
assert np.isclose(r.dV.sum(), s.dV.sum())
@@ -222,6 +228,7 @@ def test_mirror():
def test_triangulate():
m = fe.Rectangle(n=3)
n = fe.mesh.triangulate(m)
+ n = m.triangulate()
rm = fe.RegionQuad(m)
rn = fe.RegionTriangle(n)
@@ -244,6 +251,7 @@ def test_triangulate():
def test_runouts():
m = fe.Rectangle(n=3)
+ n = m.add_runouts(values=[0.0], axis=0, centerpoint=[0, 0])
n = fe.mesh.runouts(m, values=[0.0], axis=0, centerpoint=[0, 0])
assert n.points[:, 1].max() == m.points[:, 1].max()
@@ -346,6 +354,16 @@ def test_read_nocells(filename="tests/mesh_no-cells.bdf"):
assert mesh[0].cells.shape == (0, 0)
+def test_mesh_methods():
+ mesh = fe.Cube()
+ mesh.collect_edges()
+ mesh.collect_faces()
+ mesh.collect_volumes()
+ mesh.add_midpoints_edges()
+ mesh.add_midpoints_faces()
+ mesh.add_midpoints_volumes()
+
+
if __name__ == "__main__":
test_meshes()
test_mirror()
@@ -356,4 +374,5 @@ def test_read_nocells(filename="tests/mesh_no-cells.bdf"):
test_grid_1d()
test_container()
test_read(filename="mesh.bdf")
+ test_mesh_methods()
test_read_nocells(filename="mesh_no-cells.bdf")