From fea17e2f17c6a2aa40057a3ef126402ac0eb42e4 Mon Sep 17 00:00:00 2001
From: Thomas Gibson <gibsonthomas1120@hotmail.com>
Date: Mon, 26 Apr 2021 10:20:51 -0500
Subject: [PATCH] Add mass matrix test on surfaces
---
test/test_new_world_grudge.py | 102 +++++++++++++++++++++++++++++++++-
1 file changed, 100 insertions(+), 2 deletions(-)
diff --git a/test/test_new_world_grudge.py b/test/test_new_world_grudge.py
index 97c2b167..b3fe5672 100644
--- a/test/test_new_world_grudge.py
+++ b/test/test_new_world_grudge.py
@@ -35,7 +35,7 @@ import meshmode.mesh.generation as mgen
from pytools.obj_array import flat_obj_array, make_obj_array
import grudge.op as op
-from grudge import DiscretizationCollection
+from grudge import DiscretizationCollection, sym, bind
import grudge.dof_desc as dof_desc
import pytest
@@ -45,7 +45,7 @@ import logging
logger = logging.getLogger(__name__)
-# {{{ mass operator trig integration
+# {{{ mass operator
@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
@pytest.mark.parametrize("quad_tag", [dof_desc.QTAG_NONE, "OVSMP"])
@@ -114,4 +114,102 @@ def test_mass_mat_trig(actx_factory, ambient_dim, quad_tag):
# }}}
+# {{{ mass operator on surface
+
+def _ellipse_surface_area(radius, aspect_ratio):
+ # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html
+ eccentricity = 1.0 - (1/aspect_ratio)**2
+
+ if abs(aspect_ratio - 2.0) < 1.0e-14:
+ # NOTE: hardcoded value so we don't need scipy for the test
+ ellip_e = 1.2110560275684594
+ else:
+ from scipy.special import ellipe # pylint: disable=no-name-in-module
+ ellip_e = ellipe(eccentricity)
+
+ return 4.0 * radius * ellip_e
+
+
+def _spheroid_surface_area(radius, aspect_ratio):
+ # https://en.wikipedia.org/wiki/Ellipsoid#Surface_area
+ a = 1.0
+ c = aspect_ratio
+
+ if a < c:
+ e = np.sqrt(1.0 - (a/c)**2)
+ return 2.0 * np.pi * radius**2 * (1.0 + (c/a) / e * np.arcsin(e))
+ else:
+ e = np.sqrt(1.0 - (c/a)**2)
+ return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e))
+
+
+@pytest.mark.parametrize("name", [
+ "2-1-ellipse", "spheroid", "box2d", "box3d"
+ ])
+def test_mass_surface_area(actx_factory, name):
+ actx = actx_factory()
+
+ # {{{ cases
+
+ if name == "2-1-ellipse":
+ from mesh_data import EllipseMeshBuilder
+ builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
+ surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio)
+ elif name == "spheroid":
+ from mesh_data import SpheroidMeshBuilder
+ builder = SpheroidMeshBuilder()
+ surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio)
+ elif name == "box2d":
+ from mesh_data import BoxMeshBuilder
+ builder = BoxMeshBuilder(ambient_dim=2)
+ surface_area = 1.0
+ elif name == "box3d":
+ from mesh_data import BoxMeshBuilder
+ builder = BoxMeshBuilder(ambient_dim=3)
+ surface_area = 1.0
+ else:
+ raise ValueError("unknown geometry name: %s" % name)
+
+ # }}}
+
+ # {{{ convergence
+
+ from pytools.convergence import EOCRecorder
+ eoc = EOCRecorder()
+
+ for resolution in builder.resolutions:
+ mesh = builder.get_mesh(resolution, builder.mesh_order)
+ dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
+ volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
+
+ logger.info("ndofs: %d", volume_discr.ndofs)
+ logger.info("nelements: %d", volume_discr.mesh.nelements)
+
+ # {{{ compute surface area
+
+ dd = dof_desc.DD_VOLUME
+ ones_volm = volume_discr.zeros(actx) + 1
+ mass_weights = op.mass_operator(dcoll, dd, ones_volm)
+ approx_surface_area = np.dot(actx.to_numpy(flatten(ones_volm)),
+ actx.to_numpy(flatten(mass_weights)))
+
+ logger.info("surface: got {:.5e} / expected {:.5e}".format(
+ approx_surface_area, surface_area))
+ area_error = abs(approx_surface_area - surface_area) / abs(surface_area)
+
+ # }}}
+
+ h_max = bind(dcoll, sym.h_max_from_volume(dcoll.ambient_dim,
+ dim=dcoll.dim, dd=dd))(actx)
+ eoc.add_data_point(h_max, area_error)
+
+ # }}}
+
+ logger.info("surface area error\n%s", str(eoc))
+
+ assert eoc.max_error() < 3e-13 or eoc.order_estimate() > builder.order
+
+# }}}
+
+
# vim: foldmethod=marker
--
GitLab