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")