diff --git a/grudge/geometry/__init__.py b/grudge/geometry/__init__.py
index 02dbc5f155e4d2f0058208b98075e17749222250..4106888f42bf4e2dd5bb94ad90622b3cc4a8955f 100644
--- a/grudge/geometry/__init__.py
+++ b/grudge/geometry/__init__.py
@@ -27,12 +27,14 @@ from grudge.geometry.metrics import (
forward_metric_derivative_mat,
inverse_metric_derivative_mat,
inverse_surface_metric_derivative,
- area_element
+ area_element,
+ normal
)
__all__ = (
"forward_metric_derivative_mat",
"inverse_metric_derivative_mat",
"inverse_surface_metric_derivative",
- "area_element"
+ "area_element",
+ "normal"
)
diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py
index a879c861bf85f447f4b03055ee97d3f59129b5da..b8a109e4075aa0835f5bdb08a27620a3bbc9bddd 100644
--- a/grudge/geometry/metrics.py
+++ b/grudge/geometry/metrics.py
@@ -271,3 +271,60 @@ def area_element(actx, dcoll, dd=None):
return actx.np.sqrt(
pseudoscalar(actx, dcoll, dd=dd).norm_squared()
)
+
+
+def surface_normal(actx, dcoll, dd=None, dim=None):
+ import grudge.dof_desc as dof_desc
+
+ dd = dof_desc.as_dofdesc(dd)
+ dim = dim or dcoll.dim
+
+ # NOTE: Don't be tempted to add a sign here. As it is, it produces
+ # exterior normals for positively oriented curves.
+
+ pder = pseudoscalar(actx, dcoll, dd=dd) / area_element(actx, dcoll, dd=dd)
+
+ # Dorst Section 3.7.2
+ return pder << pder.I.inv()
+
+
+def mv_normal(actx, dcoll, dd):
+ """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`."""
+ import grudge.dof_desc as dof_desc
+
+ dd = dof_desc.as_dofdesc(dd)
+ if not dd.is_trace():
+ raise ValueError("may only request normals on boundaries")
+
+ dim = dcoll.dim
+ ambient_dim = dcoll.ambient_dim
+
+ if dim == ambient_dim - 1:
+ return surface_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"
+ # is (e.g the "normal" to a 1D curve in 3D space is actually a
+ # "normal plane")
+ #
+ # The trick done here is that we take the surface normal, move it to the
+ # face and then take a cross product with the face tangent to get the
+ # correct exterior face normal vector.
+ assert dim == ambient_dim - 2
+
+ from grudge.operators import project
+ import grudge.dof_desc as dof_desc
+
+ volm_normal = project(dcoll, dof_desc.DD_VOLUME, dd,
+ surface_normal(actx, dcoll,
+ dd=dof_desc.DD_VOLUME,
+ dim=dim + 1))
+ pder = pseudoscalar(actx, dcoll, dd=dd)
+
+ mv = -(volm_normal ^ pder) << volm_normal.I.inv()
+
+ return mv / actx.np.sqrt(mv.norm_squared())
+
+
+def normal(actx, dcoll, dd):
+ return mv_normal(actx, dcoll, dd).as_vector()
diff --git a/grudge/operators.py b/grudge/operators.py
index f157810995f1021c8892046bb30a1934d5979217..22d290129316e02346fa2dfb601c6960e46c42a7 100644
--- a/grudge/operators.py
+++ b/grudge/operators.py
@@ -22,9 +22,11 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
+from numbers import Number
from pytools import (
- keyed_memoize_in, memoize_in, memoize_on_first_arg
+ memoize_in,
+ memoize_on_first_arg,
+ # keyed_memoize_in
)
from pytools.obj_array import obj_array_vectorize
from meshmode.array_context import (
@@ -48,6 +50,60 @@ class MassOperatorTag(HasElementwiseMatvecTag):
# }}}
+def project(dcoll, src, tgt, vec):
+ """Project from one discretization to another, e.g. from the
+ volume to the boundary, or from the base to the an overintegrated
+ quadrature discretization.
+
+ :arg src: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one
+ :arg tgt: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ """
+ src = dof_desc.as_dofdesc(src)
+ tgt = dof_desc.as_dofdesc(tgt)
+ if src == tgt:
+ return vec
+
+ if isinstance(vec, np.ndarray):
+ return obj_array_vectorize(
+ lambda el: project(dcoll, src, tgt, el), vec)
+
+ if isinstance(vec, Number):
+ return vec
+
+ return dcoll.connection_from_dds(src, tgt)(vec)
+
+
+# {{{ geometric properties
+
+def nodes(dcoll, dd=None):
+ r"""Return the nodes of a discretization.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ Defaults to the base volume discretization.
+ :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+ """
+ if dd is None:
+ return dcoll._volume_discr.nodes()
+ else:
+ return dcoll.discr_from_dd(dd).nodes()
+
+
+@memoize_on_first_arg
+def normal(dcoll, dd):
+ """Get unit normal to specified surface discretization, *dd*.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
+ :returns: an object array of :class:`~meshmode.dof_array.DOFArray`.
+ """
+ from grudge.geometry import normal
+ surface_discr = dcoll.discr_from_dd(dd)
+ actx = surface_discr._setup_actx
+ return actx.freeze(normal(actx, dcoll, dd))
+
+# }}}
+
+
# {{{ Derivative operator
def reference_derivative_matrices(actx, element_group):
@@ -201,7 +257,7 @@ def reference_stiffness_transpose_matrix(actx, out_element_group, in_element_gro
# in_grp.discretization_key()))
def get_ref_stiffness_transpose_mat(out_grp, in_grp):
if in_grp == out_grp:
- mmat = reference_mass_matrix(actx, in_grp)
+ mmat = reference_mass_matrix(actx, out_grp, in_grp)
return [dmat.T.dot(mmat.T)
for dmat in reference_derivative_matrices(actx, in_grp)]