From 4d38b905bb3a37c48704a6a7a6186daf6235beeb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 26 May 2021 16:19:22 -0500 Subject: [PATCH] Refactor metrics.normal to handle all normal computations - Rename, document the supporting functions for improved clarity - Drop dim argument from (former) surface_normal --- examples/advection/surface.py | 7 ++----- grudge/geometry/__init__.py | 2 -- grudge/geometry/metrics.py | 39 ++++++++++++++++++++--------------- test/test_grudge.py | 10 ++++----- 4 files changed, 28 insertions(+), 30 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index da0d4666..6603b91e 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -194,12 +194,9 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False): return adv_operator.operator(t, u) # check velocity is tangential - from grudge.geometry import surface_normal + from grudge.geometry import normal - surf_normal = surface_normal( - actx, dcoll, dim=dim-1, - dd=dof_desc.DD_VOLUME - ).as_vector(dtype=object) + surf_normal = normal(actx, dcoll, dd=dof_desc.DD_VOLUME) error = op.norm(dcoll, c.dot(surf_normal), 2) logger.info("u_dot_n: %.5e", error) diff --git a/grudge/geometry/__init__.py b/grudge/geometry/__init__.py index f3d96ead..fec1e595 100644 --- a/grudge/geometry/__init__.py +++ b/grudge/geometry/__init__.py @@ -35,7 +35,6 @@ from grudge.geometry.metrics import ( pseudoscalar, area_element, - surface_normal, normal, second_fundamental_form, @@ -55,7 +54,6 @@ __all__ = ( "pseudoscalar", "area_element", - "surface_normal", "normal", "second_fundamental_form", diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py index 5fc01869..3090cd56 100644 --- a/grudge/geometry/metrics.py +++ b/grudge/geometry/metrics.py @@ -21,7 +21,6 @@ Geometry terms Normal vectors -------------- -.. autofunction:: surface_normal .. autofunction:: normal Curvature tensors @@ -513,21 +512,18 @@ def area_element( # {{{ Surface normal vectors -def surface_normal( - actx: ArrayContext, dcoll: DiscretizationCollection, dim=None, dd=None - ) -> MultiVector: - r"""Computes surface normals at each nodal location. +def rel_mv_normal( + actx: ArrayContext, dcoll: DiscretizationCollection, dd=None) -> MultiVector: + r"""Computes surface normals at each nodal location as a + :class:`~pymbolic.geometric_algebra.MultiVector` relative to the 'volume' + discretization's surface pseudoscalar. - :arg dim: an integer denoting the spatial dimension. Defaults to the - :attr:`grudge.DiscretizationCollection.dim`. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. """ import grudge.dof_desc as dof_desc dd = dof_desc.as_dofdesc(dd) - - if dim is None: - dim = dcoll.discr_from_dd(dd).dim + dim = dcoll.discr_from_dd(dd).dim # NOTE: Don't be tempted to add a sign here. As it is, it produces # exterior normals for positively oriented curves. @@ -539,10 +535,15 @@ def surface_normal( return pder << pder.I.inv() -def _compute_mv_normal( - actx: ArrayContext, dcoll: DiscretizationCollection, dd +def mv_normal( + actx: ArrayContext, dcoll: DiscretizationCollection, dd, ) -> MultiVector: """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`. + This supports both volume discretizations + (where ambient == topological dimension) and surface discretizations + (where ambient == topological dimension + 1). In the latter case, extra + processing ensures that the returned normal is in the local tangent space + of the element at the point where the normal is being evaluated. :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization. :returns: a :class:`~pymbolic.geometric_algebra.MultiVector` @@ -558,7 +559,7 @@ def _compute_mv_normal( ambient_dim = dcoll.ambient_dim if dim == ambient_dim - 1: - return surface_normal(actx, dcoll, dim=dim, dd=dd) + return rel_mv_normal(actx, dcoll, dd=dd) # NOTE: In the case of (d - 2)-dimensional curves, we don't really have # enough information on the face to decide what an "exterior face normal" @@ -575,9 +576,8 @@ def _compute_mv_normal( volm_normal = MultiVector( project(dcoll, dof_desc.DD_VOLUME, dd, - surface_normal( + rel_mv_normal( actx, dcoll, - dim=dim + 1, dd=dof_desc.DD_VOLUME ).as_vector(dtype=object)) ) @@ -590,12 +590,17 @@ def _compute_mv_normal( def normal(actx: ArrayContext, dcoll: DiscretizationCollection, dd): """Get the unit normal to the specified surface discretization, *dd*. + This supports both volume discretizations + (where ambient == topological dimension) and surface discretizations + (where ambient == topological dimension + 1). In the latter case, extra + processing ensures that the returned normal is in the local tangent space + of the element at the point where the normal is being evaluated. :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization. :returns: an object array of :class:`~meshmode.dof_array.DOFArray` containing the unit normals at each nodal location. """ - return _compute_mv_normal(actx, dcoll, dd).as_vector(dtype=object) + return mv_normal(actx, dcoll, dd).as_vector(dtype=object) # }}} @@ -626,7 +631,7 @@ def second_fundamental_form( if dim is None: dim = dcoll.ambient_dim - 1 - normal = surface_normal(actx, dcoll, dim=dim, dd=dd).as_vector(dtype=object) + normal = rel_mv_normal(actx, dcoll, dd=dd).as_vector(dtype=object) if dim == 1: second_ref_axes = [((0, 2),)] diff --git a/test/test_grudge.py b/test/test_grudge.py index 6d219205..913511e1 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -356,7 +356,7 @@ def test_face_normal_surface(actx_factory, mesh_name): # {{{ Compute surface and face normals from meshmode.discretization.connection import FACE_RESTR_INTERIOR - from grudge.geometry import surface_normal + from grudge.geometry import normal dv = dof_desc.DD_VOLUME df = dof_desc.as_dofdesc(FACE_RESTR_INTERIOR) @@ -366,8 +366,7 @@ def test_face_normal_surface(actx_factory, mesh_name): surf_normal = op.project( dcoll, dv, df, - surface_normal(actx, dcoll, - dim=dim, dd=dv).as_vector(dtype=object) + normal(actx, dcoll, dd=dv) ) surf_normal = surf_normal / actx.np.sqrt(sum(surf_normal**2)) @@ -598,11 +597,10 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False): f_num = f(thaw(op.nodes(dcoll, dd=dd), actx)) f_quad_num = f(thaw(op.nodes(dcoll, dd=dq), actx)) - from grudge.geometry import surface_normal, summed_curvature + from grudge.geometry import normal, summed_curvature kappa = summed_curvature(actx, dcoll, dim=dim, dd=dq) - normal = surface_normal(actx, dcoll, - dim=dim, dd=dq).as_vector(dtype=object) + normal = normal(actx, dcoll, dd=dq) face_normal = thaw(op.normal(dcoll, df), actx) face_f = op.project(dcoll, dd, df, f_num) -- GitLab