diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index c4c63b99e80a89ccbe515a2cba99f8edfd6845c4..f13d28973a8274f3431af229011ad8fc56fa5bb7 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -274,20 +274,31 @@ class NumReferenceDerivative(DiscretizationProperty): mapper_method = intern("map_num_reference_derivative") -def parametrization_derivative_matrix(ambient_dim, dim, where=None): - """Return a :class:`pymbolic.geometric_algebra.MultiVector` representing - the derivative of the reference-to-global parametrization. +def reference_jacobian(field, field_dim, dim, where=None): + """Return a :class:`np.array` representing the Jacobian of a vector field with + respect to the reference coordinates. """ + jac = np.zeros((field_dim, dim), np.object) - par_grad = np.zeros((ambient_dim, dim), np.object) - for i in range(ambient_dim): + for i in range(field_dim): + field_component = field[i] for j in range(dim): - par_grad[i, j] = NumReferenceDerivative( - frozenset([j]), - NodeCoordinateComponent(i, where), - where) + jac[i, j] = NumReferenceDerivative( + frozenset([j]), + field_component, + where) + + return jac + + +def parametrization_derivative_matrix(ambient_dim, dim, where=None): + """Return a :class:`np.array` representing the derivative of the + reference-to-global parametrization. + """ - return par_grad + return reference_jacobian( + [NodeCoordinateComponent(i, where) for i in range(ambient_dim)], + ambient_dim, dim) def parametrization_derivative(ambient_dim, dim, where=None): @@ -343,9 +354,25 @@ def normal(ambient_dim, dim=None, where=None): cse_scope.DISCRETIZATION) -def mean_curvature(where): - raise NotImplementedError() +def mean_curvature(ambient_dim, dim=None, where=None): + """(Numerical) mean curvature.""" + + 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 = cse( + parametrization_derivative_matrix(ambient_dim, dim, where), + "pd_matrix", cse_scope.DISCRETIZATION) + + xpp, ypp = cse( + reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where), + "p2d_matrix", cse_scope.DISCRETIZATION) + return (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2) # FIXME: make sense of this in the context of GA # def xyz_to_local_matrix(dim, where=None): diff --git a/test/extra_curve_data.py b/test/extra_curve_data.py index e6d667e41ade1e444b97647e58f03f4b6152922b..c6679953f523ad5aac6bb56ed4307c41fec97bf8 100644 --- a/test/extra_curve_data.py +++ b/test/extra_curve_data.py @@ -158,3 +158,10 @@ horseshoe = ( Segment((-5, 1), (0, 1)) + Arc((0, 1), (0.5, 0.5), (0, 0)) ) + +# unit square +unit_square = ( + Segment((1, -1), (1, 1)) + + Segment((1, 1), (-1, 1)) + + Segment((-1, 1), (-1, -1)) + + Segment((-1, -1), (1, -1))) diff --git a/test/test_symbolic.py b/test/test_symbolic.py new file mode 100644 index 0000000000000000000000000000000000000000..83f78d64d2291bd3dfa676e3e0c2a77c16409dc0 --- /dev/null +++ b/test/test_symbolic.py @@ -0,0 +1,137 @@ +from __future__ import division, print_function + +__copyright__ = "Copyright (C) 2017 Matt Wala" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import pytest + +import numpy as np +import pyopencl as cl + +import logging +logger = logging.getLogger(__name__) + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + make_curve_mesh) + +from functools import partial +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + +# {{{ discretization getters + +def get_ellipse_with_ref_mean_curvature(cl_ctx, aspect=1): + nelements = 20 + order = 16 + + 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)) + + with cl.CommandQueue(cl_ctx) as queue: + nodes = discr.nodes().get(queue=queue) + + 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) + + discr = Discretization(cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) + + return discr, 0 + + +def get_unit_sphere_with_ref_mean_curvature(cl_ctx): + order = 8 + + from meshmode.mesh.generation import generate_icosphere + mesh = generate_icosphere(1, order) + + discr = Discretization(cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) + + return discr, 1 + +# }}} + + +@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), + ]) +def test_mean_curvature(ctx_getter, 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_getter() + + discr, ref_mean_curvature = discr_and_ref_mean_curvature_getter(ctx) + + 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) + + +# You can test individual routines by typing +# $ python test_tools.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker