diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index de679e68c428c62f536209a1a3b885b3f052c465..f7332a450f593ff094640d87af9556e68754e422 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 2f5633d34025210c9e37e761db1596a0201470ce..22a2a2fd23dc792951a10364c291406b09af11e7 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 @@ -49,85 +50,88 @@ from meshmode.discretization.poly_element import \ # {{{ discretization getters -def get_ellipse_with_ref_mean_curvature(cl_ctx, aspect=1): - nelements = 20 - order = 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) - a = 1 - b = 1/aspect - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) + InterpolatoryQuadratureSimplexGroupFactory(order)) with cl.CommandQueue(cl_ctx) as queue: nodes = discr.nodes().get(queue=queue) + a = 1 + b = 1/aspect t = np.arctan2(nodes[1] * aspect, nodes[0]) 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 = 8 - order = 8 - - from extra_curve_data import unit_square - - mesh = make_curve_mesh( - unit_square, - np.linspace(0, 1, nelements+1), - order) +def get_torus_with_ref_mean_curvature(cl_ctx, h): + order = 4 + r_inner = 1.0 + r_outer = 3.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)) - - return discr, 0 + InterpolatoryQuadratureSimplexGroupFactory(order)) + with cl.CommandQueue(cl_ctx) as queue: + nodes = discr.nodes().get(queue=queue) -def get_unit_sphere_with_ref_mean_curvature(cl_ctx): - order = 8 - - from meshmode.mesh.generation import generate_icosphere - mesh = generate_icosphere(1, order) + # copied from meshmode.mesh.generation.generate_torus + a = r_outer + b = r_inner - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(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 discr, 1 + return discr, (a + 2.0 * b * cosv) / (2 * b * (a + b * cosv)) # }}} # {{{ test_mean_curvature -@pytest.mark.parametrize(("discr_name", "discr_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, - discr_and_ref_mean_curvature_getter): - if discr_name == "unit sphere": - pytest.skip("not implemented in 3D yet") - - import pytential.symbolic.primitives as prim +def test_mean_curvature(ctx_factory, discr_name, resolutions, + discr_and_ref_mean_curvature_getter, visualize=False): ctx = ctx_factory() + queue = cl.CommandQueue(ctx) - discr, ref_mean_curvature = discr_and_ref_mean_curvature_getter(ctx) + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() - with cl.CommandQueue(ctx) as queue: - from pytential import bind + for r in resolutions: + discr, ref_mean_curvature = \ + discr_and_ref_mean_curvature_getter(ctx, r) mean_curvature = bind( discr, - prim.mean_curvature(discr.ambient_dim))(queue) + sym.mean_curvature(discr.ambient_dim))(queue).get(queue) + + h = 1.0 / r + h_error = la.norm(mean_curvature - ref_mean_curvature, np.inf) + eoc.add_data_point(h, h_error) + print(eoc) - assert np.allclose(mean_curvature.get(), ref_mean_curvature) + order = min([g.order for g in discr.groups]) + assert eoc.order_estimate() > order - 1.1 # }}} @@ -174,18 +178,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 ) ]