From cf53a72488ec30b44c936bf8a0d58fe2bae4dd93 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 11 Jun 2017 00:07:03 -0500 Subject: [PATCH 1/7] Symbolic: Implement mean_curvature() for curves in 2D. --- pytential/symbolic/primitives.py | 47 ++++++++--- test/extra_curve_data.py | 7 ++ test/test_symbolic.py | 141 +++++++++++++++++++++++++++++++ 3 files changed, 183 insertions(+), 12 deletions(-) create mode 100644 test/test_symbolic.py diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index c4c63b99..61ae32da 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,21 @@ 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 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 = reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where) + return MultiVector(make_obj_array( + [(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 e6d667e4..c6679953 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 00000000..5e697d1d --- /dev/null +++ b/test/test_symbolic.py @@ -0,0 +1,141 @@ +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) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + 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).as_vector(np.object) + + assert np.allclose(mean_curvature[0].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 -- GitLab From f769d1a135235d60b43059232a357c0b08379ca4 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 11 Jun 2017 00:56:46 -0500 Subject: [PATCH 2/7] mean_curvature(): Wrap derivatives in CSEs. --- pytential/symbolic/primitives.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 61ae32da..0c0d1fec 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -365,10 +365,14 @@ def mean_curvature(ambient_dim, dim=None, where=None): "only know how to calculate curvature for a curve in 2D") xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where) - xpp, ypp = reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where) + xp = cse(xp[0], cse_scope.DISCRETIZATION) + yp = cse(yp[0], cse_scope.DISCRETIZATION) + xpp, ypp = reference_jacobian([xp, yp], ambient_dim, dim, where) + xpp = cse(xpp[0], cse_scope.DISCRETIZATION) + ypp = cse(ypp[0], cse_scope.DISCRETIZATION) return MultiVector(make_obj_array( - [(xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)])) + [(xp*ypp - yp*xpp) / (xp**2 + yp**2)**(3/2)])) # FIXME: make sense of this in the context of GA # def xyz_to_local_matrix(dim, where=None): -- GitLab From 85e934fd31efae6d0341addb1892496a42db1016 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 11 Jun 2017 01:03:26 -0500 Subject: [PATCH 3/7] Revert "mean_curvature(): Wrap derivatives in CSEs." This reverts commit f769d1a135235d60b43059232a357c0b08379ca4. --- pytential/symbolic/primitives.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 0c0d1fec..61ae32da 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -365,14 +365,10 @@ def mean_curvature(ambient_dim, dim=None, where=None): "only know how to calculate curvature for a curve in 2D") xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where) - xp = cse(xp[0], cse_scope.DISCRETIZATION) - yp = cse(yp[0], cse_scope.DISCRETIZATION) - xpp, ypp = reference_jacobian([xp, yp], ambient_dim, dim, where) - xpp = cse(xpp[0], cse_scope.DISCRETIZATION) - ypp = cse(ypp[0], cse_scope.DISCRETIZATION) + xpp, ypp = reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where) return MultiVector(make_obj_array( - [(xp*ypp - yp*xpp) / (xp**2 + yp**2)**(3/2)])) + [(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): -- GitLab From e530297698d10c3ea86ef6b71f6ec13bfd51620b Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 11 Jun 2017 01:05:57 -0500 Subject: [PATCH 4/7] mean_curvature: Wrap matrices in CSEs. --- pytential/symbolic/primitives.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 61ae32da..2102b8a1 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -364,8 +364,13 @@ def mean_curvature(ambient_dim, dim=None, where=None): raise NotImplementedError( "only know how to calculate curvature for a curve in 2D") - xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where) - xpp, ypp = reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where) + 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 MultiVector(make_obj_array( [(xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)])) -- GitLab From 9b29076ba17410df057f0f4e3ee844a44ecb53ee Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 11 Jun 2017 01:51:44 -0500 Subject: [PATCH 5/7] Fix redundant import. --- test/test_symbolic.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 5e697d1d..69e18d8c 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -79,10 +79,6 @@ def get_square_with_ref_mean_curvature(cl_ctx): np.linspace(0, 1, nelements+1), order) - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - discr = Discretization(cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(order)) -- GitLab From 73d30afe1a8234c21e7fbf5064991df736442414 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 14 Jun 2017 22:17:37 -0500 Subject: [PATCH 6/7] mean_curvature(): Return an array. --- pytential/symbolic/primitives.py | 3 +-- test/test_symbolic.py | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 2102b8a1..ef3deb3e 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -372,8 +372,7 @@ def mean_curvature(ambient_dim, dim=None, where=None): reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where), "p2d_matrix", cse_scope.DISCRETIZATION) - return MultiVector(make_obj_array( - [(xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)])) + 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/test_symbolic.py b/test/test_symbolic.py index 69e18d8c..83f78d64 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -118,9 +118,9 @@ def test_mean_curvature(ctx_getter, discr_name, discr_and_ref_mean_curvature_get from pytential import bind mean_curvature = bind( discr, - prim.mean_curvature(discr.ambient_dim))(queue).as_vector(np.object) + prim.mean_curvature(discr.ambient_dim))(queue) - assert np.allclose(mean_curvature[0].get(), ref_mean_curvature) + assert np.allclose(mean_curvature.get(), ref_mean_curvature) # You can test individual routines by typing -- GitLab From 3828c1d4fa531968ff01599bbcedfe2310fa5de1 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 14 Jun 2017 22:18:49 -0500 Subject: [PATCH 7/7] mean_curvature(): Fix check for 2D curves. --- pytential/symbolic/primitives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index ef3deb3e..f13d2897 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -360,7 +360,7 @@ def mean_curvature(ambient_dim, dim=None, where=None): if dim is None: dim = ambient_dim - 1 - if dim != 1 and ambient_dim != 2: + if not (dim == 1 and ambient_dim == 2): raise NotImplementedError( "only know how to calculate curvature for a curve in 2D") -- GitLab