hex: add tests for interior facet integration of geometric quantities (#3878)

Author: ksagiyam

tests/firedrake/regression/test_integral_hex.py

@@ -53,3 +53,42 @@ def test_integral_hex_interior_facet_solve(mesh_from_file):
     solve(a == L, sol, bcs=[bc])
     err = assemble((sol - f)**2 * dx)**0.5
     assert err < 1.e-14
+
+
+def make_nonuniform_box_mesh():
+    mesh = BoxMesh(2, 1, 1, 2., 1., 1., hexahedral=True)
+    coordV = mesh.coordinates.function_space()
+    coords = Function(coordV).assign(mesh.coordinates)
+    bc = DirichletBC(coordV.sub(0), 3., 2)
+    bc.apply(coords)
+    return Mesh(coords)
+
+
+@pytest.mark.parametrize('GQ_expected', [(CellSize, sqrt(6.)),
+                                         (CellVolume, 2.),
+                                         (FacetArea, 1.)])
+def test_integral_hex_interior_facet_geometric_quantities(GQ_expected):
+    GQ, expected = GQ_expected
+    mesh = make_nonuniform_box_mesh()
+    x, y, z = SpatialCoordinate(mesh)
+    e = y('+') * z('-')**2
+    E = assemble(e * dS)
+    assert abs(E - 1. / 6.) < 1.e-14
+    a = GQ(mesh)('-')
+    A = assemble(a * dS)
+    assert abs(A - expected) < 1.e-14
+    EA = assemble(e * a * dS)
+    assert abs(EA - E * A) < 1.e-14
+
+
+def test_integral_hex_interior_facet_facet_avg():
+    mesh = make_nonuniform_box_mesh()
+    x, y, z = SpatialCoordinate(mesh)
+    e = y('+') * z('-')**2
+    E = assemble(e * dS)
+    assert abs(E - 1. / 6.) < 1.e-14
+    a = facet_avg(y('-')**3 * z('+')**4)
+    A = assemble(a * dS)
+    assert abs(A - 1. / 20.) < 1.e-14
+    EA = assemble(e * a * dS)
+    assert abs(EA - E * A) < 1.e-14

tsfc/fem.py

@@ -47,15 +47,18 @@
 class ContextBase(ProxyKernelInterface):
     """Common UFL -> GEM translation context."""
 
-    keywords = ('ufl_cell',
-                'fiat_cell',
-                'integral_type',
-                'integration_dim',
-                'entity_ids',
-                'argument_multiindices',
-                'facetarea',
-                'index_cache',
-                'scalar_type')
+    keywords = (
+        'ufl_cell',
+        'fiat_cell',
+        'integral_type',
+        'integration_dim',
+        'entity_ids',
+        'argument_multiindices',
+        'facetarea',
+        'index_cache',
+        'scalar_type',
+        'use_canonical_quadrature_point_ordering',
+    )
 
     def __init__(self, interface, **kwargs):
         ProxyKernelInterface.__init__(self, interface)
@@ -113,6 +116,9 @@ def translator(self):
 
     @cached_property
     def use_canonical_quadrature_point_ordering(self):
+        # Directly set use_canonical_quadrature_point_ordering = False in context
+        # for translation of special nodes, e.g., CellVolume, FacetArea, CellOrigin, and CellVertices,
+        # as quadrature point ordering is not relevant for those node types.
         return isinstance(self.fiat_cell, UFCHexahedron) and self.integral_type in ['exterior_facet', 'interior_facet']
 
 
@@ -162,6 +168,7 @@ def jacobian_at(self, point):
             expr = NegativeRestricted(expr)
         config = {"point_set": PointSingleton(point)}
         config.update(self.config)
+        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
         context = PointSetContext(**config)
         expr = self.preprocess(expr, context)
         return map_expr_dag(context.translator, expr)
@@ -174,6 +181,7 @@ def detJ_at(self, point):
             expr = NegativeRestricted(expr)
         config = {"point_set": PointSingleton(point)}
         config.update(self.config)
+        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
         context = PointSetContext(**config)
         expr = self.preprocess(expr, context)
         return map_expr_dag(context.translator, expr)
@@ -222,6 +230,7 @@ def physical_edge_lengths(self):
         expr = ufl.as_vector([ufl.sqrt(ufl.dot(expr[i, :], expr[i, :])) for i in range(num_edges)])
         config = {"point_set": PointSingleton(cell.make_points(sd, 0, sd+1)[0])}
         config.update(self.config)
+        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
         context = PointSetContext(**config)
         expr = self.preprocess(expr, context)
         return map_expr_dag(context.translator, expr)
@@ -244,6 +253,7 @@ def physical_points(self, point_set, entity=None):
         if entity is not None:
             config.update({name: getattr(self.interface, name)
                            for name in ["integration_dim", "entity_ids"]})
+        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
         context = PointSetContext(**config)
         expr = self.preprocess(expr, context)
         mapped = map_expr_dag(context.translator, expr)
@@ -539,7 +549,7 @@ def translate_cellvolume(terminal, mt, ctx):
 
     config = {name: getattr(ctx, name)
               for name in ["ufl_cell", "index_cache", "scalar_type"]}
-    config.update(interface=interface, quadrature_degree=degree)
+    config.update(interface=interface, quadrature_degree=degree, use_canonical_quadrature_point_ordering=False)
     expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True)
     return expr
 
@@ -553,7 +563,7 @@ def translate_facetarea(terminal, mt, ctx):
     config = {name: getattr(ctx, name)
               for name in ["ufl_cell", "integration_dim", "scalar_type",
                            "entity_ids", "index_cache"]}
-    config.update(interface=ctx, quadrature_degree=degree)
+    config.update(interface=ctx, quadrature_degree=degree, use_canonical_quadrature_point_ordering=False)
     expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True)
     return expr
 
@@ -567,7 +577,7 @@ def translate_cellorigin(terminal, mt, ctx):
 
     config = {name: getattr(ctx, name)
               for name in ["ufl_cell", "index_cache", "scalar_type"]}
-    config.update(interface=ctx, point_set=point_set)
+    config.update(interface=ctx, point_set=point_set, use_canonical_quadrature_point_ordering=False)
     context = PointSetContext(**config)
     return context.translator(expression)
 
@@ -580,7 +590,7 @@ def translate_cell_vertices(terminal, mt, ctx):
 
     config = {name: getattr(ctx, name)
               for name in ["ufl_cell", "index_cache", "scalar_type"]}
-    config.update(interface=ctx, point_set=ps)
+    config.update(interface=ctx, point_set=ps, use_canonical_quadrature_point_ordering=False)
     context = PointSetContext(**config)
     expr = context.translator(ufl_expr)
 

tsfc/kernel_interface/common.py

@@ -215,7 +215,7 @@ def compile_gem(self, ctx):
         # Let the kernel interface inspect the optimised IR to register
         # what kind of external data is required (e.g., cell orientations,
         # cell sizes, etc.).
-        oriented, needs_cell_sizes, tabulations = self.register_requirements(expressions)
+        oriented, needs_cell_sizes, tabulations, need_facet_orientation = self.register_requirements(expressions)
 
         # Extract Variables that are actually used
         active_variables = gem.extract_type(expressions, gem.Variable)
@@ -226,7 +226,7 @@ def compile_gem(self, ctx):
             impero_c = impero_utils.compile_gem(assignments, index_ordering, remove_zeros=True)
         except impero_utils.NoopError:
             impero_c = None
-        return impero_c, oriented, needs_cell_sizes, tabulations, active_variables
+        return impero_c, oriented, needs_cell_sizes, tabulations, active_variables, need_facet_orientation
 
     def fem_config(self):
         """Return a dictionary used with fem.compile_ufl.
@@ -429,6 +429,7 @@ def check_requirements(ir):
     in one pass."""
     cell_orientations = False
     cell_sizes = False
+    facet_orientation = False
     rt_tabs = {}
     for node in traversal(ir):
         if isinstance(node, gem.Variable):
@@ -438,7 +439,9 @@ def check_requirements(ir):
                 cell_sizes = True
             elif node.name.startswith("rt_"):
                 rt_tabs[node.name] = node.shape
-    return cell_orientations, cell_sizes, tuple(sorted(rt_tabs.items()))
+            elif node.name == "facet_orientation":
+                facet_orientation = True
+    return cell_orientations, cell_sizes, tuple(sorted(rt_tabs.items())), facet_orientation
 
 
 def prepare_constant(constant, number):

tsfc/kernel_interface/firedrake_loopy.py

@@ -11,7 +11,7 @@
 
 import loopy as lp
 
-from tsfc import kernel_args, fem
+from tsfc import kernel_args
 from finat.element_factory import create_element
 from tsfc.kernel_interface.common import KernelBuilderBase as _KernelBuilderBase, KernelBuilderMixin, get_index_names, check_requirements, prepare_coefficient, prepare_arguments, prepare_constant
 from tsfc.loopy import generate as generate_loopy
@@ -190,7 +190,7 @@ def set_coefficient_numbers(self, coefficient_numbers):
     def register_requirements(self, ir):
         """Inspect what is referenced by the IR that needs to be
         provided by the kernel interface."""
-        self.oriented, self.cell_sizes, self.tabulations = check_requirements(ir)
+        self.oriented, self.cell_sizes, self.tabulations, _ = check_requirements(ir)
 
     def set_output(self, o):
         """Produce the kernel return argument"""
@@ -368,7 +368,7 @@ def construct_kernel(self, name, ctx, log=False):
         :arg log: bool if the Kernel should be profiled with Log events
         :returns: :class:`Kernel` object
         """
-        impero_c, oriented, needs_cell_sizes, tabulations, active_variables = self.compile_gem(ctx)
+        impero_c, oriented, needs_cell_sizes, tabulations, active_variables, need_facet_orientation = self.compile_gem(ctx)
         if impero_c is None:
             return self.construct_empty_kernel(name)
         info = self.integral_data_info
@@ -418,8 +418,10 @@ def construct_kernel(self, name, ctx, log=False):
         elif info.integral_type in ["interior_facet", "interior_facet_vert"]:
             int_loopy_arg = lp.GlobalArg("facet", numpy.uint32, shape=(2,))
             args.append(kernel_args.InteriorFacetKernelArg(int_loopy_arg))
-        # Will generalise this in the submesh PR.
-        if fem.PointSetContext(**self.fem_config()).use_canonical_quadrature_point_ordering:
+        # The submesh PR will introduce a robust mechanism to check if a Variable
+        # is actually used in the final form of the expression, so there will be
+        # no need to get "need_facet_orientation" from self.compile_gem().
+        if need_facet_orientation:
             if info.integral_type == "exterior_facet":
                 ext_ornt_loopy_arg = lp.GlobalArg("facet_orientation", gem.uint_type, shape=(1,))
                 args.append(kernel_args.ExteriorFacetOrientationKernelArg(ext_ornt_loopy_arg))