Implement the switch statement for sumfact facet integrals

python/dune/perftool/pdelab/argument.py

@@ -176,14 +176,14 @@ def name_accumulation_variable(restrictions=None):
             if measure == "cell":
                 restrictions = (Restriction.NONE,)
             else:
-                restrictions = (Restriction.OUTSIDE,)
+                restrictions = (Restriction.NEGATIVE,)
         return name_residual(*restrictions)
     if ft == 'jacobian':
         if restrictions is None:
             if measure == "cell":
                 restrictions = (Restriction.NONE, Restriction.NONE)
             else:
-                restrictions = (Restriction.OUTSIDE, Restriction.OUTSIDE)
+                restrictions = (Restriction.NEGATIVE, Restriction.NEGATIVE)
         return name_jacobian(*restrictions)
     assert False
 

python/dune/perftool/pdelab/geometry.py

@@ -232,6 +232,10 @@ def name_dimension():
     return formdata.geometric_dimension
 
 
+def world_dimension():
+    return name_dimension()
+
+
 def name_intersection_dimension():
     formdata = get_global_context_value('formdata')
     return formdata.geometric_dimension - 1

python/dune/perftool/pdelab/signatures.py

@@ -56,12 +56,12 @@ def assembly_routine_signature():
     formdata = get_global_context_value("formdata")
 
     templates, args = {('residual', 'cell'): (alpha_volume_templates, alpha_volume_args),
-              ('residual', 'exterior_facet'): (alpha_boundary_templates, alpha_boundary_args),
-              ('residual', 'interior_facet'): (alpha_skeleton_templates, alpha_skeleton_args),
-              ('jacobian', 'cell'): (jacobian_volume_templates, jacobian_volume_args),
-              ('jacobian', 'exterior_facet'): (jacobian_boundary_templates, jacobian_boundary_args),
-              ('jacobian', 'interior_facet'): (jacobian_skeleton_templates, jacobian_skeleton_args),
-              }.get((form_type, integral_type), (None, None))
+                       ('residual', 'exterior_facet'): (alpha_boundary_templates, alpha_boundary_args),
+                       ('residual', 'interior_facet'): (alpha_skeleton_templates, alpha_skeleton_args),
+                       ('jacobian', 'cell'): (jacobian_volume_templates, jacobian_volume_args),
+                       ('jacobian', 'exterior_facet'): (jacobian_boundary_templates, jacobian_boundary_args),
+                       ('jacobian', 'interior_facet'): (jacobian_skeleton_templates, jacobian_skeleton_args),
+                       }.get((form_type, integral_type), (None, None))
 
     if templates is None:
         # Check if form is linear
@@ -69,16 +69,45 @@ def assembly_routine_signature():
         linear = is_linear(formdata.original_form)
 
         templates, args = {('jacobian_apply', 'cell', True): (jacobian_apply_volume_templates, jacobian_apply_volume_args),
-                  ('jacobian_apply', 'exterior_facet', True): (jacobian_apply_boundary_templates, jacobian_apply_boundary_args),
-                  ('jacobian_apply', 'interior_facet', True): (jacobian_apply_skeleton_templates, jacobian_apply_skeleton_args),
-                  ('jacobian_apply', 'cell', False): (nonlinear_jacobian_apply_volume_templates, nonlinear_jacobian_apply_volume_args),
-                  ('jacobian_apply', 'exterior_facet', False): (nonlinear_jacobian_apply_boundary_templates, nonlinear_jacobian_apply_boundary_args),
-                  ('jacobian_apply', 'interior_facet', False): (nonlinear_jacobian_apply_skeleton_templates, nonlinear_jacobian_apply_skeleton_args),
-                  }.get((form_type, integral_type, linear), None)
+                           ('jacobian_apply', 'exterior_facet', True): (jacobian_apply_boundary_templates, jacobian_apply_boundary_args),
+                           ('jacobian_apply', 'interior_facet', True): (jacobian_apply_skeleton_templates, jacobian_apply_skeleton_args),
+                           ('jacobian_apply', 'cell', False): (nonlinear_jacobian_apply_volume_templates, nonlinear_jacobian_apply_volume_args),
+                           ('jacobian_apply', 'exterior_facet', False): (nonlinear_jacobian_apply_boundary_templates, nonlinear_jacobian_apply_boundary_args),
+                           ('jacobian_apply', 'interior_facet', False): (nonlinear_jacobian_apply_skeleton_templates, nonlinear_jacobian_apply_skeleton_args),
+                           }.get((form_type, integral_type, linear), None)
 
     return construct_signature(templates(), args(), kernel_name())
 
 
+def assembly_routine_args():
+    integral_type = get_global_context_value("integral_type")
+    form_type = get_global_context_value("form_type")
+    formdata = get_global_context_value("formdata")
+
+    args = {('residual', 'cell'): alpha_volume_args,
+            ('residual', 'exterior_facet'): alpha_boundary_args,
+            ('residual', 'interior_facet'): alpha_skeleton_args,
+            ('jacobian', 'cell'): jacobian_volume_args,
+            ('jacobian', 'exterior_facet'): jacobian_boundary_args,
+            ('jacobian', 'interior_facet'): jacobian_skeleton_args,
+            }.get((form_type, integral_type), None)
+
+    if args is None:
+        # Check if form is linear
+        from dune.perftool.pdelab.driver import is_linear
+        linear = is_linear(formdata.original_form)
+
+        args = {('jacobian_apply', 'cell', True): jacobian_apply_volume_args,
+                ('jacobian_apply', 'exterior_facet', True): jacobian_apply_boundary_args,
+                ('jacobian_apply', 'interior_facet', True): jacobian_apply_skeleton_args,
+                ('jacobian_apply', 'cell', False): nonlinear_jacobian_apply_volume_args,
+                ('jacobian_apply', 'exterior_facet', False): nonlinear_jacobian_apply_boundary_args,
+                ('jacobian_apply', 'interior_facet', False): nonlinear_jacobian_apply_skeleton_args,
+                }.get((form_type, integral_type, linear), None)
+
+    return args()
+
+
 def construct_signature(types, args, name):
     templates = "template<{}>".format(", ".join("typename {}".format(t) for t in set(types)))
     func = "void {}({}) const".format(name, ", ".join("{}{}& {}".format("const " if c else "", t, a) for t, (c, a) in zip(types, args)))

python/dune/perftool/sumfact/__init__.py

@@ -8,6 +8,7 @@ from dune.perftool.sumfact.basis import (lfs_inames,
                                          pymbolic_trialfunction,
                                          pymbolic_trialfunction_gradient,
                                          )
+import dune.perftool.sumfact.switch
 
 from dune.perftool.pdelab import PDELabInterface
 

python/dune/perftool/sumfact/switch.py

+""" boundary and skeleton integrals come in variants in sum factorization - implement the switch! """
+
+from dune.perftool.generation import (backend,
+                                      get_global_context_value,
+                                      global_context,
+                                      )
+from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.pdelab.localoperator import generate_kernel
+from dune.perftool.pdelab.signatures import (assembly_routine_args,
+                                             assembly_routine_signature,
+                                             kernel_name,
+                                             )
+from dune.perftool.cgen.clazz import ClassMember
+
+
+@backend(interface="generate_kernels_per_integral", name="sumfact")
+def generate_kernels_per_integral(integrals):
+    dim = get_global_context_value("formdata").geometric_dimension
+    measure = get_global_context_value("integral_type")
+
+    if measure == "cell":
+        yield generate_kernel(integrals)
+
+    if measure == "exterior_facet":
+        # Generate all necessary kernels
+        for facedir in range(dim):
+            for facemod in range(2):
+                with global_context(facedir_s=facedir, facemod_s=facemod):
+                    yield generate_kernel(integrals)
+
+        # Generate switch statement
+        yield generate_exterior_facet_switch()
+
+    if measure == "interior_facet":
+        # Generate all necessary kernels
+        for facedir_s in range(dim):
+            for facemod_s in range(2):
+                for facedir_n in range(dim):
+                    for facemod_n in range(2):
+                        with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
+                            yield generate_kernel(integrals)
+
+        # Generate switch statement
+        yield generate_interior_facet_switch()
+
+
+def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=None):
+    with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
+        return kernel_name()
+
+
+def generate_exterior_facet_switch():
+    # Extract the signature
+    signature = assembly_routine_signature()
+    args = ", ".join(tuple(a for c, a in assembly_routine_args()))
+    dim = world_dimension()
+
+    # Construct the switch statement
+    block = []
+    block.append("{")
+    block.append("  size_t variant = ig.indexInInside();")
+    block.append("  switch(variant)")
+    block.append("  {")
+
+    for facedir_s in range(dim):
+        for facemod_s in range(2):
+            block.append("    case {}: {}({}); break;".format(dim * facedir_s + facemod_s,
+                                                              get_kernel_name(facedir_s=facedir_s,
+                                                                              facemod_s=facemod_s,
+                                                                              ),
+                                                              args))
+
+    block.append("  }")
+    block.append("}")
+
+    return ClassMember(signature + block)
+
+
+def generate_interior_facet_switch():
+    # Extract the signature
+    signature = assembly_routine_signature()
+    args = ", ".join(tuple(a for c, a in assembly_routine_args()))
+    dim = world_dimension()
+
+    # Construct the switch statement
+    block = []
+    block.append("{")
+    block.append("  size_t variant = ig.indexInOutside() + 6 * ig.indexInInside();")
+    block.append("  switch(variant)")
+    block.append("  {")
+
+    for facedir_s in range(dim):
+        for facemod_s in range(2):
+            for facedir_n in range(dim):
+                for facemod_n in range(2):
+                    block.append("    case {}: {}({}); break;".format((dim * facedir_s + facemod_s) * (2 * dim) + dim * facedir_n + facemod_n,
+                                                                      get_kernel_name(facedir_s=facedir_s,
+                                                                                      facemod_s=facemod_s,
+                                                                                      facedir_n=facedir_n,
+                                                                                      facemod_n=facemod_n,
+                                                                                      ),
+                                                                      args))
+
+    block.append("  }")
+    block.append("}")
+
+    return ClassMember(signature + block)