From 17809abab7d750cb500e42d7bb518d1724cbb392 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 25 Apr 2019 10:51:41 -0500 Subject: [PATCH 1/5] primitives: implement mean curvature in 3d --- pytential/symbolic/primitives.py | 27 +++++++++++++++++---------- test/test_symbolic.py | 13 ++++++------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index de679e68..f7332a45 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -540,17 +540,24 @@ def mean_curvature(ambient_dim, dim=None, where=None): if dim is None: dim = ambient_dim - 1 - if not (dim == 1 and ambient_dim == 2): - raise NotImplementedError( - "only know how to calculate curvature for a curve in 2D") - - xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where) - - xpp, ypp = cse( - reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where), - "p2d_matrix", cse_scope.DISCRETIZATION) + if ambient_dim == 2 and dim == 1: + # https://en.wikipedia.org/wiki/Curvature#Local_expressions + xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where) + + xpp, ypp = cse( + reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where), + "p2d_matrix", cse_scope.DISCRETIZATION) + + kappa = (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2) + elif ambient_dim == 3 and dim == 2: + # https://en.wikipedia.org/wiki/Mean_curvature#Surfaces_in_3D_space + s_op = shape_operator(ambient_dim, dim=dim, where=where) + kappa = -0.5 * sum(s_op[i, i] for i in range(s_op.shape[0])) + else: + raise NotImplementedError('not available in {}D for {}D surfaces' + .format(ambient_dim, dim)) - return (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2) + return kappa def first_fundamental_form(ambient_dim, dim=None, where=None): diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 2f5633d3..fae9d857 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -90,7 +90,7 @@ def get_square_with_ref_mean_curvature(cl_ctx): def get_unit_sphere_with_ref_mean_curvature(cl_ctx): - order = 8 + order = 16 from meshmode.mesh.generation import generate_icosphere mesh = generate_icosphere(1, order) @@ -113,21 +113,20 @@ def get_unit_sphere_with_ref_mean_curvature(cl_ctx): ]) def test_mean_curvature(ctx_factory, discr_name, discr_and_ref_mean_curvature_getter): - if discr_name == "unit sphere": - pytest.skip("not implemented in 3D yet") - - import pytential.symbolic.primitives as prim ctx = ctx_factory() - discr, ref_mean_curvature = discr_and_ref_mean_curvature_getter(ctx) + import pytential.symbolic.primitives as prim with cl.CommandQueue(ctx) as queue: from pytential import bind mean_curvature = bind( discr, prim.mean_curvature(discr.ambient_dim))(queue) - assert np.allclose(mean_curvature.get(), ref_mean_curvature) + if discr_name == 'unit sphere': + assert np.allclose(mean_curvature.get(), ref_mean_curvature, rtol=1.0e-4) + else: + assert np.allclose(mean_curvature.get(), ref_mean_curvature) # }}} -- GitLab From b10cd238c67a67124d83ae95604cf4fd672a0142 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 1 May 2019 15:57:59 -0500 Subject: [PATCH 2/5] test_symbolic: modify curvature test --- test/test_symbolic.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index fae9d857..d537afef 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -90,15 +90,16 @@ def get_square_with_ref_mean_curvature(cl_ctx): def get_unit_sphere_with_ref_mean_curvature(cl_ctx): - order = 16 + order = 18 + radius = 4.0 from meshmode.mesh.generation import generate_icosphere - mesh = generate_icosphere(1, order) + mesh = generate_icosphere(radius, order) discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) + InterpolatoryQuadratureSimplexGroupFactory(order)) - return discr, 1 + return discr, 1.0 / radius # }}} @@ -123,10 +124,7 @@ def test_mean_curvature(ctx_factory, discr_name, discr, prim.mean_curvature(discr.ambient_dim))(queue) - if discr_name == 'unit sphere': - assert np.allclose(mean_curvature.get(), ref_mean_curvature, rtol=1.0e-4) - else: - assert np.allclose(mean_curvature.get(), ref_mean_curvature) + assert np.allclose(mean_curvature.get(), ref_mean_curvature) # }}} -- GitLab From 41eec6beb2895f1db9ee0a803ecfae0a29969f0c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 2 May 2019 21:09:19 -0500 Subject: [PATCH 3/5] attempt to test mean curvature h-convergence --- test/test_symbolic.py | 76 ++++++++++++++++++++++++------------------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index d537afef..6cc5c388 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -25,6 +25,7 @@ THE SOFTWARE. import pytest import numpy as np +import numpy.linalg as la import pyopencl as cl import pyopencl.array # noqa import pyopencl.clmath # noqa @@ -50,31 +51,29 @@ from meshmode.discretization.poly_element import \ # {{{ discretization getters def get_ellipse_with_ref_mean_curvature(cl_ctx, aspect=1): - nelements = 20 - order = 16 + nelements = 16 + order = 4 mesh = make_curve_mesh( partial(ellipse, aspect), np.linspace(0, 1, nelements+1), order) - a = 1 - b = 1/aspect + def ref_mean_curvature(queue, discr): + a = 1 + b = 1/aspect - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) - - with cl.CommandQueue(cl_ctx) as queue: nodes = discr.nodes().get(queue=queue) + t = np.arctan2(nodes[1] * aspect, nodes[0]) - t = np.arctan2(nodes[1] * aspect, nodes[0]) + return a*b / ((a*np.sin(t))**2 + (b*np.cos(t))**2)**(3/2) - return discr, a*b / ((a*np.sin(t))**2 + (b*np.cos(t))**2)**(3/2) + return mesh, order, ref_mean_curvature def get_square_with_ref_mean_curvature(cl_ctx): - nelements = 8 - order = 8 + nelements = 4 + order = 4 from extra_curve_data import unit_square @@ -83,48 +82,58 @@ def get_square_with_ref_mean_curvature(cl_ctx): np.linspace(0, 1, nelements+1), order) - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) - - return discr, 0 + return mesh, order, lambda queue, discr: 0 def get_unit_sphere_with_ref_mean_curvature(cl_ctx): - order = 18 + order = 8 radius = 4.0 from meshmode.mesh.generation import generate_icosphere mesh = generate_icosphere(radius, order) - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) - return discr, 1.0 / radius + return mesh, order, lambda queue, discr: 1.0 / radius # }}} # {{{ test_mean_curvature -@pytest.mark.parametrize(("discr_name", "discr_and_ref_mean_curvature_getter"), [ +@pytest.mark.parametrize(("discr_name", "mesh_and_ref_mean_curvature_getter"), [ ("unit_circle", get_ellipse_with_ref_mean_curvature), ("2-to-1 ellipse", partial(get_ellipse_with_ref_mean_curvature, aspect=2)), ("square", get_square_with_ref_mean_curvature), ("unit sphere", get_unit_sphere_with_ref_mean_curvature), ]) def test_mean_curvature(ctx_factory, discr_name, - discr_and_ref_mean_curvature_getter): + mesh_and_ref_mean_curvature_getter): ctx = ctx_factory() - discr, ref_mean_curvature = discr_and_ref_mean_curvature_getter(ctx) + queue = cl.CommandQueue(ctx) + + from meshmode.mesh.refinement import RefinerWithoutAdjacency + mesh, order, ref_mean_curvature_fn = mesh_and_ref_mean_curvature_getter(ctx) + refiner = RefinerWithoutAdjacency(mesh) + + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() + + for _ in range(4): + discr = Discretization(ctx, refiner._current_mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) + refiner.refine_uniformly() - import pytential.symbolic.primitives as prim - with cl.CommandQueue(ctx) as queue: - from pytential import bind mean_curvature = bind( discr, - prim.mean_curvature(discr.ambient_dim))(queue) + sym.mean_curvature(discr.ambient_dim))(queue).get(queue) + ref_mean_curvature = ref_mean_curvature_fn(queue, discr) + + h = 1.0 / discr.nnodes + h_error = la.norm(mean_curvature - ref_mean_curvature, np.inf) + + eoc.add_data_point(h, h_error) - assert np.allclose(mean_curvature.get(), ref_mean_curvature) + print(eoc) # }}} @@ -171,18 +180,17 @@ def test_tangential_onb(ctx_factory): def test_expr_pickling(): from sumpy.kernel import LaplaceKernel, AxisTargetDerivative import pickle - import pytential ops_for_testing = [ - pytential.sym.d_dx( + sym.d_dx( 2, - pytential.sym.D( - LaplaceKernel(2), pytential.sym.var("sigma"), qbx_forced_limit=-2 + sym.D( + LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=-2 ) ), - pytential.sym.D( + sym.D( AxisTargetDerivative(0, LaplaceKernel(2)), - pytential.sym.var("sigma"), + sym.var("sigma"), qbx_forced_limit=-2 ) ] -- GitLab From 3ba1b7950cbce199e1b88dfefa1faaf2af39b1ec Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 2 May 2019 22:31:31 -0500 Subject: [PATCH 4/5] mean_curvature: fix h-convergence --- test/test_symbolic.py | 90 +++++++++++++++++++++---------------------- 1 file changed, 44 insertions(+), 46 deletions(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 6cc5c388..1987d8ae 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -50,91 +50,89 @@ from meshmode.discretization.poly_element import \ # {{{ discretization getters -def get_ellipse_with_ref_mean_curvature(cl_ctx, aspect=1): - nelements = 16 +def get_ellipse_with_ref_mean_curvature(cl_ctx, nelements, aspect=1): order = 4 - mesh = make_curve_mesh( partial(ellipse, aspect), np.linspace(0, 1, nelements+1), order) - def ref_mean_curvature(queue, discr): - a = 1 - b = 1/aspect + discr = Discretization(cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) + with cl.CommandQueue(cl_ctx) as queue: nodes = discr.nodes().get(queue=queue) - t = np.arctan2(nodes[1] * aspect, nodes[0]) - return a*b / ((a*np.sin(t))**2 + (b*np.cos(t))**2)**(3/2) + a = 1 + b = 1/aspect + t = np.arctan2(nodes[1] * aspect, nodes[0]) - return mesh, order, ref_mean_curvature + return discr, a*b / ((a*np.sin(t))**2 + (b*np.cos(t))**2)**(3/2) -def get_square_with_ref_mean_curvature(cl_ctx): - nelements = 4 +def get_torus_with_ref_mean_curvature(cl_ctx, h): order = 4 + r_inner = 1.0 + r_outer = 3.0 - from extra_curve_data import unit_square - - mesh = make_curve_mesh( - unit_square, - np.linspace(0, 1, nelements+1), - order) - - return mesh, order, lambda queue, discr: 0 - + from meshmode.mesh.generation import generate_torus + mesh = generate_torus(r_outer, r_inner, + n_outer=h, n_inner=h, order=order) + discr = Discretization(cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) -def get_unit_sphere_with_ref_mean_curvature(cl_ctx): - order = 8 - radius = 4.0 + with cl.CommandQueue(cl_ctx) as queue: + nodes = discr.nodes().get(queue=queue) - from meshmode.mesh.generation import generate_icosphere + # copied from meshmode.mesh.generation.generate_torus + a = r_outer + b = r_inner - mesh = generate_icosphere(radius, order) + u = np.arctan2(nodes[1], nodes[0]) + rvec = np.array([np.cos(u), np.sin(u), np.zeros_like(u)]) + rvec = np.sum(nodes * rvec, axis=0) - a + cosv = np.cos(np.arctan2(nodes[2], rvec)) - return mesh, order, lambda queue, discr: 1.0 / radius + return discr, (a + 2.0 * b * cosv) / (2 * b * (a + b * cosv)) # }}} # {{{ test_mean_curvature -@pytest.mark.parametrize(("discr_name", "mesh_and_ref_mean_curvature_getter"), [ - ("unit_circle", get_ellipse_with_ref_mean_curvature), - ("2-to-1 ellipse", partial(get_ellipse_with_ref_mean_curvature, aspect=2)), - ("square", get_square_with_ref_mean_curvature), - ("unit sphere", get_unit_sphere_with_ref_mean_curvature), +@pytest.mark.parametrize(("discr_name", + "resolutions", + "discr_and_ref_mean_curvature_getter"), [ + ("unit_circle", [16, 32, 64], + get_ellipse_with_ref_mean_curvature), + ("2-to-1 ellipse", [16, 32, 64], + partial(get_ellipse_with_ref_mean_curvature, aspect=2)), + ("torus", [8, 10, 12, 16], + get_torus_with_ref_mean_curvature), ]) -def test_mean_curvature(ctx_factory, discr_name, - mesh_and_ref_mean_curvature_getter): +def test_mean_curvature(ctx_factory, discr_name, resolutions, + discr_and_ref_mean_curvature_getter, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - from meshmode.mesh.refinement import RefinerWithoutAdjacency - mesh, order, ref_mean_curvature_fn = mesh_and_ref_mean_curvature_getter(ctx) - refiner = RefinerWithoutAdjacency(mesh) - from pytools.convergence import EOCRecorder eoc = EOCRecorder() - for _ in range(4): - discr = Discretization(ctx, refiner._current_mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) - refiner.refine_uniformly() - + for r in resolutions: + discr, ref_mean_curvature = \ + discr_and_ref_mean_curvature_getter(ctx, r) mean_curvature = bind( discr, sym.mean_curvature(discr.ambient_dim))(queue).get(queue) - ref_mean_curvature = ref_mean_curvature_fn(queue, discr) - h = 1.0 / discr.nnodes + h = 1.0 / r h_error = la.norm(mean_curvature - ref_mean_curvature, np.inf) - eoc.add_data_point(h, h_error) - print(eoc) + order = min([g.order for g in discr.groups]) + assert eoc.order_estimate() > order - 0.75 + # }}} -- GitLab From 83c3d1570bead0703b11f7aeacbc348ebda2e07c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 3 May 2019 15:34:20 -0500 Subject: [PATCH 5/5] relax order assertion --- test/test_symbolic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 1987d8ae..22a2a2fd 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -131,7 +131,7 @@ def test_mean_curvature(ctx_factory, discr_name, resolutions, print(eoc) order = min([g.order for g in discr.groups]) - assert eoc.order_estimate() > order - 0.75 + assert eoc.order_estimate() > order - 1.1 # }}} -- GitLab