diff --git a/grudge/op.py b/grudge/op.py
index 3cebfce9ca389b23e90a5fb44839360d5feba8cb..c784537a9ac10b8e4ee06b106135f6dca7e2df54 100644
--- a/grudge/op.py
+++ b/grudge/op.py
@@ -761,37 +761,15 @@ def _apply_inverse_mass_operator(
ref_mass_inverse = reference_inverse_mass_matrix(actx,
element_group=grp)
- # NOTE: Some discretizations can have both affine and non-affine
- # groups. For example, discretizations on hybrid simplex-hex meshes.
- if not grp.is_affine:
- group_data.append(
- # Based on https://arxiv.org/pdf/1608.03836.pdf
- # true_Minv ~ ref_Minv * ref_M * (1/jac_det) * ref_Minv
- actx.einsum("ik,km,em,mj,ej->ei",
- # FIXME: Should we manually create a temporary for
- # the mass inverse?
- ref_mass_inverse,
- reference_mass_matrix(
- actx,
- out_element_group=grp,
- in_element_group=grp
- ),
- jac_inv,
- # FIXME: Should we manually create a temporary for
- # the mass inverse?
- ref_mass_inverse,
- vec_i,
- tagged=(HasElementwiseMatvecTag(),))
- )
- else:
- group_data.append(
- actx.einsum("ij,ej,ej->ei",
- ref_mass_inverse,
- jac_inv,
- vec_i,
- arg_names=("mass_inv_mat", "jac_det_inv", "vec"),
- tagged=(HasElementwiseMatvecTag(),))
- )
+ group_data.append(
+ # Based on https://arxiv.org/pdf/1608.03836.pdf
+ # true_Minv ~ ref_Minv * ref_M * (1/jac_det) * ref_Minv
+ actx.einsum("ei,ij,ej->ei",
+ jac_inv,
+ ref_mass_inverse,
+ vec_i,
+ tagged=(HasElementwiseMatvecTag(),))
+ )
return DOFArray(actx, data=tuple(group_data))
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index 84c11feb883db70b63bcac293c6b672c6e1af605..288cc30cf50b4ebb2dbf7eaac537ef271e6b3ca5 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -591,9 +591,6 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
self.ambient_dim = dcoll.ambient_dim
self.dim = dcoll.dim
- volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
- self.use_wadg = not all(grp.is_affine for grp in volume_discr.groups)
-
map_common_subexpression_uncached = \
IdentityMapper.map_common_subexpression
@@ -639,17 +636,10 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
jac_in * self.rec(expr.field))
elif isinstance(expr.op, op.InverseMassOperator):
- if self.use_wadg:
- # based on https://arxiv.org/pdf/1608.03836.pdf
- return op.RefInverseMassOperator(dd_in, dd_out)(
- op.RefMassOperator(dd_in, dd_out)(
- 1.0/jac_in * op.RefInverseMassOperator(dd_in, dd_out)(
- self.rec(expr.field))
- )
- )
- else:
- return op.RefInverseMassOperator(dd_in, dd_out)(
- 1/jac_in * self.rec(expr.field))
+ # based on https://arxiv.org/pdf/1608.03836.pdf
+ return (
+ 1.0/jac_in * op.RefInverseMassOperator(dd_in, dd_out)(
+ self.rec(expr.field)))
elif isinstance(expr.op, op.FaceMassOperator):
jac_in_surf = sym.area_element(
diff --git a/test/mesh_data.py b/test/mesh_data.py
index b6effe4c0d303f3d0b0f35286bf0e21ef5a926aa..0671cd87d0beac099a3c50625d0cc68d878c60e7 100644
--- a/test/mesh_data.py
+++ b/test/mesh_data.py
@@ -1,4 +1,5 @@
import numpy as np
+import meshmode.mesh.generation as mgen
class MeshBuilder:
@@ -29,8 +30,7 @@ class Curve2DMeshBuilder(MeshBuilder):
resolutions = [16, 32, 64, 128]
def get_mesh(self, resolution, mesh_order):
- from meshmode.mesh.generation import make_curve_mesh
- return make_curve_mesh(
+ return mgen.make_curve_mesh(
self.curve_fn, # pylint: disable=no-member
np.linspace(0.0, 1.0, resolution + 1),
mesh_order)
@@ -42,8 +42,7 @@ class EllipseMeshBuilder(Curve2DMeshBuilder):
@property
def curve_fn(self):
- from meshmode.mesh.generation import ellipse
- return lambda t: self.radius * ellipse(self.aspect_ratio, t)
+ return lambda t: self.radius * mgen.ellipse(self.aspect_ratio, t)
class StarfishMeshBuilder(Curve2DMeshBuilder):
@@ -52,8 +51,7 @@ class StarfishMeshBuilder(Curve2DMeshBuilder):
@property
def curve_fn(self):
- from meshmode.mesh.generation import NArmedStarfish
- return NArmedStarfish(self.narms, self.amplitude)
+ return mgen.NArmedStarfish(self.narms, self.amplitude)
class SphereMeshBuilder(MeshBuilder):
@@ -122,10 +120,18 @@ class BoxMeshBuilder(MeshBuilder):
if not isinstance(resolution, (list, tuple)):
resolution = (resolution,) * self.ambient_dim
- from meshmode.mesh.generation import generate_regular_rect_mesh
- mesh = generate_regular_rect_mesh(
+ return mgen.generate_regular_rect_mesh(
a=self.a, b=self.b,
nelements_per_axis=resolution,
order=mesh_order)
- return mesh
+
+class WarpedRectMeshBuilder(MeshBuilder):
+ resolutions = [4, 6, 8]
+
+ def __init__(self, dim):
+ self.dim = dim
+
+ def get_mesh(self, resolution, mesh_order):
+ return mgen.generate_warped_rect_mesh(
+ dim=self.dim, order=4, nelements_side=6)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index a10611eccfe51bcb95c48c5057b89ce2d369c1a2..4963642ecf46cb300f719686f87b7e63049afd73 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -263,26 +263,35 @@ def test_mass_surface_area(actx_factory, name):
# }}}
-# {{{ mass inverse on surfaces
+# {{{ mass inverse
-@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"])
-def test_surface_mass_operator_inverse(actx_factory, name):
+@pytest.mark.parametrize("name", [
+ "2-1-ellipse",
+ "spheroid",
+ "warped_rect2",
+ "warped_rect3",
+ ])
+def test_mass_operator_inverse(actx_factory, name):
actx = actx_factory()
# {{{ cases
+ import mesh_data
if name == "2-1-ellipse":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
+ # curve
+ builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
elif name == "spheroid":
- from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
+ # surface
+ builder = mesh_data.SpheroidMeshBuilder()
+ elif name.startswith("warped_rect"):
+ builder = mesh_data.WarpedRectMeshBuilder(dim=int(name[-1]))
+
else:
raise ValueError("unknown geometry name: %s" % name)
# }}}
- # {{{ convergence
+ # {{{ inv(m) @ m == id
from pytools.convergence import EOCRecorder
eoc = EOCRecorder()
@@ -316,14 +325,14 @@ def test_surface_mass_operator_inverse(actx_factory, name):
eoc.add_data_point(h_max, inv_error)
- # }}}
-
logger.info("inverse mass error\n%s", str(eoc))
# NOTE: both cases give 1.0e-16-ish at the moment, but just to be on the
# safe side, choose a slightly larger tolerance
assert eoc.max_error() < 1.0e-14
+ # }}}
+
# }}}