Skip to content

Commit

Permalink
Merge pull request #425 from adtzlr/mesh-tools-as-methods
Browse files Browse the repository at this point in the history
Enhance `Mesh`: Mesh tools as methods
  • Loading branch information
adtzlr authored Apr 5, 2023
2 parents b526111 + af0846f commit 62a40fc
Show file tree
Hide file tree
Showing 7 changed files with 246 additions and 58 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
2 changes: 1 addition & 1 deletion src/felupe/mesh/_convert.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
102 changes: 102 additions & 0 deletions src/felupe/mesh/_discrete_geometry.py
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
"""

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)
15 changes: 8 additions & 7 deletions src/felupe/mesh/_helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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:]

Expand Down Expand Up @@ -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)

Expand Down
161 changes: 112 additions & 49 deletions src/felupe/mesh/_mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
"""

Expand All @@ -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 = "<felupe mesh object>"
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
Expand Down Expand Up @@ -160,13 +148,88 @@ def copy(self):

return deepcopy(self)

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}"]
@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
3 changes: 2 additions & 1 deletion src/felupe/mesh/_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@

from ..math import rotation_matrix
from ._helpers import mesh_or_data
from ._mesh import Mesh


@mesh_or_data
Expand Down Expand Up @@ -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)])
Expand Down
Loading

0 comments on commit 62a40fc

Please sign in to comment.