# -*- 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 <http://www.gnu.org/licenses/>.
"""
from functools import wraps
import numpy as np
from ..view import ViewMesh
from ._convert import (
add_midpoints_edges,
add_midpoints_faces,
add_midpoints_volumes,
cell_types,
collect_edges,
collect_faces,
collect_volumes,
convert,
subdivide,
)
from ._discrete_geometry import DiscreteGeometry
from ._dual import dual
from ._tools import (
expand,
fill_between,
flip,
merge_duplicate_cells,
merge_duplicate_points,
mirror,
revolve,
rotate,
runouts,
translate,
triangulate,
)
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)
[docs]
class Mesh(DiscreteGeometry):
r"""A mesh with points, cells and optional a specified cell type.
Parameters
----------
points : ndarray
Array with point coordinates of shape ``(npoints, dim)``.
cells : ndarray
Array with cell connectivities of shape ``(ncells, points_per_cell)``.
cell_type : str or None, optional
An optional string in VTK-convention that specifies the cell type (default is
None). Necessary when a mesh is saved to a file.
Notes
-----
The ``points``-array is of shape ``(npoints, dim)``, as denoted in Eq. :eq:`points`.
.. math::
:label: points
\hat{\boldsymbol{X}} = \begin{bmatrix}
\hat{\boldsymbol{X}_1}^T \\
\hat{\boldsymbol{X}_2}^T \\
\dots \\
\hat{\boldsymbol{X}_n}^T \\
\end{bmatrix}
Exemplaric point coordinates of a single point at the origin in :math:`\mathbb{R}^3`
are shown in Eq. :eq:`point`.
.. math::
:label: point
\hat{\boldsymbol{X}_1}^T = \begin{bmatrix} 0.0 & 0.0 & 0.0 \end{bmatrix}
Topologic cell connectivies are stored inside the ``cells``-array, see Eq.
:eq:`cells`,
.. math::
:label: cells
\boldsymbol{t} = \begin{bmatrix}
\boldsymbol{t}_1^T \\
\boldsymbol{t}_2^T \\
\vdots \\
\boldsymbol{t}_n^T \\
\end{bmatrix}
with an exemplaric topologic cell-connectivity of a :class:`~felupe.Quad`, as shown
in Eq. :eq:`cell`.
.. math::
:label: cell
\boldsymbol{t}_1^T = \begin{bmatrix} 0 & 1 & 2 & 3 \end{bmatrix}
Examples
--------
.. pyvista-plot::
:context:
:force_static:
>>> import numpy as np
>>> import felupe as fem
>>>
>>> points = np.array(
... [[0.0, 0.0], [0.5, 0.1], [1.0, 0.2], [0.0, 1.0], [0.5, 0.9], [1.0, 0.8]]
... )
>>> cells = np.array([[0, 1, 4, 3], [1, 2, 5, 4]])
>>> mesh = fem.Mesh(points, cells, cell_type="quad")
>>> mesh.plot().show()
A mesh provides several useful attributes, like the dimension, the number of
points, number of cells, number of deegrees of freedom, points connected to cells
and points not connected to cells.
.. pyvista-plot::
:context:
:force_static:
>>> mesh.dim
2
>>> mesh.npoints
6
>>> mesh.ncells
2
>>> mesh.ndof
12
>>> mesh.points_with_cells
array([0, 1, 2, 3, 4, 5])
See Also
--------
felupe.MeshContainer : A container which operates on a list of meshes with identical
dimensions.
felupe.Rectangle : A rectangular 2d-mesh with quads between ``a`` and ``b`` with
``n`` points per axis.
felupe.Cube : A cube shaped 3d-mesh with hexahedrons between ``a`` and ``b`` with
``n`` points per axis.
felupe.Grid : A grid shaped 3d-mesh with hexahedrons. Basically a wrapper for
:func:`numpy.meshgrid` with default ``indexing="ij"``.
felupe.Circle : A circular shaped 2d-mesh with quads and ``n`` points on the
circumferential edge of a 45-degree section. 90-degree ``sections`` are placed
at given angles in degree.
felupe.mesh.Triangle : A triangular shaped 2d-mesh with quads and ``n`` points at
the edges of the three sub-quadrilaterals.
"""
def __init__(self, points, cells, cell_type=None):
self.points = np.array(points)
self.cells = np.array(cells)
self.cell_type = cell_type
super().__init__(points=points, cells=cells, cell_type=cell_type)
self.__mesh__ = Mesh
# alias
self.sweep = self.merge_duplicate_points
self.save = self.write
self.as_pyvista = self.as_unstructured_grid
def __repr__(self):
header = "<felupe Mesh object>"
points = f" Number of points: {len(self.points)}"
cells_header = " Number of cells:"
cells = [f" {self.cell_type}: {self.ncells}"]
return "\n".join([header, points, cells_header, *cells])
def __str__(self):
return self.__repr__()
[docs]
def disconnect(self, points_per_cell=None, calc_points=True):
"""Return a new instance of a Mesh with disconnected cells. Optionally, the
points-per-cell may be specified (must be lower or equal the number of points-
per-cell of the original Mesh). If the Mesh is to be used as a *dual* Mesh, then
the point-coordinates do not have to be re-created because they are not used.
See Also
--------
felupe.Mesh.dual : Create a new dual mesh with given points per cell.
"""
return self.dual(
points_per_cell=points_per_cell,
disconnect=True,
calc_points=calc_points,
offset=0,
npoints=None,
)
[docs]
def as_meshio(self, **kwargs):
"""Export the mesh as :class:`meshio.Mesh`.
Parameters
----------
**kwargs : dict, optional
Additional keyword-arguments for ``meshio.Mesh(points, cells, **kwargs)``.
Returns
-------
meshio.Mesh
The mesh as :class:`meshio.Mesh`.
"""
import meshio
points = np.pad(self.points, ((0, 0), (0, 3 - self.points.shape[1])))
return meshio.Mesh(points=points, cells={self.cell_type: self.cells}, **kwargs)
[docs]
def as_unstructured_grid(self, cell_type=None, **kwargs):
"""Export the mesh as :class:`pyvista.UnstructuredGrid`.
Parameters
----------
cell_type : pyvista.CellType or None, optional
Cell-type of PyVista (default is None).
**kwargs : dict, optional
Additional keyword-arguments for :class:`pyvista.UnstructuredGrid`.
Returns
-------
pyvista.UnstructuredGrid
The mesh as :class:`pyvista.UnstructuredGrid`.
"""
import pyvista as pv
if cell_type is None:
cell_type = dict(cell_types())[self.cell_type]
points = np.pad(self.points, ((0, 0), (0, 3 - self.points.shape[1])))
cells = np.pad(
self.cells, ((0, 0), (1, 0)), constant_values=self.cells.shape[1]
)
cell_types_int = cell_type * np.ones(self.ncells, dtype=int)
return pv.UnstructuredGrid(cells, cell_types_int, points)
[docs]
def write(self, filename="mesh.vtk", **kwargs):
"""Write the mesh to a file.
Parameters
----------
filename : str, optional
The filename of the mesh (default is ``mesh.vtk``).
**kwargs : dict, optional
Additional keyword arguments for :meth:`meshio.Mesh.write`.
Notes
-----
.. note::
For XDMF-export please ensure to have ``h5py`` (as an optional dependency of
``meshio``) installed.
Examples
--------
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=3)
>>> mesh.write(filename="mesh.vtk")
See Also
--------
felupe.mesh.read : Read a mesh from a file using :func:`meshio.read`.
felupe.Mesh.write : Write the mesh to a file.
"""
self.as_meshio(**kwargs).write(filename)
[docs]
def view(self, point_data=None, cell_data=None, cell_type=None):
"""View the mesh with optional given dicts of point- and cell-data items.
Parameters
----------
point_data : dict or None, optional
Additional point-data dict (default is None).
cell_data : dict or None, optional
Additional cell-data dict (default is None).
cell_type : pyvista.CellType or None, optional
Cell-type of PyVista (default is None).
Returns
-------
felupe.ViewMesh
A object which provides visualization methods for :class:`felupe.Mesh`.
See Also
--------
felupe.ViewMesh : Visualization methods for :class:`~felupe.Mesh`.
"""
return ViewMesh(
self,
point_data=point_data,
cell_data=cell_data,
cell_type=cell_type,
)
[docs]
def plot(self, *args, **kwargs):
"""Plot the mesh.
See Also
--------
felupe.Scene.plot : Plot method of a scene.
"""
return self.view().plot(*args, show_undeformed=False, **kwargs)
[docs]
def screenshot(
self,
*args,
filename="mesh.png",
transparent_background=None,
scale=None,
**kwargs,
):
"""Take a screenshot of the mesh.
See Also
--------
pyvista.Plotter.screenshot : Take a screenshot of a PyVista plotter.
"""
return self.plot(*args, off_screen=True, **kwargs).screenshot(
filename=filename,
transparent_background=transparent_background,
scale=scale,
)
[docs]
def imshow(self, *args, ax=None, **kwargs):
"""Take a screenshot of the mesh, show the image data in a figure and return the
ax.
"""
if ax is None:
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.imshow(self.screenshot(*args, filename=None, **kwargs))
ax.set_axis_off()
return ax
[docs]
def add_points(self, points):
"Add additional points to the mesh in-place."
self.update(points=np.vstack([self.points, points]))
[docs]
def clear_points_without_cells(self):
"Clear the list of points without cells."
self.points_without_cells = self.points_without_cells[:0]
[docs]
def get_point_ids(self, value, fun=np.isclose, mode=np.all, **kwargs):
"""Return point ids for points which are close to a given value.
Parameters
----------
value : float, list or ndarray
Scalar value or point coordinates.
fun : callable, optional
The function used to compare the points to the given value, i.e. with a
function signature ``fun(mesh.points, value, **kwargs)``. Default is
:func:`numpy.isclose`.
mode : callable, optional
A callable used to combine the search results, either :func:`numpy.any` or
:func:`numpy.all`.
**kwargs : dict, optional
Additional keyword arguments for ``fun(mesh.points, value, **kwargs)``.
Returns
-------
ndarray
Array with point ids.
Examples
--------
Get point ids at given coordinates for a mesh with duplicate points.
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=11)
>>> mesh.update(points=np.vstack([mesh.points, mesh.points]))
>>> point_ids = mesh.get_point_ids([0, 1, 1])
>>> point_ids
array([1320, 2651])
>>> mesh.points[point_ids]
array([[0., 1., 1.],
[0., 1., 1.]])
"""
return np.argwhere(mode(fun(self.points, value, **kwargs), axis=1))[:, 0]
[docs]
def get_cell_ids(self, point_ids):
"""Return cell ids which have the given point ids in their connectivity.
Parameters
----------
point_ids : list or ndarray
Array with point ids which are used to search for cells.
Returns
-------
ndarray
Array with cell ids which have the given point ids in their connectivity.
Examples
--------
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=11)
>>> point_ids = mesh.get_point_ids([0, 1, 1])
>>> point_ids
array([1320])
>>> cell_ids = mesh.get_cell_ids(point_ids)
>>> cell_ids
array([990])
"""
return np.argwhere(np.isin(self.cells, point_ids).any(axis=1))[:, 0]
[docs]
def get_cell_ids_neighbours(self, cell_ids):
"""Return cell ids which share points with given cell ids.
Parameters
----------
cell_ids : list or ndarray
Array with cell ids which are used to search for neighbour cells.
Returns
-------
ndarray
Array with cell ids which are next to the given cells.
Examples
--------
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=11)
>>> point_ids = mesh.get_point_ids([0, 1, 1])
>>> point_ids
array([1320])
>>> cell_ids = mesh.get_cell_ids(point_ids)
>>> cell_ids
array([990])
Find the cell ids which share at least one point with the given cell id(s).
>>> cell_ids_neighbours = mesh.get_cell_ids_neighbours(cell_ids)
>>> cell_ids_neighbours
array([880, 881, 890, 891, 980, 981, 990, 991])
"""
return self.get_cell_ids(self.cells[cell_ids])
[docs]
def get_point_ids_shared(self, cell_ids_neighbours):
"""Return shared point ids for given cell ids.
Parameters
----------
cells_neighbours : list or ndarray
Array with cell ids.
Returns
-------
ndarray
Array with point ids which are connected to all given cell neighbours.
Examples
--------
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=11)
>>> point_ids = mesh.get_point_ids([0, 1, 1])
>>> point_ids
array([1320])
>>> cell_ids = mesh.get_cell_ids(point_ids)
>>> cell_ids
array([990])
Find the cell ids which share at least one point with the given cell id(s).
>>> cell_ids_neighbours = mesh.get_cell_ids_neighbours(cell_ids)
>>> cell_ids_neighbours
array([880, 881, 890, 891, 980, 981, 990, 991])
Find the shared point ids for the list of cell ids.
>>> point_ids_shared = mesh.get_point_ids_shared(cell_ids_neighbours)
>>> point_ids_shared
array([1189])
"""
neighbours = self.cells[cell_ids_neighbours]
cell = neighbours[0]
return cell[[np.isin(neighbours, point).any(axis=1).all() for point in cell]]
[docs]
def get_point_ids_corners(self):
"""Return point ids which are located at (xmin, ymin), (xmax, ymin), etc.
Returns
-------
ndarray
Array with point ids which are located at the corners.
Examples
--------
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=11)
>>> point_ids_corners = mesh.get_point_ids_corners()
>>> point_ids_corners
array([ 0, 1210, 10, 1220, 110, 1320, 120, 1330])
"""
xmin = np.min(self.points, axis=0)
xmax = np.max(self.points, axis=0)
points = np.vstack(
[x.ravel() for x in np.meshgrid(*np.vstack([xmin, xmax]).T)]
).T
return np.concatenate([self.get_point_ids(point) for point in points])
[docs]
def modify_corners(self, point_ids=None):
"""Modify the corners of a regular rectangle (quad) or cube (hexahedron)
inplace. Only the cells array is modified, the points array remains unchanged.
Parameters
----------
point_ids : ndarray or None, optional
Array with point ids located at the corners which are modified (default is
None). If None, all corners are modified.
Notes
-----
Description of the algorithm:
1. Get corner point ids.
For each corner point:
2. Get attached cell and find cell neighbours.
3. Get the shared point of the cell and its neighbours.
4. Get pair-wise shared points which are located on an edge.
5. Replace the shared points with the corner point.
6. Delete the cell attached to the corner point.
Examples
--------
.. pyvista-plot::
:context:
:force_static:
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(b=(3, 1), n=(16, 6))
>>> mesh.plot().show()
.. pyvista-plot::
:context:
:force_static:
>>> mesh = mesh.modify_corners() # inplace
>>> mesh.plot().show()
"""
if self.cell_type not in ["quad", "hexahedron"]:
message = [
"Cell type not supported.",
"Must be either 'quad' or 'hexahedron'",
f"but given cell type is '{self.cell_type}'.",
]
raise TypeError(" ".join(message))
if point_ids is None:
point_ids = self.get_point_ids_corners()
for point_id in point_ids:
cell_id = self.get_cell_ids(point_id)[0]
cell_id_with_neighbours = self.get_cell_ids_neighbours(cell_id)
cell_ids_neighbours = cell_id_with_neighbours[
cell_id_with_neighbours != cell_id
]
point_id_shared = self.get_point_ids_shared(cell_id_with_neighbours)[0]
point_ids_shared_individual = [
self.get_point_ids_shared([cell_id, neighbour])
for neighbour in cell_ids_neighbours
]
if self.cell_type == "hexahedron":
edges = np.argwhere(
np.isclose(self.points, self.points[point_id]).sum(axis=1) >= 2
)[:, 0]
point_ids_shared_individual = [
p[np.isin(p, edges)] for p in point_ids_shared_individual
]
point_ids_shared_individual = [
shared[shared != point_id_shared]
for shared in point_ids_shared_individual
]
for shared, neighbour in zip(
point_ids_shared_individual, cell_ids_neighbours
):
point_id_replace = np.argwhere(np.isin(self.cells[neighbour], shared))
if len(point_id_replace) > 0:
self.cells[neighbour, point_id_replace[0][0]] = point_id
self.cells = np.delete(self.cells, cell_id, axis=0)
self.update(cells=self.cells)
return self
[docs]
def dual(
self,
points_per_cell=None,
disconnect=True,
calc_points=False,
offset=0,
npoints=None,
):
"""Create a new dual mesh with given points per cell.
Parameters
----------
points_per_cell : int or None, optional
Number of points per cell, must be equal or lower than ``cells.shape[1]`` (
default is None). If None, all points per cell are considered for the dual
mesh.
disconnect : bool, optional
A flag to disconnect the mesh (each cell has its own points). Default is
True.
calc_points : bool, optional
A flag to calculate the point coordinates for the dual mesh (default is
False). If False, the points array is filled with zeros.
offset : int, optional
An offset to be added to the cells array (default is 0).
npoints : int or None, optional
Number of points for the dual mesh. If the given number of points is greater
than ``npoints * points_per_cell``, then the missing points are added to the
points array (filled with zeros). Default is None.
Returns
-------
Mesh
The dual mesh.
Notes
-----
.. note::
The points array of the dual mesh always has a shape of
``(npoints * points_per_cell, dim)``.
Examples
--------
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=5).add_midpoints_edges()
>>> region = fem.RegionQuadraticQuad(mesh=mesh)
>>>
>>> mesh_dual = mesh.dual(points_per_cell=1, disconnect=False)
>>> region_dual = fem.RegionConstantQuad(
... mesh_dual, quadrature=region.quadrature, grad=False
... )
>>>
>>> displacement = fem.FieldPlaneStrain(region, dim=2)
>>> pressure = fem.Field(region_dual)
>>> field = fem.FieldContainer([displacement, pressure])
See Also
--------
felupe.mesh.dual : Create a new dual mesh with given points per cell.
"""
return as_mesh(
dual(
self,
points_per_cell=points_per_cell,
disconnect=disconnect,
calc_points=calc_points,
offset=offset,
npoints=npoints,
)
)
[docs]
def expand(self, n=11, z=1, axis=-1, expand_dim=True):
"""Expand a 0d-Point to a 1d-Line, a 1d-Line to a 2d-Quad or a 2d-Quad to a
3d-Hexahedron Mesh.
Parameters
----------
n : int, optional
Number of n-point repetitions or (n-1)-cell repetitions, default is 11. Must
be greater than 0.
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.
axis : int, optional
Axis of expansion (default is -1).
mask : ndarray or None, optional
A boolean mask to select points which are rotated (default is None).
expand_dim : bool, optional
Expand the dimension of the point coordinates (default is True).
Returns
-------
Mesh
The expanded mesh.
Examples
--------
Expand a rectangle to a cube.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> rect = fem.Rectangle(n=4)
>>> cube = rect.expand(n=7, z=2)
>>>
>>> cube.plot().show()
>>> cube
<felupe Mesh object>
Number of points: 112
Number of cells:
hexahedron: 54
See Also
--------
felupe.mesh.expand : Expand a 0d-Point to a 1d-Line, a 1d-Line to a 2d-Quad or a
2d-Quad to a 3d-Hexahedron Mesh.
"""
return as_mesh(expand(self, n=n, z=z, axis=axis, expand_dim=expand_dim))
[docs]
def rotate(self, angle_deg, axis, center=None, mask=None):
"""Rotate a Mesh.
Parameters
----------
angle_deg : int
Rotation angle in degree.
axis : int
Rotation axis.
center : list or ndarray or None, optional
Center point coordinates (default is None).
mask : ndarray or None, optional
A boolean mask to select points which are rotated (default is None).
Returns
-------
Mesh
The rotated mesh.
Examples
--------
Rotate a rectangle in the xy-plane by 35 degree.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> rect = fem.Rectangle(b=(3, 1), n=(10, 4))
>>> mesh = rect.rotate(angle_deg=35, axis=2, center=[1.5, 0.5])
>>> mesh.plot().show()
>>> mesh
<felupe Mesh object>
Number of points: 40
Number of cells:
quad: 27
See Also
--------
felupe.mesh.rotate : Rotate a Mesh.
"""
return as_mesh(
rotate(self, angle_deg=angle_deg, axis=axis, center=center, mask=mask)
)
[docs]
def revolve(self, n=11, phi=180, axis=0, expand_dim=True):
"""Revolve a 2d-Quad to a 3d-Hexahedron Mesh.
Parameters
----------
n : int, optional
Number of n-point revolutions (or (n-1) cell revolutions),
default is 11.
phi : float or ndarray, optional
Revolution angle in degree (default is 180).
axis : int, optional
Revolution axis (default is 0).
expand_dim : bool, optional
Expand the dimension of the point coordinates (default is True).
Returns
-------
Mesh
The revolved mesh.
Examples
--------
Revolve a cylinder from a rectangle.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> rect = fem.Rectangle(a=(0, 4), b=(3, 5), n=(10, 4))
>>> mesh = rect.revolve(n=11, phi=180, axis=0)
>>> mesh.plot().show()
>>> mesh
<felupe Mesh object>
Number of points: 440
Number of cells:
hexahedron: 270
See Also
--------
felupe.mesh.revolve : Revolve a 2d-Quad to a 3d-Hexahedron Mesh.
"""
return as_mesh(revolve(self, n=n, phi=phi, axis=axis, expand_dim=expand_dim))
[docs]
def merge_duplicate_points(self, decimals=None):
"""Merge duplicate points and update cells of a Mesh.
Parameters
----------
decimals : int or None, optional
Number of decimals for point coordinate comparison (default is None).
Returns
-------
Mesh
The mesh with merged duplicate points and updated cells.
Notes
-----
.. warning::
This function re-sorts points.
.. note::
This function does not merge duplicate cells.
Examples
--------
Two quad meshes to be merged overlap some points. Merge these duplicated
points and update the cells.
.. pyvista-plot::
:context:
>>> import felupe as fem
>>>
>>> rect1 = fem.Rectangle(n=11)
>>> rect2 = fem.Rectangle(a=(0.9, 0), b=(1.9, 1), n=11)
>>> rect2
<felupe Mesh object>
Number of points: 121
Number of cells:
quad: 100
Each mesh contains 121 points and 100 cells. These two meshes are now stored in a
:class:`~felupe.MeshContainer`.
.. pyvista-plot::
:context:
>>> container = fem.MeshContainer([rect1, rect2])
>>> container
<felupe mesh container object>
Number of points: 242
Number of cells:
quad: 100
quad: 100
The meshes of the mesh container are :func:`stacked <felupe.mesh.stack>`.
.. pyvista-plot::
:context:
>>> stack = fem.mesh.stack(container.meshes)
>>> stack
<felupe Mesh object>
Number of points: 242
Number of cells:
quad: 200
After merging the duplicated points and cells, the number of points is reduced
but the number of cells is unchanged.
.. pyvista-plot::
:context:
:force_static:
>>> mesh = stack.merge_duplicate_points()
>>> mesh
<felupe Mesh object>
Number of points: 220
Number of cells:
quad: 200
>>> mesh.plot(opacity=0.6).show()
.. note::
The :class:`~felupe.MeshContainer` may be directly created with
``merge=True``. This enforces :func:`~felupe.mesh.merge_duplicate_points`
for the shared points array of the container.
See Also
--------
felupe.mesh.merge_duplicate_points : Merge duplicated points and update cells of
a Mesh.
felupe.MeshContainer : A container which operates on a list of meshes with
identical dimensions.
"""
return as_mesh(merge_duplicate_points(self, decimals=decimals))
@wraps(merge_duplicate_cells)
def merge_duplicate_cells(self):
"""Merge duplicate cells of a Mesh.
Returns
-------
Mesh
The mesh with merged duplicate cells.
Notes
-----
.. warning::
This function re-sorts cells.
.. note::
This function does not merge duplicate points.
Examples
--------
Two quad meshes to be merged overlap some cells. Merge these duplicated
points and update the cells.
.. pyvista-plot::
:context:
>>> import felupe as fem
>>>
>>> rect1 = fem.Rectangle(n=11)
>>> rect2 = fem.Rectangle(a=(0.9, 0), b=(1.9, 1), n=11)
>>> rect2
<felupe Mesh object>
Number of points: 121
Number of cells:
quad: 100
Each mesh contains 121 points and 100 cells. These two meshes are now stored in
a :class:`~felupe.MeshContainer`.
.. pyvista-plot::
:context:
>>> container = fem.MeshContainer([rect1, rect2])
>>> container
<felupe mesh container object>
Number of points: 242
Number of cells:
quad: 100
quad: 100
The meshes of the mesh container are :func:`stacked <felupe.mesh.stack>`.
.. pyvista-plot::
:context:
>>> stack = fem.mesh.stack(container.meshes)
>>> stack
<felupe Mesh object>
Number of points: 242
Number of cells:
quad: 200
After merging the duplicated points and cells, the number of points is reduced
but the number of cells is unchanged.
.. pyvista-plot::
:context:
:force_static:
>>> mesh = stack.merge_duplicate_points()
>>> mesh
<felupe Mesh object>
Number of points: 220
Number of cells:
quad: 200
mesh.plot(opacity=0.6).show()
.. note::
The :class:`~felupe.MeshContainer` may be directly created with
``merge=True``. This enforces :func:`~felupe.mesh.merge_duplicate_points`
for the shared points array of the container.
The duplicate cells are merged in a second step.
.. pyvista-plot::
:context:
>>> merged = mesh.merge_duplicate_cells()
>>> merged
<felupe Mesh object>
Number of points: 220
Number of cells:
quad: 190
See Also
--------
felupe.mesh.merge_duplicate_points : Merge duplicate points of a Mesh.
felupe.mesh.merge_duplicate_cells : Merge duplicate cells of a Mesh.
felupe.MeshContainer : A container which operates on a list of meshes with
identical dimensions.
"""
return as_mesh(merge_duplicate_cells(self))
[docs]
def fill_between(self, other_mesh, n=11):
"""Fill a 2d-Quad Mesh between two 1d-Line Meshes, embedded in 2d-space, or a
3d-Hexahedron Mesh between two 2d-Quad Meshes, embedded in 3d-space, by expansion.
Both meshes must have equal number of points and cells. The cells-array is taken
from the first mesh.
Parameters
----------
other_mesh : felupe.Mesh
The other line- or quad-mesh.
n : int or ndarray
Number of n-point repetitions or (n-1)-cell repetitions,
(default is 11). If an array is given, then its values are used for the
relative positions in a reference configuration (-1, 1) between the two meshes.
Returns
-------
felupe.Mesh
The expanded mesh.
Examples
--------
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> inner = fem.mesh.revolve(fem.Point(1)).expand(z=0.4).translate(0.2, axis=2)
>>> outer = fem.mesh.revolve(fem.Point(2), phi=160).rotate(
... axis=2, angle_deg=20
... ).expand(z=1.2)
>>> mesh = inner.fill_between(outer, n=6)
>>>
>>> mesh.plot().show()
See Also
--------
felupe.mesh.fill_between : Fill a 2d-Quad Mesh between two 1d-Line Meshes, embedded
in 2d-space, or a 3d-Hexahedron Mesh between two 2d-Quad Meshes, embedded in
3d-space, by expansion.
"""
return as_mesh(fill_between(self, other_mesh=other_mesh, n=n))
[docs]
def flip(self, mask=None):
r"""Ensure positive cell volumes for `tria`, `tetra`, `quad` and `hexahedron`
cell types.
Parameters
----------
mask: list, ndarray or None, optional
Boolean mask for selected cells to flip (default is None). If None, all
cells are selected to be flipped.
Returns
-------
Mesh
The mesh with a rearranged cells array to ensure positive cell volumes.
Examples
--------
A quad mesh with negative cell volumes occurs if one coordinate axis is
multiplied by -1. The error pops up if a region is created with this mesh.
>>> import numpy as np
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=3)
>>> mesh.update(points=mesh.points * np.array([[-1, 1]]))
>>> region = fem.RegionQuad(mesh)
The sum of the differential volumes :math:`V = \sum_c \sum_q dV_{qc}` is
evaluated to -1.0.
>>> region.dV.sum()
-1.0
Let's try to fix the mesh.
>>> mesh.cells
array([[0, 1, 4, 3],
[1, 2, 5, 4],
[3, 4, 7, 6],
[4, 5, 8, 7]])
The cells array is rearranged to ensure positive cell volumes.
>>> mesh_fixed = mesh.flip()
>>> mesh_fixed.cells
array([[3, 4, 1, 0],
[4, 5, 2, 1],
[6, 7, 4, 3],
[7, 8, 5, 4]])
A region now correctly evaluates the total volume of the mesh to 1.0.
>>> region_fixed = fem.RegionQuad(mesh_fixed)
>>> region_fixed.dV.sum()
1.0
See Also
--------
felupe.mesh.flip : Ensure positive cell volumes for `tria`, `tetra`, `quad` and
`hexahedron` cell types.
"""
return as_mesh(flip(self, mask=mask))
[docs]
def mirror(self, normal=[1, 0, 0], centerpoint=[0, 0, 0], axis=None):
"""Mirror points by plane normal and ensure positive cell volumes for
`tria`, `tetra`, `quad` and `hexahedron` cell types.
Parameters
----------
normal: list or ndarray, optional
Mirror-plane normal vector (default is [1, 0, 0]).
centerpoint: list or ndarray, optional
Center-point coordinates on the mirror plane (default is [0, 0, 0]).
axis: int or None, optional
Mirror axis (default is None).
Returns
-------
Mesh
The mirrored mesh.
Examples
--------
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Circle(sections=[0, 90, 180], n=5)
>>> mesh.plot().show()
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Circle(sections=[0, 90, 180], n=5)
>>> mesh.mirror(normal=[0, 1, 0]).plot().show()
See Also
--------
felupe.Mesh.mirror : Mirror points by plane normal and ensure positive cell
volumes for `tria`, `tetra`, `quad` and `hexahedron` cell types.
"""
return as_mesh(mirror(self, normal=normal, centerpoint=centerpoint, axis=axis))
[docs]
def translate(self, move, axis):
"""Translate (move) a Mesh along a given axis.
Parameters
----------
move : float
Translation along given axis.
axis : int
Translation axis.
Returns
-------
Mesh
The translated mesh.
Examples
--------
>>> import felupe as fem
>>>
>>> mesh = fem.Circle(n=6)
>>> mesh.points.min(axis=0), mesh.points.max(axis=0)
(array([-1., -1.]), array([1., 1.]))
>>> translated = mesh.translate(0.3, axis=1)
>>> translated.points.min(axis=0), translated.points.max(axis=0)
(array([-1. , -0.7]), array([1. , 1.3]))
See Also
--------
felupe.mesh.translate : Translate (move) a Mesh along a given axis.
"""
return as_mesh(translate(self, move=move, axis=axis))
[docs]
def triangulate(self, mode=3):
"""Triangulate a quad or a hex mesh.
Parameters
----------
mode: int, optional
Choose a mode how to convert hexahedrons to tets [1]_ (default is 3).
Returns
-------
Mesh
The triangulated mesh.
Examples
--------
Use ``mode=0`` to convert a mesh of hexahedrons into tetrahedrons [1]_.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> mesh.triangulate(mode=0).plot().show()
Use ``mode=3`` to convert a mesh of hexahedrons into tetrahedrons [1]_.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> mesh.triangulate(mode=3).plot().show()
References
----------
.. [1] Dompierre, J., Labbé, P., Vallet, M. G., & Camarero, R. (1999).
How to Subdivide Pyramids, Prisms, and Hexahedra into Tetrahedra.
IMR, 99, 195.
See Also
--------
felupe.mesh.triangulate : Triangulate a quad or a hex mesh.
"""
return as_mesh(triangulate(self, mode=mode))
[docs]
def add_runouts(
self,
values=[0.1, 0.1],
centerpoint=[0, 0, 0],
axis=0,
exponent=5,
mask=slice(None),
normalize=False,
):
"""Add simple rubber-runouts for realistic rubber-metal structures.
Parameters
----------
values : list or ndarray, optional
Relative amount of runouts (per coordinate) perpendicular to the axis
(default is 10% per coordinate, i.e. [0.1, 0.1]).
centerpoint : list or ndarray, optional
Center-point coordinates (default is [0, 0, 0]).
axis : int or None, optional
Axis (default is 0).
exponent : int, optional
Positive exponent to control the shape of the runout. The higher
the exponent, the steeper the transition (default is 5).
mask : list or None, optional
List of points to be considered (default is None).
normalize : bool, optional
Normalize the runouts to create indents, i.e. maintain the original shape at the
ends (default is False).
Returns
-------
Mesh
The mesh with the modified point coordinates.
Examples
--------
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> rect = fem.Rectangle(a=(-3, -1), b=(3, 1), n=(31, 11))
>>> mesh = rect.add_runouts(axis=1, values=[0.2], normalize=True)
>>>
>>> mesh.plot().show()
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> cube = fem.Cube(a=(-3, -2, -1), b=(3, 2, 1), n=(31, 21, 11))
>>> mesh = cube.add_runouts(axis=2, values=[0.1, 0.3], normalize=True)
>>>
>>> mesh.plot().show()
See Also
--------
felupe.mesh.runouts : Add simple rubber-runouts for realistic rubber-metal
structures.
"""
return as_mesh(
runouts(
self,
values=values,
centerpoint=centerpoint,
axis=axis,
exponent=exponent,
mask=mask,
normalize=normalize,
)
)
[docs]
def convert(
self,
order=0,
calc_points=False,
calc_midfaces=False,
calc_midvolumes=False,
):
"""Convert a mesh to a given order. Only conversions to ``order=0`` and
``order=2`` are supported. This function supports meshes with cell types
``"triangle"``, ``"tetra"``, ``"quad"`` and ``"hexahedron"``.
Parameters
----------
order : int, optional
The order of the converted mesh (default is 0). If 0, the points-array will
be of shape ``(ncells, dim)``. If 0 and ``calc_points`` is True, the mean of
all points per cell is evaluated. If 0 and ``calc_points`` is False, the
points array is filled with zeros. If 2, at least midpoints on cell edges
are added to the mesh. If 2 and ``calc_midfaces`` is True, midpoints on cell
faces are also added. If 2 and ``calc_midvolumes`` is True, midpoints on
cell volumes are also added. Raises an error if not 0 or 2.
calc_points : bool, optional
Flag to return the mean of all points per cell if ``order=0`` (default is
False). If False, the points-array is filled with zeros.
calc_midfaces : bool, optional
Flag to add midpoints on cell faces if ``order=2`` (default is False).
calc_midvolumes : bool, optional
Flag to add midpoints on cell volumes if ``order=2`` (default is False).
Returns
-------
Mesh
The converted mesh.
Examples
--------
Convert a mesh of hexahedrons to quadratic hexahedrons by inserting midpoints on
the cell edges.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=6)
>>> mesh2 = mesh.convert(order=2)
>>>
>>> mesh2.plot(plotter=mesh.plot(), style="points", color="black").show()
>>> mesh2
<felupe Mesh object>
Number of points: 96
Number of cells:
quad8: 25
See Also
--------
felupe.mesh.add_midpoints_edges : Add midpoints on edges for given points and cells
and update cell_type accordingly.
felupe.mesh.add_midpoints_faces : Add midpoints on faces for given points and cells
and update cell_type accordingly.
felupe.mesh.add_midpoints_volumes : Add midpoints on volumes for given points and
cells and update cell_type accordingly.
felupe.Mesh.add_midpoints_edges : Add midpoints on edges for given points and cells
and update cell_type accordingly.
felupe.Mesh.add_midpoints_faces : Add midpoints on faces for given points and cells
and update cell_type accordingly.
felupe.Mesh.add_midpoints_volumes : Add midpoints on volumes for given points and
cells and update cell_type accordingly.
"""
return as_mesh(
convert(
self,
order=order,
calc_points=calc_points,
calc_midfaces=calc_midfaces,
calc_midvolumes=calc_midvolumes,
)
)
@wraps(collect_edges)
def collect_edges(self):
return collect_edges(self)
@wraps(collect_faces)
def collect_faces(self):
return collect_faces(self)
@wraps(collect_volumes)
def collect_volumes(self):
return collect_volumes(self)
[docs]
def add_midpoints_edges(self, cell_type=None):
"""Add midpoints on edges for given points and cells and update cell_type
accordingly.
Parameters
----------
cell_type: str or None, optional
A string in VTK-convention that specifies the new cell type (default is None).
If None, the cell type is chosen automatically.
Returns
-------
Mesh
A new mesh with inserted midpoints on cell edges.
Examples
--------
Convert a mesh of hexahedrons to quadratic hexahedrons by inserting midpoints on
the cell edges.
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=6)
>>> mesh_with_midpoints_edges = mesh.add_midpoints_edges()
>>>
>>> mesh_with_midpoints_edges.plot(
... plotter=mesh.plot(), style="points", color="black"
... ).show()
>>> mesh_with_midpoints_edges
<felupe Mesh object>
Number of points: 96
Number of cells:
quad8: 25
See Also
--------
felupe.mesh.add_midpoints_edges : Add midpoints on edges for given points and
cells and update cell_type accordingly.
"""
return add_midpoints_edges(self, cell_type_new=cell_type)
[docs]
def add_midpoints_faces(self, cell_type=None):
"""Add midpoints on faces for given points and cells and update cell_type
accordingly.
Parameters
----------
cell_type: str or None, optional
A string in VTK-convention that specifies the new cell type (default is
None). If None, the cell type is chosen automatically.
Returns
-------
Mesh
A new mesh with inserted midpoints on cell faces.
Examples
--------
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=6)
>>> mesh_with_midpoints_faces = mesh.add_midpoints_faces(cell_type="quad")
>>>
>>> mesh_with_midpoints_faces.plot(
... plotter=mesh.plot(), style="points", color="black"
... ).show()
>>> mesh_with_midpoints_faces
<felupe Mesh object>
Number of points: 61
Number of cells:
quad: 25
See Also
--------
felupe.mesh.add_midpoints_faces : Add midpoints on faces for given points and
cells and update cell_type accordingly.
"""
return add_midpoints_faces(self, cell_type_new=cell_type)
[docs]
def add_midpoints_volumes(self, cell_type=None):
"""Add midpoints on volumes for given points and cells and update cell_type
accordingly.
Parameters
----------
cell_type: str or None, optional
A string in VTK-convention that specifies the new cell type (default is None).
If None, the cell type is chosen automatically.
Returns
-------
Mesh
A new mesh with inserted midpoints on cell volumes.
Examples
--------
.. pyvista-plot::
:force_static:
>>> import felupe as fem
>>> import pyvista as pv
>>>
>>> mesh = fem.Cube(n=6)
>>> mesh_with_midpoints_volumes = fem.mesh.add_midpoints_volumes(
... mesh, cell_type_new="hexahedron9"
... )
>>> plotter = pv.Plotter()
>>> actor = plotter.add_points(mesh_with_midpoints_volumes.points, color="black")
>>> mesh.plot(opacity=0.5, plotter=plotter).show()
>>> mesh_with_midpoints_volumes
<felupe Mesh object>
Number of points: 341
Number of cells:
hexahedron9: 125
See Also
--------
felupe.mesh.add_midpoints_volumes : Add midpoints on volumes for given points
and cells and update cell_type accordingly.
"""
return add_midpoints_volumes(self, cell_type_new=cell_type)
[docs]
def subdivide(self):
"""Subdivide cells by adding midpoints on edges, faces and volumes.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
Take a rectangle mesh and subdivide it two times.
.. pyvista-plot::
:context:
:force_static:
>>> import felupe as fem
>>>
>>> rect = fem.Rectangle(n=4).modify_corners()
>>> rect.plot().show()
.. pyvista-plot::
:context:
:force_static:
>>> mesh = rect.subdivide().subdivide()
>>> mesh.plot().show()
See Also
--------
felupe.Mesh.subdivide : Subdivide cells by adding midpoints on edges, faces and
volumes.
"""
return subdivide(self)
[docs]
def cells_in(self, other_mesh, decimals=8):
"""Return a boolean cells-mask for cells that are also present in another mesh.
Parameters
----------
other_mesh : fem.Mesh
The mesh to compare with.
decimals : int, optional
The number of decimals to consider for point comparison. Default is 8.
Returns
-------
np.ndarray
A boolean array where True indicates that cells are also in the other mesh.
Examples
--------
.. pyvista-plot::
:context:
:force_static:
>>> import felupe as fem
>>>
>>> rect_top = fem.Rectangle(n=(4, 2)).translate(1.0, axis=1)
>>> rect_bottom = fem.Rectangle(n=(4, 2)).translate(1e-6, axis=0)
>>>
>>> meshes = [rect_top, rect_bottom]
>>> mesh = fem.MeshContainer(meshes, merge=True, decimals=6).stack()
>>>
>>> cells_mask = mesh.cells_in(rect_top, decimals=6)
>>> cells_mask
array([ True, True, True, False, False, False])
"""
def view_points(points, decimals):
"Convert point coordinates to integers for comparison."
points = np.round(points * 10**decimals).astype(np.int64)
return points.view([("x", np.int64), ("y", np.int64)]).ravel()
# generate views on points
view = view_points(self.points, decimals)
view_other = view_points(
other_mesh.points[other_mesh.cells].reshape(-1, 2), decimals
)
# shared label mapping
inverse = np.unique(np.concatenate([view, view_other]), return_inverse=True)[1]
labels, local_labels = np.split(inverse, [len(view)])
# point and cell masks
point_mask = np.isin(labels, local_labels)
cell_mask = point_mask[self.cells].all(axis=1)
if cell_mask.sum() < len(other_mesh.cells):
raise ValueError("Cells not found, try to use lower values for decimals.")
if cell_mask.sum() > len(other_mesh.cells):
raise ValueError("More cells found, try to use higher values for decimals.")
return cell_mask
