examples/advection/surface.py

@@ -40,6 +40,7 @@ from pytools.obj_array import make_obj_array
 
 import grudge.dof_desc as dof_desc
 import grudge.op as op
+import grudge.geometry as geo
 
 import logging
 logger = logging.getLogger(__name__)
@@ -237,7 +238,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
 
         df = dof_desc.DOFDesc(FACE_RESTR_INTERIOR)
         face_discr = dcoll.discr_from_dd(df)
-        face_normal = actx.thaw(dcoll.normal(dd=df))
+        face_normal = geo.normal(actx, dcoll, dd=df)
 
         from meshmode.discretization.visualization import make_visualizer
         vis = make_visualizer(actx, face_discr)

examples/geometry.py

@@ -33,6 +33,7 @@ import pyopencl.tools as cl_tools
 from grudge.array_context import PyOpenCLArrayContext
 
 from grudge import DiscretizationCollection, shortcuts
+import grudge.geometry as geo
 
 
 def main(write_output=True):
@@ -52,7 +53,7 @@ def main(write_output=True):
 
     nodes = actx.thaw(dcoll.nodes())
     bdry_nodes = actx.thaw(dcoll.nodes(dd=BTAG_ALL))
-    bdry_normals = actx.thaw(dcoll.normal(dd=BTAG_ALL))
+    bdry_normals = geo.normal(actx, dcoll, dd=BTAG_ALL)
 
     if write_output:
         vis = shortcuts.make_visualizer(dcoll)

examples/hello-grudge.py

@@ -12,6 +12,7 @@ import numpy as np
 import pyopencl as cl
 from grudge.discretization import DiscretizationCollection
 import grudge.op as op
+import grudge.geometry as geo
 from meshmode.mesh.generation import generate_box_mesh
 from meshmode.array_context import PyOpenCLArrayContext
 from grudge.dof_desc import BoundaryDomainTag, FACE_RESTR_INTERIOR
@@ -42,7 +43,7 @@ def left_boundary_condition(x, t):
 def flux(dcoll, u_tpair):
     dd = u_tpair.dd
     velocity = np.array([2 * np.pi])
-    normal = actx.thaw(dcoll.normal(dd))
+    normal = geo.normal(actx, dcoll, dd)
 
     v_dot_n = np.dot(velocity, normal)
     u_upwind = actx.np.where(v_dot_n > 0,

examples/wave/wave-op-mpi.py

@@ -47,6 +47,7 @@ from grudge.discretization import DiscretizationCollection
 from grudge.shortcuts import make_visualizer, compiled_lsrk45_step
 
 import grudge.op as op
+import grudge.geometry as geo
 
 import logging
 logger = logging.getLogger(__name__)
@@ -77,7 +78,7 @@ def wave_flux(actx, dcoll, c, w_tpair):
     v = w_tpair.v
     dd = w_tpair.dd
 
-    normal = actx.thaw(dcoll.normal(dd))
+    normal = geo.normal(actx, dcoll, dd)
 
     flux_weak = WaveState(
         u=v.avg @ normal,

examples/wave/wave-op-var-velocity.py

@@ -40,6 +40,7 @@ from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, DOFDesc
 from grudge.shortcuts import make_visualizer, rk4_step
 
 import grudge.op as op
+import grudge.geometry as geo
 
 import logging
 logger = logging.getLogger(__name__)
@@ -54,7 +55,7 @@ def wave_flux(actx, dcoll, c, w_tpair):
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = actx.thaw(dcoll.normal(dd))
+    normal = geo.normal(actx, dcoll, dd)
 
     flux_weak = flat_obj_array(
             np.dot(v.avg, normal),

grudge/discretization.py

@@ -765,8 +765,15 @@ class DiscretizationCollection:
         :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
         :returns: an object array of frozen :class:`~meshmode.dof_array.DOFArray`\ s.
         """
-        from grudge.geometry import normal
+        warn("DiscretizationCollection.normal is deprecated and "
+             "will stop working in 2024. "
+             "Use grudge.geometry.normal instead. "
+             "DiscretizationCollection.normal may provide non-P^0 normals "
+             "even when this would otherwise be possible, see "
+             "https://github.com/inducer/grudge/issues/314 for details.",
+             DeprecationWarning, stacklevel=2)
 
+        from grudge.geometry import normal
         return self._setup_actx.freeze(normal(self._setup_actx, self, dd))
 
     # }}}

grudge/models/advection.py

@@ -29,6 +29,7 @@ THE SOFTWARE.
 import numpy as np
 import grudge.op as op
 import types
+import grudge.geometry as geo
 
 from grudge.models import HyperbolicOperator
 
@@ -41,7 +42,7 @@ def advection_weak_flux(dcoll, flux_type, u_tpair, velocity):
     """
     actx = u_tpair.int.array_context
     dd = u_tpair.dd
-    normal = actx.thaw(dcoll.normal(dd))
+    normal = geo.normal(actx, dcoll, dd)
     v_dot_n = np.dot(velocity, normal)
 
     flux_type = flux_type.lower()
@@ -90,7 +91,7 @@ class StrongAdvectionOperator(AdvectionOperatorBase):
     def flux(self, u_tpair):
         actx = u_tpair.int.array_context
         dd = u_tpair.dd
-        normal = actx.thaw(self.dcoll.normal(dd))
+        normal = geo.normal(actx, self.dcoll, dd)
         v_dot_normal = np.dot(self.v, normal)
 
         return u_tpair.int * v_dot_normal - self.weak_flux(u_tpair)
@@ -283,7 +284,7 @@ def v_dot_n_tpair(actx, dcoll, velocity, trace_dd):
     from grudge.trace_pair import TracePair
     from meshmode.discretization.connection import FACE_RESTR_INTERIOR
 
-    normal = actx.thaw(dcoll.normal(trace_dd.with_discr_tag(None)))
+    normal = geo.normal(actx, dcoll, trace_dd.with_discr_tag(None))
     v_dot_n = velocity.dot(normal)
     i = op.project(dcoll, trace_dd.with_discr_tag(None), trace_dd, v_dot_n)
 

grudge/models/em.py

@@ -38,6 +38,7 @@ from pytools import memoize_method, levi_civita
 from pytools.obj_array import flat_obj_array, make_obj_array
 
 import grudge.op as op
+import grudge.geometry as geo
 import numpy as np
 
 
@@ -239,7 +240,7 @@ class MaxwellOperator(HyperbolicOperator):
         """
 
         actx = get_container_context_recursively(wtpair)
-        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
+        normal = geo.normal(actx, self.dcoll, wtpair.dd)
 
         if self.fixed_material:
             e, h = self.split_eh(wtpair)
@@ -340,7 +341,7 @@ class MaxwellOperator(HyperbolicOperator):
         """
 
         actx = get_container_context_recursively(w)
-        absorb_normal = actx.thaw(self.dcoll.normal(dd=self.absorb_tag))
+        absorb_normal = geo.normal(actx, self.dcoll, dd=self.absorb_tag)
 
         e, h = self.split_eh(w)
 

grudge/models/euler.py

@@ -64,6 +64,7 @@ from grudge.trace_pair import TracePair
 from pytools.obj_array import make_obj_array
 
 import grudge.op as op
+import grudge.geometry as geo
 
 
 # {{{ Array containers for the Euler model
@@ -203,7 +204,7 @@ class InviscidWallBC(InviscidBCObject):
             state: ConservedEulerField, t=0):
         actx = state.array_context
         dd_base = as_dofdesc("vol").with_discr_tag(DISCR_TAG_BASE)
-        nhat = actx.thaw(dcoll.normal(dd_bc))
+        nhat = geo.normal(actx, dcoll, dd_bc)
         interior = op.project(dcoll, dd_base, dd_bc, state)
 
         return TracePair(
@@ -270,7 +271,7 @@ def euler_numerical_flux(
         exterior=euler_volume_flux(dcoll, q_rr, gamma=gamma)
     )
     num_flux = flux_tpair.avg
-    normal = actx.thaw(dcoll.normal(dd_intfaces))
+    normal = geo.normal(actx, dcoll, dd_intfaces)
 
     if lf_stabilization:
         from arraycontext import outer

grudge/models/wave.py

@@ -35,6 +35,7 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE
 from pytools.obj_array import flat_obj_array
 
 import grudge.op as op
+import grudge.geometry as geo
 
 
 # {{{ constant-velocity
@@ -91,7 +92,7 @@ class WeakWaveOperator(HyperbolicOperator):
         u = wtpair[0]
         v = wtpair[1:]
         actx = u.int.array_context
-        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
+        normal = geo.normal(actx, self.dcoll, wtpair.dd)
 
         central_flux_weak = -self.c*flat_obj_array(
                 np.dot(v.avg, normal),
@@ -134,7 +135,7 @@ class WeakWaveOperator(HyperbolicOperator):
         neu_bc = flat_obj_array(neu_u, -neu_v)
 
         # radiation BCs -------------------------------------------------------
-        rad_normal = actx.thaw(dcoll.normal(dd=self.radiation_tag))
+        rad_normal = geo.normal(actx, dcoll, dd=self.radiation_tag)
 
         rad_u = op.project(dcoll, "vol", self.radiation_tag, u)
         rad_v = op.project(dcoll, "vol", self.radiation_tag, v)
@@ -244,7 +245,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         u = wtpair[1]
         v = wtpair[2:]
         actx = u.int.array_context
-        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
+        normal = geo.normal(actx, self.dcoll, wtpair.dd)
 
         flux_central_weak = -0.5 * flat_obj_array(
             np.dot(v.int*c.int + v.ext*c.ext, normal),
@@ -296,7 +297,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         neu_bc = flat_obj_array(neu_c, neu_u, -neu_v)
 
         # radiation BCs -------------------------------------------------------
-        rad_normal = actx.thaw(dcoll.normal(dd=self.radiation_tag))
+        rad_normal = geo.normal(actx, dcoll, dd=self.radiation_tag)
 
         rad_c = op.project(dcoll, "vol", self.radiation_tag, c)
         rad_u = op.project(dcoll, "vol", self.radiation_tag, u)

test/test_grudge.py

@@ -42,6 +42,7 @@ from grudge import DiscretizationCollection, make_discretization_collection
 
 import grudge.dof_desc as dof_desc
 import grudge.op as op
+import grudge.geometry as geo
 
 
 import pytest
@@ -340,7 +341,7 @@ def test_face_normal_surface(actx_factory, mesh_name):
     )
     surf_normal = surf_normal / actx.np.sqrt(sum(surf_normal**2))
 
-    face_normal_i = actx.thaw(dcoll.normal(df))
+    face_normal_i = geo.normal(actx, dcoll, df)
     face_normal_e = dcoll.opposite_face_connection(
             dof_desc.BoundaryDomainTag(
                 dof_desc.FACE_RESTR_INTERIOR, dof_desc.VTAG_ALL)
@@ -463,7 +464,7 @@ def test_2d_gauss_theorem(actx_factory):
     int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_volm))
 
     prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
-    normal = actx.thaw(dcoll.normal(BTAG_ALL))
+    normal = geo.normal(actx, dcoll, BTAG_ALL)
     int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal))
 
     assert abs(int_1 - int_2) < 1e-13
@@ -572,7 +573,7 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
 
         kappa = summed_curvature(actx, dcoll, dd=dq)
         normal = normal(actx, dcoll, dd=dq)
-        face_normal = actx.thaw(dcoll.normal(df))
+        face_normal = geo.normal(actx, dcoll, df)
         face_f = op.project(dcoll, dd, df, f_num)
 
         # operators
@@ -1064,7 +1065,7 @@ def test_empty_boundary(actx_factory):
             a=(-0.5,)*dim, b=(0.5,)*dim,
             nelements_per_axis=(8,)*dim, order=4)
     dcoll = DiscretizationCollection(actx, mesh, order=4)
-    normal = dcoll.normal(BTAG_NONE)
+    normal = geo.normal(actx, dcoll, BTAG_NONE)
     from meshmode.dof_array import DOFArray
     for component in normal:
         assert isinstance(component, DOFArray)

test/test_op.py

@@ -27,7 +27,7 @@ import meshmode.mesh.generation as mgen
 
 from pytools.obj_array import make_obj_array
 
-from grudge import op, DiscretizationCollection
+from grudge import op, geometry as geo, DiscretizationCollection
 from grudge.dof_desc import DOFDesc
 
 import pytest
@@ -97,7 +97,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested,
         def get_flux(u_tpair):
             dd = u_tpair.dd
             dd_allfaces = dd.with_dtag("all_faces")
-            normal = actx.thaw(dcoll.normal(dd))
+            normal = geo.normal(actx, dcoll, dd)
             u_avg = u_tpair.avg
             if vectorize:
                 if nested:
@@ -223,7 +223,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested,
         def get_flux(u_tpair):
             dd = u_tpair.dd
             dd_allfaces = dd.with_dtag("all_faces")
-            normal = actx.thaw(dcoll.normal(dd))
+            normal = geo.normal(actx, dcoll, dd)
             flux = u_tpair.avg @ normal
             return op.project(dcoll, dd, dd_allfaces, flux)