From bab3926cff44c4b38a3e86a7853832dd7896ed63 Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sun, 2 Jan 2022 14:52:31 -0600
Subject: [PATCH] use generate_sphere to generate a spheroid
gmsh generation was failing in v4.9 and couldn't find a setup
that worked in v4.8 and v4.9. Possible bug?
---
test/mesh_data.py | 45 +++++++++++----------------------------------
test/test_grudge.py | 25 +++++++++++++------------
2 files changed, 24 insertions(+), 46 deletions(-)
diff --git a/test/mesh_data.py b/test/mesh_data.py
index 0671cd87..0ccc369a 100644
--- a/test/mesh_data.py
+++ b/test/mesh_data.py
@@ -61,8 +61,8 @@ class SphereMeshBuilder(MeshBuilder):
radius = 1.0
def get_mesh(self, resolution, mesh_order):
- from meshmode.mesh.generation import generate_icosphere
- return generate_icosphere(self.radius, order=mesh_order,
+ from meshmode.mesh.generation import generate_sphere
+ return generate_sphere(self.radius, order=mesh_order,
uniform_refinement_rounds=resolution)
@@ -70,41 +70,18 @@ class SpheroidMeshBuilder(MeshBuilder):
ambient_dim = 3
mesh_order = 4
- resolutions = [1.0, 0.11, 0.05]
- # resolutions = [1.0, 0.11, 0.05, 0.03, 0.015]
-
- @property
- def radius(self):
- return 0.25
-
- @property
- def diameter(self):
- return 2.0 * self.radius
+ resolutions = [0, 1, 2, 3]
- @property
- def aspect_ratio(self):
- return 2.0
+ radius = 1.0
+ aspect_ratio = 2.0
def get_mesh(self, resolution, mesh_order):
- from meshmode.mesh.io import ScriptSource
- source = ScriptSource("""
- SetFactory("OpenCASCADE");
- Sphere(1) = {{0, 0, 0, {r}}};
- Dilate {{ {{0, 0, 0}}, {{ {r}, {r}, {rr} }} }} {{ Volume{{1}}; }}
- """.format(r=self.diameter, rr=self.aspect_ratio * self.diameter),
- "geo"
- )
-
- from meshmode.mesh.io import generate_gmsh
- mesh = generate_gmsh(source, 2, order=mesh_order,
- other_options=[
- "-optimize_ho",
- "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution
- ],
- target_unit="MM")
-
- from meshmode.mesh.processing import perform_flips
- return perform_flips(mesh, np.ones(mesh.nelements))
+ from meshmode.mesh.generation import generate_sphere
+ mesh = generate_sphere(self.radius, order=mesh_order,
+ uniform_refinement_rounds=resolution)
+
+ from meshmode.mesh.processing import affine_map
+ return affine_map(mesh, A=np.diag([1.0, 1.0, self.aspect_ratio]))
class BoxMeshBuilder(MeshBuilder):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 65f31336..9de657e1 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -33,6 +33,7 @@ pytest_generate_tests = pytest_generate_tests_for_array_contexts(
from arraycontext.container.traversal import thaw
+from meshmode import _acf # noqa: F401
from meshmode.dof_array import flat_norm
import meshmode.mesh.generation as mgen
@@ -509,7 +510,7 @@ def test_2d_gauss_theorem(actx_factory):
assert abs(int_1 - int_2) < 1e-13
-@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
+@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "2-1-spheroid"])
def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
r"""Check the surface divergence theorem.
@@ -531,9 +532,9 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
if mesh_name == "2-1-ellipse":
from mesh_data import EllipseMeshBuilder
builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
- elif mesh_name == "spheroid":
+ elif mesh_name == "2-1-spheroid":
from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
+ builder = SpheroidMeshBuilder(radius=3.1, aspect_ratio=2.0)
elif mesh_name == "circle":
from mesh_data import EllipseMeshBuilder
builder = EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0)
@@ -627,17 +628,17 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
op_global = op.nodal_sum(dcoll, dd, stiff - (stiff_t + kterm))
op_local = op.elementwise_sum(dcoll, dd, stiff - (stiff_t + kterm + flux))
- err_global = abs(op_global)
- err_local = op.norm(dcoll, op_local, np.inf)
- logger.info("errors: global %.5e local %.5e", err_global, err_local)
-
# compute max element size
from grudge.dt_utils import h_max_from_volume
- h_max = h_max_from_volume(dcoll)
+ h_max = actx.to_numpy(h_max_from_volume(dcoll))
+ err_global = actx.to_numpy(abs(op_global))
+ err_local = actx.to_numpy(op.norm(dcoll, op_local, np.inf))
+ logger.info("errors: h_max %.5e global %.5e local %.5e",
+ h_max, err_global, err_local)
- eoc_global.add_data_point(actx.to_numpy(h_max), actx.to_numpy(err_global))
- eoc_local.add_data_point(actx.to_numpy(h_max), actx.to_numpy(err_local))
+ eoc_global.add_data_point(h_max, err_global)
+ eoc_local.add_data_point(h_max, err_local)
if visualize:
from grudge.shortcuts import make_visualizer
@@ -651,8 +652,8 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
# }}}
order = min(builder.order, builder.mesh_order) - 0.5
- logger.info("\n%s", str(eoc_global))
- logger.info("\n%s", str(eoc_local))
+ logger.info("eoc_global:\n%s", str(eoc_global))
+ logger.info("eoc_local:\n%s", str(eoc_local))
assert eoc_global.max_error() < 1.0e-12 \
or eoc_global.order_estimate() > order - 0.5
--
GitLab