From 915fec99c0412723810a6924fb56f4f8ce2f3902 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 5 Mar 2024 10:15:02 +0200 Subject: [PATCH 1/6] mesh: make Mesh into a dataclass --- doc/conf.py | 2 +- meshmode/discretization/__init__.py | 3 +- meshmode/mesh/__init__.py | 747 +++++++++++++++++++--------- meshmode/mesh/processing.py | 63 ++- test/test_mesh.py | 2 +- 5 files changed, 543 insertions(+), 274 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index 64cd3677..840609f7 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -28,7 +28,7 @@ "arraycontext": ("https://documen.tician.de/arraycontext", None), "fenics": ("https://fenics.readthedocs.io/projects/fiat/en/latest", None), "FInAT": ("https://finat.github.io/FInAT/", None), - "firedrake": ("https://firedrakeproject.org", None), + "firedrake": ("https://www.firedrakeproject.org", None), "gmsh_interop": ("https://documen.tician.de/gmsh_interop", None), "h5py": ("https://docs.h5py.org/en/stable", None), "loopy": ("https://documen.tician.de/loopy", None), diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index a14d02d8..ba8df8dd 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -479,8 +479,7 @@ def copy(self, group_factory: Optional[ElementGroupFactory] = None, real_dtype: Optional[np.dtype] = None) -> "Discretization": """Creates a new object of the same type with all arguments that are not - *None* replaced. The copy is not recursive (e.g. it does not call - :meth:`meshmode.mesh.Mesh.copy`). + *None* replaced. The copy is not recursive. """ return type(self)( diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index e019a6e8..d0bce19f 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -21,26 +21,29 @@ """ from abc import ABC, abstractmethod -from dataclasses import dataclass, field -from typing import Any, ClassVar, Hashable, Optional, Sequence, Tuple, Type +from dataclasses import InitVar, dataclass, field, replace +from typing import ( + Any, ClassVar, Hashable, Iterable, Literal, Optional, Sequence, Tuple, Type, + Union) from warnings import warn import numpy as np import numpy.linalg as la import modepy as mp -from pytools import Record, memoize_method +from pytools import memoize_method from meshmode.mesh.tools import AffineMap __doc__ = """ - .. autoclass:: MeshElementGroup .. autoclass:: SimplexElementGroup .. autoclass:: TensorProductElementGroup .. autoclass:: Mesh +.. autofunction:: make_mesh +.. autofunction:: check_mesh_consistency .. autoclass:: NodalAdjacency .. autoclass:: FacialAdjacencyGroup @@ -730,52 +733,490 @@ def as_python(self): # {{{ mesh -class Mesh(Record): - r""" - .. attribute:: ambient_dim +DTypeLike = Union[np.dtype, np.generic] +NodalAdjacencyLike = Union[ + Literal[False], Iterable[np.ndarray], NodalAdjacency +] +FacialAdjacencyLike = Union[ + Literal[False], Iterable[Iterable[FacialAdjacencyGroup]] +] + + +def check_mesh_consistency( + mesh: "Mesh", + *, + node_vertex_consistency_tolerance: Optional[ + Union[Literal[False], float]] = None, + skip_element_orientation_test: bool = False, + ) -> None: + """Check the mesh for consistency between the vertices, nodes, and their + adjacency. + + This function checks: + + * The node to vertex consistency, by interpolation. + * The :class:`~numpy.dtype` of the various arrays matching the ones in + :class:`Mesh`. + * The nodal adjacency shapes and dtypes. + * The facial adjacency shapes and dtypes. + * The mesh orientation using + :func:`~meshmode.mesh.processing.find_volume_mesh_element_orientations`. + """ + if node_vertex_consistency_tolerance is not False: + assert _test_node_vertex_consistency( + mesh, node_vertex_consistency_tolerance) + + for g in mesh.groups: + if g.vertex_indices is not None: + assert g.vertex_indices.dtype == mesh.vertex_id_dtype + + nodal_adjacency = mesh._nodal_adjacency + if nodal_adjacency: + assert nodal_adjacency.neighbors_starts.shape == (mesh.nelements + 1,) + assert len(nodal_adjacency.neighbors.shape) == 1 + assert nodal_adjacency.neighbors_starts.dtype == mesh.element_id_dtype + assert nodal_adjacency.neighbors.dtype == mesh.element_id_dtype + + facial_adjacency_groups = mesh._facial_adjacency_groups + if facial_adjacency_groups: + assert len(facial_adjacency_groups) == len(mesh.groups) + + for fagrp_list in facial_adjacency_groups: + for fagrp in fagrp_list: + nfagrp_elements, = fagrp.elements.shape + assert fagrp.element_faces.dtype == mesh.face_id_dtype + assert fagrp.element_faces.shape == (nfagrp_elements,) + + if isinstance(fagrp, InteriorAdjacencyGroup): + assert fagrp.neighbors.dtype == mesh.element_id_dtype + assert fagrp.neighbors.shape == (nfagrp_elements,) + assert fagrp.neighbor_faces.dtype == mesh.face_id_dtype + assert fagrp.neighbor_faces.shape == (nfagrp_elements,) + + from meshmode.mesh.processing import test_volume_mesh_element_orientations + + if mesh.dim == mesh.ambient_dim and not skip_element_orientation_test: + assert test_volume_mesh_element_orientations(mesh) + + +def make_mesh( + groups: Iterable[MeshElementGroup], + *, + vertices: Optional[np.ndarray] = None, + nodal_adjacency: Optional[NodalAdjacencyLike] = None, + facial_adjacency_groups: Optional[FacialAdjacencyLike] = None, + is_conforming: Optional[bool] = None, + # dtypes + vertex_id_dtype: DTypeLike = "int32", + element_id_dtype: DTypeLike = "int32", + face_id_dtype: DTypeLike = "int8", + # tests + skip_tests: bool = False, + node_vertex_consistency_tolerance: Optional[float] = None, + skip_element_orientation_test: bool = False, + ) -> "Mesh": + """Construct a new mesh from a given list of *groups*. + + This constructor performs additional checks on the mesh once constructed and + should be preferred to calling the constructor of the :class:`Mesh` class + directly. + + :arg vertices: an array of vertices that match the given element *groups*. + These can be *None* for meshes where adjacency is not required + (e.g. non-conforming meshes). + :arg nodal_adjacency: a definition of the nodal adjacency of the mesh. + This argument can take one of four values: + + * *False*, in which case the information is marked as unavailable for + this mesh and will not be computed at all. This should be used if the + vertex adjacency does not convey the full picture, e.g if there are + hanging nodes in the geometry. + * *None*, in which case the nodal adjacency will be deduced from the + vertex adjacency on demand (this requires the *vertices*). + * a tuple of ``(element_neighbors_starts, element_neighbors)`` from which + a :class:`NodalAdjacency` object can be constructed. + * a :class:`NodalAdjacency` object. + + :arg facial_adjacency_groups: a definition of the facial adjacency for + each group in the mesh. This argument can take one of three values: + + * *False*, in which case the information is marked as unavailable for + this mesh and will not be computed. + * *None*, in which case the facial adjacency will be deduced from the + vertex adjacency on demand (this requires *vertices*). + * an iterable of :class:`FacialAdjacencyGroup` objects. + + :arg is_conforming: *True* if the mesh is known to be conforming. + + :arg vertex_id_dtype: an integer :class:`~numpy.dtype` for the vertex indices. + :arg element_id_dtype: an integer :class:`~numpy.dtype` for the element indices + (relative to an element group). + :arg face_id_dtype: an integer :class:`~numpy.dtype` for the face indices + (relative to an element). + + :arg skip_tests: a flag used to skip any mesh consistency checks. This can + be set to *True* in special situation, e.g. when loading a broken mesh + that will be fixed later. + :arg node_vertex_consistency_tolerance: If *False*, do not check for + consistency between vertex and nodal data. If *None*, a default tolerance + based on the :class:`~numpy.dtype` of the *vertices* array will be used. + :arg skip_element_orientation_test: If *False*, check that element + orientation is positive in volume meshes (i.e. ones where ambient and + topological dimension match). + """ + vertex_id_dtype = np.dtype(vertex_id_dtype) + if vertex_id_dtype.kind not in {"i", "u"}: + raise ValueError( + f"'vertex_id_dtype' expected to be an integer kind: {vertex_id_dtype}" + ) + + element_id_dtype = np.dtype(element_id_dtype) + if element_id_dtype.kind not in {"i", "u"}: + raise ValueError( + f"'element_id_dtype' expected to be an integer kind: {element_id_dtype}" + ) + + face_id_dtype = np.dtype(face_id_dtype) + if face_id_dtype.kind not in {"i", "u"}: + raise ValueError( + f"'face_id_dtype' expected to be an integer kind: {face_id_dtype}" + ) + + if vertices is None: + if is_conforming is not None: + warn("No vertices provided and 'is_conforming' is set to " + f"'{is_conforming}'. Setting to 'None' instead, since no " + "adjacency can be known.", + UserWarning, stacklevel=2) + + is_conforming = None + + if not is_conforming: + if nodal_adjacency is None: + nodal_adjacency = False + + if facial_adjacency_groups is None: + facial_adjacency_groups = False + + if ( + nodal_adjacency is not False + and nodal_adjacency is not None + and not isinstance(nodal_adjacency, NodalAdjacency)): + nb_starts, nbs = nodal_adjacency + nodal_adjacency = ( + NodalAdjacency(neighbors_starts=nb_starts, neighbors=nbs)) + + if ( + facial_adjacency_groups is not False + and facial_adjacency_groups is not None): + facial_adjacency_groups = _complete_facial_adjacency_groups( + facial_adjacency_groups, + element_id_dtype, + face_id_dtype) + facial_adjacency_groups = tuple([ + tuple(grps) for grps in facial_adjacency_groups + ]) + + mesh = Mesh( + groups=tuple(groups), + vertices=vertices, + is_conforming=is_conforming, + vertex_id_dtype=vertex_id_dtype, + element_id_dtype=element_id_dtype, + face_id_dtype=face_id_dtype, + _nodal_adjacency=nodal_adjacency, + _facial_adjacency_groups=facial_adjacency_groups, + factory_constructed=True + ) + + if __debug__ and not skip_tests: + check_mesh_consistency( + mesh, + node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, + skip_element_orientation_test=skip_element_orientation_test) + + return mesh + + +# TODO: should be `init=True` once everything is ported to `make_mesh` +@dataclass(frozen=True, init=False, eq=False) +class Mesh: + """ + .. autoproperty:: ambient_dim + .. autoproperty:: dim + .. autoproperty:: nvertices + .. autoproperty:: nelements + .. autoproperty:: base_element_nrs + .. autoproperty:: base_node_nrs + .. autoproperty:: vertex_dtype - .. attribute:: dim + .. autoattribute :: groups + .. autoattribute:: vertices + .. autoattribute:: is_conforming - .. attribute:: vertices + .. autoattribute:: vertex_id_dtype + .. autoattribute:: element_id_dtype + .. autoattribute:: face_id_dtype - *None* or an array of vertex coordinates with shape - *(ambient_dim, nvertices)*. If *None*, vertices are not - known for this mesh. + .. autoproperty:: nodal_adjacency + .. autoproperty:: facial_adjacency_groups - .. attribute:: nvertices + .. autoattribute:: _nodal_adjacency + .. autoattribute:: _facial_adjacency_groups - .. attribute:: groups + .. automethod:: __eq__ + """ - A list of :class:`MeshElementGroup` instances. + groups: Tuple[MeshElementGroup, ...] + """A tuple of :class:`MeshElementGroup` instances.""" - .. attribute:: nelements + vertices: Optional[np.ndarray] + """*None* or an array of vertex coordinates with shape + *(ambient_dim, nvertices)*. If *None*, vertices are not known for this mesh + and no adjacency information can be constructed. + """ + + is_conforming: Optional[bool] + """*True* if it is known that all element interfaces are conforming. *False* + if it is known that some element interfaces are non-conforming. *None* if + neither of the two is known. + """ - .. attribute:: base_element_nrs + vertex_id_dtype: np.dtype + """The :class:`~numpy.dtype` used to index into the vertex array.""" - An array of size ``(len(groups),)`` of starting element indices for - each group in the mesh. + element_id_dtype: np.dtype + """The :class:`~numpy.dtype` used to index into the element array (relative + to each group). + """ + + face_id_dtype: np.dtype + """The :class:`~numpy.dtype` used to index element faces (relative to each + element). + """ + + # TODO: Once the @property(nodal_adjacency) is past its deprecation period + # and removed, these can deprecated in favor of non-underscored variants. + + _nodal_adjacency: Union[None, Literal[False], NodalAdjacency] + """A description of the nodal adjacency of the mesh. This can be *False* if + no adjacency is known or should be computed, *None* to compute the adjacency + on demand or a given :class:`NodalAdjacency` instance. + + This attribute caches the values of :attr:`nodal_adjacency`. + """ + + _facial_adjacency_groups: \ + Union[None, Literal[False], Tuple[Tuple[FacialAdjacencyGroup, ...], ...]] + """A description of the facial adjacency of the mesh. This can be *False* if + no adjacency is known or should be computed, *None* to compute the adjacency + on demand or a list of :class:`FacialAdjacencyGroup` instances. + + This attribute caches the values of :attr:`facial_adjacency_groups`. + """ + + # TODO: remove once porting to `make_mesh` is complete. + skip_tests: InitVar[bool] = False + node_vertex_consistency_tolerance: InitVar[ + Optional[Union[Literal[False], float]]] = None + skip_element_orientation_test: InitVar[bool] = False + factory_constructed: InitVar[bool] = False + + def __init__( + self, + vertices: Optional[np.ndarray], + groups: Iterable[MeshElementGroup], + is_conforming: Optional[bool] = None, + vertex_id_dtype: DTypeLike = "int32", + element_id_dtype: DTypeLike = "int32", + face_id_dtype: DTypeLike = "int8", + # cached variables + nodal_adjacency: Optional[NodalAdjacencyLike] = None, + facial_adjacency_groups: Optional[FacialAdjacencyLike] = None, + _nodal_adjacency: Optional[NodalAdjacencyLike] = None, + _facial_adjacency_groups: Optional[FacialAdjacencyLike] = None, + # init vars + skip_tests: bool = False, + node_vertex_consistency_tolerance: Optional[float] = None, + skip_element_orientation_test: bool = False, + factory_constructed: bool = False, + ) -> None: + if _nodal_adjacency is None: + if nodal_adjacency is not None: + warn("Passing 'nodal_adjacency' is deprecated and will be removed " + "in 2025. Use the underscore '_nodal_adjacency' instead to " + "match the dataclass field.", + DeprecationWarning, stacklevel=2) + + _nodal_adjacency = nodal_adjacency + + if _facial_adjacency_groups is None: + if facial_adjacency_groups is not None: + warn("Passing 'facial_adjacency_groups' is deprecated and will be " + "removed in 2025. Use the underscore '_facial_adjacency_groups'" + " instead to match the dataclass field.", + DeprecationWarning, stacklevel=2) + + _facial_adjacency_groups = facial_adjacency_groups + + if not factory_constructed: + warn(f"Calling '{type(self).__name__}(...)' constructor is deprecated. " + "Use the 'make_mesh(...)' factory function instead. The input " + "handling in the constructor will be removed in 2025.", + DeprecationWarning, stacklevel=2) + + vertex_id_dtype = np.dtype(vertex_id_dtype) + element_id_dtype = np.dtype(element_id_dtype) + face_id_dtype = np.dtype(face_id_dtype) + + if vertices is None: + is_conforming = None + + if not is_conforming: + if _nodal_adjacency is None: + _nodal_adjacency = False + if _facial_adjacency_groups is None: + _facial_adjacency_groups = False + + if ( + _nodal_adjacency is not False + and _nodal_adjacency is not None): + if not isinstance(_nodal_adjacency, NodalAdjacency): + nb_starts, nbs = _nodal_adjacency + _nodal_adjacency = NodalAdjacency( + neighbors_starts=nb_starts, + neighbors=nbs) + + del nb_starts + del nbs + + if ( + _facial_adjacency_groups is not False + and _facial_adjacency_groups is not None): + _facial_adjacency_groups = _complete_facial_adjacency_groups( + _facial_adjacency_groups, + element_id_dtype, + face_id_dtype) + + object.__setattr__(self, "groups", tuple(groups)) + object.__setattr__(self, "vertices", vertices) + object.__setattr__(self, "is_conforming", is_conforming) + object.__setattr__(self, "vertex_id_dtype", vertex_id_dtype) + object.__setattr__(self, "element_id_dtype", element_id_dtype) + object.__setattr__(self, "face_id_dtype", face_id_dtype) + object.__setattr__(self, "_nodal_adjacency", _nodal_adjacency) + object.__setattr__(self, "_facial_adjacency_groups", + _facial_adjacency_groups) + + if __debug__ and not factory_constructed and not skip_tests: + check_mesh_consistency( + self, + node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, + skip_element_orientation_test=skip_element_orientation_test) + + def copy(self, **kwargs: Any) -> "Mesh": + warn(f"'{type(self).__name__}.copy' is deprecated and will be removed in " + f"2025. '{type(self).__name__}' is a dataclass and can use the " + "standard 'replace' function.", + DeprecationWarning, stacklevel=2) + + if "nodal_adjacency" in kwargs: + kwargs["_nodal_adjacency"] = kwargs.pop("nodal_adjacency") + + if "facial_adjacency_groups" in kwargs: + kwargs["_facial_adjacency_groups"] = ( + kwargs.pop("facial_adjacency_groups")) + + return replace(self, **kwargs) + + @property + def ambient_dim(self) -> int: + """Ambient dimension in which the mesh is embedded.""" + from pytools import single_valued + return single_valued(grp.nodes.shape[0] for grp in self.groups) + + @property + def dim(self) -> int: + """Dimension of the elements in the mesh.""" + from pytools import single_valued + return single_valued(grp.dim for grp in self.groups) - .. attribute:: base_node_nrs + @property + def nvertices(self) -> int: + """Number of vertices in the mesh, if available.""" + if self.vertices is None: + from meshmode import DataUnavailable + raise DataUnavailable("vertices") - An array of size ``(len(groups),)`` of starting node indices for + return self.vertices.shape[-1] + + @property + def nelements(self) -> int: + """Number of elements in the mesh (sum over all the :attr:`~Mesh.groups`).""" + return sum(grp.nelements for grp in self.groups) + + @property + @memoize_method + def base_element_nrs(self) -> np.ndarray: + """An array of size ``(len(groups),)`` of starting element indices for each group in the mesh. + """ + return np.cumsum([0] + [grp.nelements for grp in self.groups[:-1]]) + + @property + @memoize_method + def base_node_nrs(self): + """An array of size ``(len(groups),)`` of starting node indices for + each group in the mesh. + """ + return np.cumsum([0] + [grp.nnodes for grp in self.groups[:-1]]) + + @property + def vertex_dtype(self): + """The :class:`~numpy.dtype` of the :attr:`~Mesh.vertices` array, if any.""" + if self.vertices is None: + from meshmode import DataUnavailable + raise DataUnavailable("vertices") + + return self.vertices.dtype + + @property + def nodal_adjacency(self) -> NodalAdjacency: + """Nodal adjacency of the mesh, if available. + + This property gets the :attr:`Mesh._nodal_adjacency` of the mesh. If the + attribute value is *None*, the adjacency is computed and cached. + + :raises DataUnavailable: if the nodal adjacency cannot be obtained. + """ + from meshmode import DataUnavailable - .. attribute:: nodal_adjacency + nodal_adjacency = self._nodal_adjacency + if nodal_adjacency is False: + raise DataUnavailable("Nodal adjacency is not available") - An instance of :class:`NodalAdjacency`. + if nodal_adjacency is None: + if not self.is_conforming: + raise DataUnavailable( + "Nodal adjacency can only be computed for conforming meshes" + ) - Referencing this attribute may raise - :exc:`meshmode.DataUnavailable`. + nodal_adjacency = _compute_nodal_adjacency_from_vertices(self) + object.__setattr__(self, "_nodal_adjacency", nodal_adjacency) - .. attribute:: facial_adjacency_groups + return nodal_adjacency - A list of lists of instances of :class:`FacialAdjacencyGroup`. + @property + def facial_adjacency_groups( + self) -> Tuple[Tuple[FacialAdjacencyGroup, ...], ...]: + r"""Facial adjacency of the mesh, if available. - ``facial_adjacency_groups[igrp]`` gives the facial adjacency relations for - group *igrp*, expressed as a list of :class:`FacialAdjacencyGroup` instances. + This function gets the :attr:`Mesh._facial_adjacency_groups` of the mesh. + If the attribute value is *None*, the adjacency is computed and cached. - Referencing this attribute may raise - :exc:`meshmode.DataUnavailable`. + Each ``facial_adjacency_groups[igrp]`` gives the facial adjacency + relations for group *igrp*, expressed as a list of + :class:`FacialAdjacencyGroup` instances. .. tikz:: Facial Adjacency Group :align: center @@ -792,7 +1233,6 @@ class Mesh(Record): \draw [line width=3pt, line cap=round, green!60!black] (4, 2) -- (6, 2); - For example for the mesh in the figure, the following data structure could be present:: @@ -811,235 +1251,46 @@ class Mesh(Record): ] ] - (Note that element groups are not necessarily geometrically contiguous - like the figure may suggest.) - - .. attribute:: vertex_id_dtype + Note that element groups are not necessarily geometrically contiguous + like the figure may suggest. - .. attribute:: element_id_dtype - - .. attribute:: is_conforming - - *True* if it is known that all element interfaces are conforming. - *False* if it is known that some element interfaces are non-conforming. - *None* if neither of the two is known. - - .. automethod:: copy - .. automethod:: __eq__ - .. automethod:: __ne__ - """ - - groups: Sequence[MeshElementGroup] - - face_id_dtype = np.int8 - - def __init__(self, vertices, groups, *, skip_tests=False, - node_vertex_consistency_tolerance=None, - skip_element_orientation_test=False, - nodal_adjacency=None, - facial_adjacency_groups=None, - vertex_id_dtype=np.int32, - element_id_dtype=np.int32, - is_conforming=None): - """ - :arg skip_tests: Skip mesh tests, in case you want to load a broken - mesh anyhow and then fix it inside of this data structure. - :arg node_vertex_consistency_tolerance: If *False*, do not check - for consistency between vertex and nodal data. If *None*, use - the (small, near FP-epsilon) default tolerance. - :arg skip_element_orientation_test: If *False*, check that - element orientation is positive in volume meshes - (i.e. ones where ambient and topological dimension match). - :arg nodal_adjacency: One of three options: - *None*, in which case this information - will be deduced from vertex adjacency. *False*, in which case - this information will be marked unavailable (such as if there are - hanging nodes in the geometry, so that vertex adjacency does not convey - the full picture), and references to - :attr:`element_neighbors_starts` and :attr:`element_neighbors` - will result in exceptions. Lastly, a tuple - :class:`NodalAdjacency` object. - :arg facial_adjacency_groups: One of three options: - *None*, in which case this information - will be deduced from vertex adjacency. *False*, in which case - this information will be marked unavailable (such as if there are - hanging nodes in the geometry, so that vertex adjacency does not convey - the full picture), and references to - :attr:`element_neighbors_starts` and :attr:`element_neighbors` - will result in exceptions. Lastly, a data structure as described in - :attr:`facial_adjacency_groups` may be passed. + :raises DataUnavailable: if the facial adjacency cannot be obtained. """ - - if vertices is None: - is_conforming = None - - if not is_conforming: - if nodal_adjacency is None: - nodal_adjacency = False - if facial_adjacency_groups is None: - facial_adjacency_groups = False - - if nodal_adjacency is not False and nodal_adjacency is not None: - if not isinstance(nodal_adjacency, NodalAdjacency): - nb_starts, nbs = nodal_adjacency - nodal_adjacency = NodalAdjacency( - neighbors_starts=nb_starts, - neighbors=nbs) - - del nb_starts - del nbs - - if ( - facial_adjacency_groups is not False - and facial_adjacency_groups is not None): - facial_adjacency_groups = _complete_facial_adjacency_groups( - facial_adjacency_groups, - np.dtype(element_id_dtype), - self.face_id_dtype) - - Record.__init__( - self, vertices=vertices, groups=groups, - _nodal_adjacency=nodal_adjacency, - _facial_adjacency_groups=facial_adjacency_groups, - vertex_id_dtype=np.dtype(vertex_id_dtype), - element_id_dtype=np.dtype(element_id_dtype), - is_conforming=is_conforming, - ) - - if not skip_tests: - if node_vertex_consistency_tolerance is not False: - assert _test_node_vertex_consistency( - self, node_vertex_consistency_tolerance) - - for g in self.groups: - if g.vertex_indices is not None: - assert g.vertex_indices.dtype == self.vertex_id_dtype - - if nodal_adjacency: - assert nodal_adjacency.neighbors_starts.shape == (self.nelements+1,) - assert len(nodal_adjacency.neighbors.shape) == 1 - - assert (nodal_adjacency.neighbors_starts.dtype - == self.element_id_dtype) - assert nodal_adjacency.neighbors.dtype == self.element_id_dtype - - if facial_adjacency_groups: - assert len(facial_adjacency_groups) == len(self.groups) - for fagrp_list in facial_adjacency_groups: - for fagrp in fagrp_list: - nfagrp_elements, = fagrp.elements.shape - assert fagrp.element_faces.dtype == self.face_id_dtype - assert fagrp.element_faces.shape == (nfagrp_elements,) - if isinstance(fagrp, InteriorAdjacencyGroup): - assert fagrp.neighbors.dtype == self.element_id_dtype - assert fagrp.neighbors.shape == (nfagrp_elements,) - assert fagrp.neighbor_faces.dtype == self.face_id_dtype - assert fagrp.neighbor_faces.shape == (nfagrp_elements,) - - from meshmode.mesh.processing import ( - test_volume_mesh_element_orientations) - - if self.dim == self.ambient_dim and not skip_element_orientation_test: - # only for volume meshes, for now - if not test_volume_mesh_element_orientations(self): - raise ValueError("negatively oriented elements found") - - def get_copy_kwargs(self, **kwargs): - def set_if_not_present(name, from_name=None): - if from_name is None: - from_name = name - if name not in kwargs: - kwargs[name] = getattr(self, from_name) - - set_if_not_present("vertices") - if "groups" not in kwargs: - kwargs["groups"] = self.groups - - set_if_not_present("nodal_adjacency", "_nodal_adjacency") - set_if_not_present("facial_adjacency_groups", "_facial_adjacency_groups") - set_if_not_present("vertex_id_dtype") - set_if_not_present("element_id_dtype") - set_if_not_present("is_conforming") - - return kwargs - - @property - def ambient_dim(self): - from pytools import single_valued - return single_valued(grp.nodes.shape[0] for grp in self.groups) - - @property - def dim(self): - from pytools import single_valued - return single_valued(grp.dim for grp in self.groups) - - @property - def nvertices(self): - if self.vertices is None: - from meshmode import DataUnavailable - raise DataUnavailable("vertices") - - return self.vertices.shape[-1] - - @property - def nelements(self): - return sum(grp.nelements for grp in self.groups) - - @property - @memoize_method - def base_element_nrs(self): - return np.cumsum([0] + [grp.nelements for grp in self.groups[:-1]]) - - @property - @memoize_method - def base_node_nrs(self): - return np.cumsum([0] + [grp.nnodes for grp in self.groups[:-1]]) - - @property - def nodal_adjacency(self): from meshmode import DataUnavailable - # pylint: disable=access-member-before-definition - if self._nodal_adjacency is False: - raise DataUnavailable("nodal_adjacency") + facial_adjacency_groups = self._facial_adjacency_groups + if facial_adjacency_groups is False: + raise DataUnavailable("Facial adjacency is not available") - elif self._nodal_adjacency is None: + if facial_adjacency_groups is None: if not self.is_conforming: - raise DataUnavailable("nodal_adjacency can only " - "be computed for known-conforming meshes") - - self._nodal_adjacency = _compute_nodal_adjacency_from_vertices(self) - - return self._nodal_adjacency - - def nodal_adjacency_init_arg(self): - """Returns a *nodal_adjacency* argument that can be - passed to a :class:`Mesh` constructor. - """ + raise DataUnavailable( + "Facial adjacency can only be computed for conforming meshes" + ) - return self._nodal_adjacency + facial_adjacency_groups = _compute_facial_adjacency_from_vertices( + self.groups, self.element_id_dtype, self.face_id_dtype) + object.__setattr__(self, "_facial_adjacency_groups", + facial_adjacency_groups) - @property - def facial_adjacency_groups(self) -> Sequence[Sequence[FacialAdjacencyGroup]]: - from meshmode import DataUnavailable + return facial_adjacency_groups - # pylint: disable=access-member-before-definition - if self._facial_adjacency_groups is False: - raise DataUnavailable("facial_adjacency_groups") + def __eq__(self, other): + """Compare two meshes for equality. - elif self._facial_adjacency_groups is None: - if not self.is_conforming: - raise DataUnavailable("facial_adjacency_groups can only " - "be computed for known-conforming meshes") + .. warning:: - self._facial_adjacency_groups = _compute_facial_adjacency_from_vertices( - self.groups, - self.element_id_dtype, - self.face_id_dtype) + This operation is very expensive, as it compares all the vertices and + groups between the two meshes. If available, the nodal and facial + adjacency information is compared as well. - return self._facial_adjacency_groups + .. warning:: - def __eq__(self, other): + Only the (uncached) :attr:`~Mesh._nodal_adjacency` and + :attr:`~Mesh._facial_adjacency_groups` are compared. This can fail + for two meshes if one called :meth:`~Mesh.nodal_adjacency` + and the other one did not, even if they would be equal. + """ return ( type(self) is type(other) and np.array_equal(self.vertices, other.vertices) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index e1d387a2..88133fe0 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -25,7 +25,7 @@ from dataclasses import dataclass, replace from functools import reduce from typing import ( - Callable, Dict, List, Mapping, Optional, Sequence, Set, Tuple, Union) + Callable, Dict, List, Literal, Mapping, Optional, Sequence, Set, Tuple, Union) import numpy as np import numpy.linalg as la @@ -451,6 +451,9 @@ def _get_mesh_part( .. versionadded:: 2017.1 """ + if mesh.vertices is None: + raise ValueError("Mesh must have vertices") + element_counts = np.zeros(mesh.nelements) for elements in part_id_to_elements.values(): element_counts[elements] += 1 @@ -608,6 +611,8 @@ def find_volume_mesh_element_orientations( :arg tolerate_unimplemented_checks: If *True*, elements for which no check is available will return *NaN*. """ + if mesh.vertices is None: + raise ValueError("Mesh must have vertices to check orientation") result: np.ndarray = np.empty(mesh.nelements, dtype=np.float64) @@ -743,6 +748,8 @@ def perform_flips( indicating by their Boolean value whether the element is to be flipped. """ + if mesh.vertices is None: + raise ValueError("Mesh must have vertices to perform flips") flip_flags = flip_flags.astype(bool) @@ -772,6 +779,8 @@ def find_bounding_box(mesh: Mesh) -> Tuple[np.ndarray, np.ndarray]: :return: a tuple *(min, max)*, each consisting of a :class:`numpy.ndarray` indicating the minimal and maximal extent of the geometry along each axis. """ + if mesh.vertices is None: + raise ValueError("Mesh must have vertices to compute bounding box") return ( np.min(mesh.vertices, axis=-1), @@ -796,34 +805,35 @@ def merge_disjoint_meshes( # {{{ assemble combined vertex array - ambient_dim = meshes[0].ambient_dim - nvertices = sum( - mesh.vertices.shape[-1] - for mesh in meshes) + if all(mesh.vertices is not None for mesh in meshes): + ambient_dim = meshes[0].ambient_dim + nvertices = sum(mesh.nvertices for mesh in meshes) - vert_dtype = np.result_type(*[mesh.vertices.dtype for mesh in meshes]) - vertices = np.empty( - (ambient_dim, nvertices), vert_dtype) + vert_dtype = np.result_type(*[mesh.vertex_dtype for mesh in meshes]) + vertices = np.empty((ambient_dim, nvertices), vert_dtype) - current_vert_base = 0 - vert_bases = [] - for mesh in meshes: - mesh_nvert = mesh.vertices.shape[-1] - vertices[:, current_vert_base:current_vert_base+mesh_nvert] = \ - mesh.vertices + current_vert_base = 0 + vert_bases = [] + for mesh in meshes: + assert mesh.vertices is not None + mesh_nvert = mesh.nvertices + vertices[:, current_vert_base:current_vert_base+mesh_nvert] = \ + mesh.vertices - vert_bases.append(current_vert_base) - current_vert_base += mesh_nvert + vert_bases.append(current_vert_base) + current_vert_base += mesh_nvert + else: + raise ValueError("All meshes must have vertices to perform merge") # }}} # {{{ assemble new groups list - nodal_adjacency = None + nodal_adjacency: Optional[Literal[False]] = None if any(mesh._nodal_adjacency is not None for mesh in meshes): nodal_adjacency = False - facial_adjacency_groups = None + facial_adjacency_groups: Optional[Literal[False]] = None if any(mesh._facial_adjacency_groups is not None for mesh in meshes): facial_adjacency_groups = False @@ -948,6 +958,9 @@ def _match_vertices( aff_map: Optional[AffineMap] = None, tol: float = 1e-12, use_tree: Optional[bool] = None) -> np.ndarray: + if mesh.vertices is None: + raise ValueError("Mesh must have vertices") + if aff_map is None: aff_map = AffineMap() @@ -1106,6 +1119,9 @@ def _match_boundary_faces( second contains faces from boundary *bdry_pair_mapping.to_btag*. The order of the faces is unspecified. """ + if mesh.vertices is None: + raise ValueError("Mesh must have vertices") + btag_m = bdry_pair_mapping.from_btag btag_n = bdry_pair_mapping.to_btag @@ -1307,9 +1323,12 @@ def map_mesh(mesh: Mesh, f: Callable[[np.ndarray], np.ndarray]) -> Mesh: "affine mappings in its facial adjacency. If the map is affine, " "use affine_map instead") - vertices = f(mesh.vertices) - if not vertices.flags.c_contiguous: - vertices = np.copy(vertices, order="C") + if mesh.vertices is not None: + vertices = f(mesh.vertices) + if not vertices.flags.c_contiguous: + vertices = np.copy(vertices, order="C") + else: + vertices = None # {{{ assemble new groups list @@ -1512,7 +1531,7 @@ def make_mesh_grid( meshes = [] for index in product(*(range(n) for n in shape)): - b = sum([i * o for i, o in zip(index, offset)]) + b = sum([i * o for i, o in zip(index, offset)], offset[0]) meshes.append(affine_map(mesh, b=b)) return merge_disjoint_meshes(meshes, skip_tests=skip_tests) diff --git a/test/test_mesh.py b/test/test_mesh.py index 026f4df7..a474d30c 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -831,7 +831,7 @@ def test_boundary_tags(): def test_volume_tags(): from meshmode.mesh.io import read_gmsh - mesh, tag_to_elements_map = read_gmsh( + mesh, tag_to_elements_map = read_gmsh( # pylint: disable=unpacking-non-sequence str(thisdir / "testmesh_multivol.msh"), return_tag_to_elements_map=True) assert len(tag_to_elements_map) == 2 From be847eaf66dfc772034cdc503a8989175e893dc8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 7 Apr 2024 20:17:29 +0300 Subject: [PATCH 2/6] firedrake: set mesh dtypes explicitly --- meshmode/interop/firedrake/mesh.py | 31 +++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py index 35565c5c..e7b9411f 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -214,7 +214,8 @@ def _get_facet_markers(dm, facets): def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, - cells_to_use=None): + cells_to_use=None, + face_id_dtype=None): """ Return facial_adjacency_groups corresponding to the given firedrake mesh topology. Note that as we do not @@ -236,6 +237,8 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, by a :mod:`meshmode` :class:`Mesh`. """ top = fdrake_mesh_topology.topology + if face_id_dtype is None: + face_id_dtype = np.int8 # We only need one group # for interconnectivity and one for boundary connectivity. @@ -283,10 +286,10 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, int_elements = int_facet_cell.flatten() int_neighbors = np.concatenate((int_facet_cell[:, 1], int_facet_cell[:, 0])) - int_element_faces = int_fac_loc_nr.flatten().astype(Mesh.face_id_dtype) + int_element_faces = int_fac_loc_nr.flatten().astype(face_id_dtype) int_neighbor_faces = np.concatenate((int_fac_loc_nr[:, 1], int_fac_loc_nr[:, 0])) - int_neighbor_faces = int_neighbor_faces.astype(Mesh.face_id_dtype) + int_neighbor_faces = int_neighbor_faces.astype(face_id_dtype) # If only using some of the cells from pyop2.datatypes import IntType if cells_to_use is not None: @@ -344,7 +347,7 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, ext_element_faces = np.array([fd_loc_fac_nr_to_mm[fac_nr] for fac_nr in top.exterior_facets.local_facet_dat.data], - dtype=Mesh.face_id_dtype) + dtype=face_id_dtype) # If only using some of the cells, throw away unused cells and # move to new cell index @@ -707,11 +710,12 @@ def import_firedrake_mesh(fdrake_mesh, cells_to_use=None, elif 1 not in face: no_one_face_ndx = iface + face_id_dtype = np.int8 with ProcessLogger(logger, "Building (possibly) unflipped " "FacialAdjacencyGroups"): - unflipped_facial_adjacency_groups = \ - _get_firedrake_facial_adjacency_groups(fdrake_mesh, - cells_to_use=cells_to_use) + unflipped_facial_adjacency_groups = ( + _get_firedrake_facial_adjacency_groups( + fdrake_mesh, cells_to_use=cells_to_use, face_id_dtype=face_id_dtype)) # applied below to take elements and element_faces # (or neighbors and neighbor_faces) and flip in any faces that need to @@ -757,10 +761,15 @@ def flip_local_face_indices(faces, elements): elements=fagrp.elements, element_faces=new_element_faces)) - return (Mesh(vertices, [group], - nodal_adjacency=nodal_adjacency, - facial_adjacency_groups=facial_adjacency_groups), - orient) + mesh = Mesh( + vertices, [group], + vertex_id_dtype=vertex_indices.dtype, + element_id_dtype=vertex_indices.dtype, + face_id_dtype=face_id_dtype, + nodal_adjacency=nodal_adjacency, + facial_adjacency_groups=facial_adjacency_groups) + + return mesh, orient # }}} From b64e16dbf6b8e388c656059f43999f23cdca1568 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 Apr 2024 11:51:00 -0500 Subject: [PATCH 3/6] make_mesh: vertices remains first argument --- meshmode/mesh/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d0bce19f..0ea56276 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -800,9 +800,9 @@ def check_mesh_consistency( def make_mesh( + vertices: Optional[np.ndarray], groups: Iterable[MeshElementGroup], *, - vertices: Optional[np.ndarray] = None, nodal_adjacency: Optional[NodalAdjacencyLike] = None, facial_adjacency_groups: Optional[FacialAdjacencyLike] = None, is_conforming: Optional[bool] = None, From c1cc5e0ef833057464cfc9269c03d794638001e8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 Apr 2024 11:41:11 -0500 Subject: [PATCH 4/6] Port internal uses of Mesh constructor to make_mesh --- examples/tp-lagrange-stl.py | 4 +-- meshmode/discretization/connection/face.py | 4 +-- meshmode/interop/firedrake/mesh.py | 4 +-- meshmode/interop/nodal_dg.py | 4 +-- meshmode/mesh/__init__.py | 4 +-- meshmode/mesh/generation.py | 30 ++++++++-------------- meshmode/mesh/io.py | 12 ++++----- meshmode/mesh/processing.py | 10 ++++---- meshmode/mesh/refinement/no_adjacency.py | 4 +-- test/test_mesh.py | 9 ++++--- test/test_meshmode.py | 9 ++++--- 11 files changed, 44 insertions(+), 50 deletions(-) diff --git a/examples/tp-lagrange-stl.py b/examples/tp-lagrange-stl.py index bd973bb4..063efa5d 100644 --- a/examples/tp-lagrange-stl.py +++ b/examples/tp-lagrange-stl.py @@ -12,7 +12,7 @@ import meshmode.mesh.generation as mgen from meshmode.array_context import PyOpenCLArrayContext from meshmode.discretization import Discretization -from meshmode.mesh import BTAG_ALL, Mesh +from meshmode.mesh import BTAG_ALL, make_mesh def main(): @@ -79,7 +79,7 @@ def main(): from meshmode.mesh.generation import make_group_from_vertices mod_grp = make_group_from_vertices(vertices, grp.vertex_indices, order=grp.order) - mod_mesh = Mesh( + mod_mesh = make_mesh( vertices=vertices, groups=[mod_grp], is_conforming=bdry_mesh.is_conforming) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 62d67cd1..4c6b6b01 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -217,7 +217,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # }}} - from meshmode.mesh import Mesh, _ModepyElementGroup + from meshmode.mesh import _ModepyElementGroup, make_mesh bdry_mesh_groups = [] connection_data = {} @@ -346,7 +346,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, unit_nodes=bdry_unit_nodes) bdry_mesh_groups.append(bdry_mesh_group) - bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups) + bdry_mesh = make_mesh(bdry_vertices, bdry_mesh_groups) bdry_discr = discr.copy( actx=actx, diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py index e7b9411f..a9a61983 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -32,7 +32,7 @@ get_affine_reference_simplex_mapping, get_finat_element_unit_nodes) from meshmode.mesh import ( BTAG_ALL, BTAG_INDUCED_BOUNDARY, BoundaryAdjacencyGroup, InteriorAdjacencyGroup, - Mesh, NodalAdjacency, SimplexElementGroup) + Mesh, NodalAdjacency, SimplexElementGroup, make_mesh) __doc__ = """ @@ -761,7 +761,7 @@ def flip_local_face_indices(faces, elements): elements=fagrp.elements, element_faces=new_element_faces)) - mesh = Mesh( + mesh = make_mesh( vertices, [group], vertex_id_dtype=vertex_indices.dtype, element_id_dtype=vertex_indices.dtype, diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py index ec539f91..5afc6270 100644 --- a/meshmode/interop/nodal_dg.py +++ b/meshmode/interop/nodal_dg.py @@ -146,7 +146,7 @@ def get_discr(self, actx) -> meshmode.discretization.Discretization: nodes = np.array([self.octave.pull(self.AXES[ax]).T for ax in range(dim)]) vertex_indices = (self.octave.pull("EToV")).astype(np.int32)-1 - from meshmode.mesh import Mesh, SimplexElementGroup + from meshmode.mesh import SimplexElementGroup, make_mesh order = int(self.octave.pull("N")) egroup = SimplexElementGroup.make_group( order, @@ -154,7 +154,7 @@ def get_discr(self, actx) -> meshmode.discretization.Discretization: nodes=nodes, unit_nodes=unit_nodes) - mesh = Mesh(vertices=vertices, groups=[egroup], is_conforming=True) + mesh = make_mesh(vertices=vertices, groups=[egroup], is_conforming=True) from meshmode.discretization import Discretization from meshmode.discretization.poly_element import ( diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 0ea56276..93185d3a 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1774,7 +1774,7 @@ def as_python(mesh, function_name="make_mesh"): import numpy as np from meshmode.mesh import ( - Mesh, + make_mesh as mm_make_mesh, MeshElementGroup, FacialAdjacencyGroup, InteriorAdjacencyGroup, @@ -1823,7 +1823,7 @@ def as_python(mesh, function_name="make_mesh"): # }}} - cg("return Mesh(vertices, groups, skip_tests=True,") + cg("return mm_make_mesh(vertices, groups, skip_tests=True,") cg(" vertex_id_dtype=np.%s," % mesh.vertex_id_dtype.name) cg(" element_id_dtype=np.%s," % mesh.element_id_dtype.name) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 9d5cebfc..9e76de67 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -31,7 +31,7 @@ import modepy as mp from pytools import deprecate_keyword, log_process -from meshmode.mesh import Mesh, MeshElementGroup +from meshmode.mesh import Mesh, MeshElementGroup, make_mesh from meshmode.mesh.refinement import Refiner @@ -429,14 +429,14 @@ def make_curve_mesh( t = t.ravel() nodes = curve_f(t).reshape(vertices.shape[0], nelements, -1) - from meshmode.mesh import Mesh, SimplexElementGroup + from meshmode.mesh import SimplexElementGroup egroup = SimplexElementGroup.make_group( order, vertex_indices=vertex_indices, nodes=nodes, unit_nodes=unit_nodes) - mesh = Mesh( + mesh = make_mesh( vertices=vertices, groups=[egroup], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -578,8 +578,7 @@ def generate_icosahedron( grp = make_group_from_vertices(vertices, vertex_indices, order, unit_nodes=unit_nodes) - from meshmode.mesh import Mesh - return Mesh( + return make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -601,8 +600,7 @@ def generate_cube_surface(r: float, order: int, *, group_cls=TensorProductElementGroup, unit_nodes=unit_nodes) - from meshmode.mesh import Mesh - return Mesh( + return make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -678,8 +676,7 @@ def generate_sphere(r: float, order: int, *, nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0)) ) - from meshmode.mesh import Mesh - return Mesh( + return make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -732,8 +729,7 @@ def generate_surface_of_revolution( grp = make_group_from_vertices(vertices, vertex_indices, order, unit_nodes=unit_nodes) - from meshmode.mesh import Mesh - mesh = Mesh( + mesh = make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -752,8 +748,7 @@ def ensure_radius(arr: np.ndarray) -> np.ndarray: grp, = mesh.groups grp = replace(grp, nodes=ensure_radius(grp.nodes)) - from meshmode.mesh import Mesh - return Mesh( + return make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True) @@ -857,9 +852,8 @@ def idx(i: int, j: int) -> int: from dataclasses import replace grp = replace(grp, vertex_indices=vertex_indices, nodes=nodes) - from meshmode.mesh import Mesh return ( - Mesh( + make_mesh( vertices, [grp], node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, is_conforming=True), @@ -1007,8 +1001,7 @@ def warp_mesh(mesh: Mesh) -> Mesh: replace(grp, nodes=map_coords(grp.nodes)) for grp in mesh.groups] - from meshmode.mesh import Mesh - return Mesh( + return make_mesh( map_coords(mesh.vertices), groups, node_vertex_consistency_tolerance=False, @@ -1391,8 +1384,7 @@ def generate_box_mesh( # }}} - from meshmode.mesh import Mesh - mesh = Mesh(vertices, [grp], + mesh = make_mesh(vertices, [grp], facial_adjacency_groups=facial_adjacency_groups, is_conforming=True) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 769d6350..8c923bba 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -165,7 +165,7 @@ def get_mesh(self, return_tag_to_elements_map=False): # }}} from meshmode.mesh import ( - Mesh, SimplexElementGroup, TensorProductElementGroup) + SimplexElementGroup, TensorProductElementGroup, make_mesh) bulk_el_types = set() @@ -268,7 +268,7 @@ def get_mesh(self, return_tag_to_elements_map=False): facial_adjacency_groups = _compute_facial_adjacency_from_vertices( groups, np.int32, np.int8, face_vertex_indices_to_tags) - mesh = Mesh( + mesh = make_mesh( vertices, groups, is_conforming=is_conforming, facial_adjacency_groups=facial_adjacency_groups, @@ -374,7 +374,7 @@ def from_meshpy(mesh_info, order=1): which may be generated by either :mod:`meshpy.triangle` or :mod:`meshpy.tet`. """ - from meshmode.mesh import Mesh + from meshmode.mesh import make_mesh from meshmode.mesh.generation import make_group_from_vertices vertices = np.array(mesh_info.points).T @@ -384,7 +384,7 @@ def from_meshpy(mesh_info, order=1): # FIXME: Should transfer boundary/volume markers - return Mesh( + return make_mesh( vertices=vertices, groups=[grp], is_conforming=True) @@ -404,7 +404,7 @@ def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=Fa An array *(nelements, nvertices)* of (mesh-wide) vertex indices. """ - from meshmode.mesh import Mesh + from meshmode.mesh import make_mesh from meshmode.mesh.generation import make_group_from_vertices grp = make_group_from_vertices(vertices, simplices, order) @@ -418,7 +418,7 @@ def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=Fa orient = find_volume_mesh_element_group_orientation(vertices, grp) grp = flip_simplex_element_group(vertices, grp, orient < 0) - return Mesh( + return make_mesh( vertices=vertices, groups=[grp], is_conforming=True) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 88133fe0..1237f3a7 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -35,7 +35,7 @@ from meshmode.mesh import ( BTAG_PARTITION, BoundaryAdjacencyGroup, FacialAdjacencyGroup, InteriorAdjacencyGroup, InterPartAdjacencyGroup, Mesh, MeshElementGroup, PartID, - _FaceIDs) + _FaceIDs, make_mesh) from meshmode.mesh.tools import AffineMap @@ -514,7 +514,7 @@ def _gather_grps(igrp: int) -> List[FacialAdjacencyGroup]: self_facial_adj_groups = [ _gather_grps(igrp) for igrp in range(len(self_mesh_groups))] - return Mesh( + return make_mesh( self_vertices, self_mesh_groups, facial_adjacency_groups=self_facial_adj_groups, @@ -764,7 +764,7 @@ def perform_flips( new_groups.append(new_grp) - return Mesh( + return make_mesh( mesh.vertices, new_groups, skip_tests=skip_tests, is_conforming=mesh.is_conforming, ) @@ -876,7 +876,7 @@ def merge_disjoint_meshes( # }}} - return Mesh( + return make_mesh( vertices, new_groups, skip_tests=skip_tests, nodal_adjacency=nodal_adjacency, @@ -936,7 +936,7 @@ def split_mesh_groups( nodes=grp.nodes[:, mask, :].copy(), )) - mesh = Mesh( + mesh = make_mesh( vertices=mesh.vertices, groups=new_groups, is_conforming=mesh.is_conforming) diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 0c44278a..537090db 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -222,8 +222,8 @@ def refine(self, refine_flags): else: new_vertices = None - from meshmode.mesh import Mesh - new_mesh = Mesh(new_vertices, new_el_groups, is_conforming=( + from meshmode.mesh import make_mesh + new_mesh = make_mesh(new_vertices, new_el_groups, is_conforming=( mesh.is_conforming and (refine_flags.all() or (~refine_flags).all()))) diff --git a/test/test_mesh.py b/test/test_mesh.py index a474d30c..414ccc72 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -41,8 +41,8 @@ from meshmode.discretization.poly_element import ( LegendreGaussLobattoTensorProductGroupFactory, default_simplex_group_factory) from meshmode.mesh import ( - BoundaryAdjacencyGroup, InteriorAdjacencyGroup, Mesh, SimplexElementGroup, - TensorProductElementGroup) + BoundaryAdjacencyGroup, InteriorAdjacencyGroup, SimplexElementGroup, + TensorProductElementGroup, make_mesh) from meshmode.mesh.tools import AffineMap @@ -272,6 +272,7 @@ def test_mesh_as_python(): print(code) exec_dict = {} exec(compile(code, "gen_code.py", "exec"), exec_dict) + exec_dict["_MODULE_SOURCE_CODE"] = code mesh_2 = exec_dict["make_mesh"]() @@ -486,7 +487,7 @@ def test_quad_single_element(visualize=False): np.array([[0, 1, 2, 3]], dtype=np.int32), 30, group_cls=TensorProductElementGroup) - Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) + make_mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) if visualize: import matplotlib.pyplot as plt plt.plot( @@ -980,7 +981,7 @@ def test_quad_mesh_2d(ambient_dim, filename, visualize=False): groups.append(g) - mesh_from_vertices = Mesh(mesh.vertices, groups=groups, is_conforming=True) + mesh_from_vertices = make_mesh(mesh.vertices, groups=groups, is_conforming=True) if visualize: from meshmode.mesh.visualization import write_vertex_vtk_file diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 8415a9b2..2257f679 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -43,7 +43,8 @@ PolynomialWarpAndBlend3DRestrictingGroupFactory, default_simplex_group_factory) from meshmode.dof_array import flat_norm from meshmode.mesh import ( - BTAG_ALL, Mesh, MeshElementGroup, SimplexElementGroup, TensorProductElementGroup) + BTAG_ALL, Mesh, MeshElementGroup, SimplexElementGroup, TensorProductElementGroup, + make_mesh) logger = logging.getLogger(__name__) @@ -553,7 +554,7 @@ def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls, vertex_indices = np.arange(shape.nvertices, dtype=np.int32).reshape(1, -1) mg = group_cls.make_group(mesh_order, vertex_indices, nodes, dim=dim) - mesh = Mesh(vertices, [mg], is_conforming=True) + mesh = make_mesh(vertices, [mg], is_conforming=True) from meshmode.discretization import Discretization vol_discr = Discretization(actx, mesh, group_factory) @@ -645,7 +646,7 @@ def test_sanity_no_elements(actx_factory, dim, mesh_order, group_cls, vertex_indices = np.empty((0, shape.nvertices), dtype=np.int32) mg = group_cls.make_group(mesh_order, vertex_indices, nodes, dim=dim) - mesh = Mesh(vertices, [mg], is_conforming=True) + mesh = make_mesh(vertices, [mg], is_conforming=True) from meshmode.discretization import Discretization vol_discr = Discretization(actx, mesh, group_factory) @@ -869,7 +870,7 @@ def test_mesh_without_vertices(actx_factory): groups = [ replace(grp, nodes=grp.nodes, vertex_indices=None) for grp in mesh.groups] - mesh = Mesh(None, groups, is_conforming=False) + mesh = make_mesh(None, groups, is_conforming=False) # try refining it from meshmode.mesh.refinement import refine_uniformly From 7ed79950ba23e25885927dd89eba808a01854e22 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 Apr 2024 14:28:13 -0500 Subject: [PATCH 5/6] _complete_facial_adjacency_groups: add types, tolerate tuples --- meshmode/mesh/__init__.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 93185d3a..357f4d8c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1689,7 +1689,10 @@ def _merge_boundary_adjacency_groups( def _complete_facial_adjacency_groups( - facial_adjacency_groups, element_id_dtype, face_id_dtype): + facial_adjacency_groups: Tuple[Tuple[FacialAdjacencyGroup, ...], ...], + element_id_dtype: np.dtype, + face_id_dtype: np.dtype + ) -> Tuple[Tuple[FacialAdjacencyGroup, ...], ...]: """ Add :class:`~meshmode.mesh.BoundaryAdjacencyGroup` instances for :class:`~meshmode.mesh.BTAG_NONE`, :class:`~meshmode.mesh.BTAG_ALL`, and @@ -1697,7 +1700,9 @@ def _complete_facial_adjacency_groups( they are not present. """ - completed_facial_adjacency_groups = facial_adjacency_groups.copy() + completed_facial_adjacency_groups = [ + list(fagrps) for fagrps in facial_adjacency_groups + ] for igrp, fagrp_list in enumerate(facial_adjacency_groups): completed_fagrp_list = completed_facial_adjacency_groups[igrp] @@ -1731,7 +1736,9 @@ def _complete_facial_adjacency_groups( igrp, bdry_grps, BTAG_REALLY_ALL, element_id_dtype, face_id_dtype)) - return completed_facial_adjacency_groups + return tuple( + tuple(fagrps) for fagrps in completed_facial_adjacency_groups + ) # }}} From bbc2a466f8093bb6c54210565f90f86091793fda Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 Apr 2024 14:28:37 -0500 Subject: [PATCH 6/6] Add types for _merge_boundary_adjacency_groups --- meshmode/mesh/__init__.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 357f4d8c..39e4021c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1652,7 +1652,12 @@ def _compute_facial_adjacency_from_vertices( # {{{ complete facial adjacency groups def _merge_boundary_adjacency_groups( - igrp, bdry_grps, merged_btag, element_id_dtype, face_id_dtype): + igrp: int, + bdry_grps: Sequence[BoundaryAdjacencyGroup], + merged_btag: BoundaryTag, + element_id_dtype: np.dtype, + face_id_dtype: np.dtype, + ) -> BoundaryAdjacencyGroup: """ Create a new :class:`~meshmode.mesh.BoundaryAdjacencyGroup` containing all of the entries from a list of existing boundary adjacency groups.