Skip to content

Commit

Permalink
Add Gmsh{2,3}DMeshBuilder
Browse files Browse the repository at this point in the history
Also, refactor test MeshBuilder to not pretend to act as source of
truth on mesh order.
  • Loading branch information
inducer committed May 28, 2024
1 parent 52a3bb2 commit 8f8b310
Show file tree
Hide file tree
Showing 4 changed files with 163 additions and 108 deletions.
96 changes: 65 additions & 31 deletions test/mesh_data.py
Original file line number Diff line number Diff line change
@@ -1,44 +1,61 @@
from abc import abstractmethod, ABC
from typing import ClassVar, Hashable, Optional, Sequence
from meshmode.mesh import Mesh
from meshmode.mesh.io import read_gmsh
import numpy as np
import meshmode.mesh.generation as mgen


class MeshBuilder:
order = 4
mesh_order = None
class MeshBuilder(ABC):
resolutions: ClassVar[Sequence[Hashable]]
ambient_dim: ClassVar[int]

def __init__(self, **kwargs):
for k, v in kwargs.items():
setattr(self, k, v)
@abstractmethod
def get_mesh(
self,
resolution: Hashable,
mesh_order: Optional[int] = None
) -> Mesh:
...

if self.mesh_order is None:
self.mesh_order = self.order

@property
def ambient_dim(self):
raise NotImplementedError
class _GmshMeshBuilder(MeshBuilder):
resolutions = [None]

@property
def resolutions(self):
raise NotImplementedError
def __init__(self, filename: str) -> None:
self._mesh_fn = filename

def get_mesh(self, resolution, mesh_order=None) -> Mesh:
assert resolution is None
assert mesh_order is None
return read_gmsh(self._mesh_fn, force_ambient_dim=self.ambient_dim)


class GmshMeshBuilder2D(_GmshMeshBuilder):
ambient_dim = 2

def get_mesh(self, resolution, mesh_order):
raise NotImplementedError

class GmshMeshBuilder3D(_GmshMeshBuilder):
ambient_dim = 3


class Curve2DMeshBuilder(MeshBuilder):
ambient_dim = 2
resolutions = [16, 32, 64, 128]

def get_mesh(self, resolution, mesh_order):
def get_mesh(self, resolution, mesh_order=None):
if mesh_order is None:
mesh_order = 4
return mgen.make_curve_mesh(
self.curve_fn, # pylint: disable=no-member
np.linspace(0.0, 1.0, resolution + 1),
mesh_order)


class EllipseMeshBuilder(Curve2DMeshBuilder):
radius = 3.1
aspect_ratio = 2.0
def __init__(self, radius=3.1, aspect_ratio=2):
self.radius = radius
self.aspect_ratio = aspect_ratio

@property
def curve_fn(self):
Expand All @@ -58,9 +75,13 @@ class SphereMeshBuilder(MeshBuilder):
ambient_dim = 3

resolutions = [0, 1, 2, 3]
radius = 1.0

def get_mesh(self, resolution, mesh_order):
radius: float

def __init__(self, radius=1):
self.radius = radius

def get_mesh(self, resolution, mesh_order=4):
from meshmode.mesh.generation import generate_sphere
return generate_sphere(self.radius, order=mesh_order,
uniform_refinement_rounds=resolution)
Expand All @@ -69,13 +90,16 @@ def get_mesh(self, resolution, mesh_order):
class SpheroidMeshBuilder(MeshBuilder):
ambient_dim = 3

mesh_order = 4
resolutions = [0, 1, 2, 3]

radius = 1.0
aspect_ratio = 2.0
radius: float
aspect_ratio: float

def __init__(self, radius=1, aspect_ratio=2):
self.radius = radius
self.aspect_ratio = aspect_ratio

def get_mesh(self, resolution, mesh_order):
def get_mesh(self, resolution, mesh_order=4):
from meshmode.mesh.generation import generate_sphere
mesh = generate_sphere(self.radius, order=mesh_order,
uniform_refinement_rounds=resolution)
Expand All @@ -84,16 +108,14 @@ def get_mesh(self, resolution, mesh_order):
return affine_map(mesh, A=np.diag([1.0, 1.0, self.aspect_ratio]))


class BoxMeshBuilder(MeshBuilder):
ambient_dim = 2

mesh_order = 1
class _BoxMeshBuilderBase(MeshBuilder):
resolutions = [4, 8, 16]
mesh_order = 1

a = (-0.5, -0.5, -0.5)
b = (+0.5, +0.5, +0.5)

def get_mesh(self, resolution, mesh_order):
def get_mesh(self, resolution, mesh_order=4):
if not isinstance(resolution, (list, tuple)):
resolution = (resolution,) * self.ambient_dim

Expand All @@ -103,12 +125,24 @@ def get_mesh(self, resolution, mesh_order):
order=mesh_order)


class BoxMeshBuilder1D(_BoxMeshBuilderBase):
ambient_dim = 1


class BoxMeshBuilder2D(_BoxMeshBuilderBase):
ambient_dim = 2


class BoxMeshBuilder3D(_BoxMeshBuilderBase):
ambient_dim = 2


class WarpedRectMeshBuilder(MeshBuilder):
resolutions = [4, 6, 8]

def __init__(self, dim):
self.dim = dim

def get_mesh(self, resolution, mesh_order):
def get_mesh(self, resolution, mesh_order=4):
return mgen.generate_warped_rect_mesh(
dim=self.dim, order=4, nelements_side=6)
31 changes: 14 additions & 17 deletions test/test_dt_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,10 @@
PytestPytatoPyOpenCLArrayContextFactory])

from grudge import DiscretizationCollection

import grudge.op as op

import mesh_data

import pytest

import logging
Expand All @@ -55,23 +56,22 @@ def test_geometric_factors_regular_refinement(actx_factory, name):
# {{{ cases

if name == "interval":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=1)
builder = mesh_data.BoxMeshBuilder1D()
elif name == "box2d":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=2)
builder = mesh_data.BoxMeshBuilder2D()
elif name == "box3d":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=3)
builder = mesh_data.BoxMeshBuilder3D()
else:
raise ValueError("unknown geometry name: %s" % name)

# }}}

order = 4

min_factors = []
for resolution in builder.resolutions:
mesh = builder.get_mesh(resolution, builder.mesh_order)
dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
mesh = builder.get_mesh(resolution, order)
dcoll = DiscretizationCollection(actx, mesh, order=order)
min_factors.append(
actx.to_numpy(
op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll))))
Expand All @@ -84,8 +84,8 @@ def test_geometric_factors_regular_refinement(actx_factory, name):
assert np.all(np.isclose(ratios, 2))

# Make sure it works with empty meshes
mesh = builder.get_mesh(0, builder.mesh_order)
dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
mesh = builder.get_mesh(0)
dcoll = DiscretizationCollection(actx, mesh, order=order)
factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841


Expand All @@ -98,14 +98,11 @@ def test_non_geometric_factors(actx_factory, name):
# {{{ cases

if name == "interval":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=1)
builder = mesh_data.BoxMeshBuilder1D()
elif name == "box2d":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=2)
builder = mesh_data.BoxMeshBuilder2D()
elif name == "box3d":
from mesh_data import BoxMeshBuilder
builder = BoxMeshBuilder(ambient_dim=3)
builder = mesh_data.BoxMeshBuilder3D()
else:
raise ValueError("unknown geometry name: %s" % name)

Expand Down
Loading

0 comments on commit 8f8b310

Please sign in to comment.