Fix neumann boundary poisson and symdiff

python/dune/perftool/pdelab/localoperator.py

@@ -267,10 +267,14 @@ class AccumulationSpace(Record):
             return (self.restriction,)
 
 
-def determine_accumulation_space(expr, number):
+def determine_accumulation_space(expr, number, measure):
     from dune.perftool.ufl.modified_terminals import extract_modified_arguments
     args = extract_modified_arguments(expr, argnumber=number)
 
+    if measure == 'exterior_facet':
+        for ma in args:
+            ma.restriction = Restriction.NEGATIVE
+
     # If this is a residual term we return a dummy object
     if len(args) == 0:
         return AccumulationSpace()
@@ -308,42 +312,44 @@ def determine_accumulation_space(expr, number):
                              )
 
 
-def boundary_predicates():
-    #     # Maybe wrap this instruction into a condiditional. This mostly happens with mixed boundary conditions
-    #     if self.subdomain_id not in ['everywhere', 'otherwise']:
-    #         # We need to reconstruct the subdomain_data parameter of the measure
-    #         # I am *totally* confused as to why this information is not at hand anyway,
-    #         # but conversation with Martin pointed me to dolfin.fem.assembly where this
-    #         # is done in preprocessing with the limitation of only one possible type of
-    #         # modified measure per integral type.
-    #
-    #         # Get the original form and inspect the present measures
-    #         from dune.perftool.generation import get_global_context_value
-    #         original_form = get_global_context_value("formdata").original_form
-    #
-    #         sd = original_form.subdomain_data()
-    #         assert len(sd) == 1
-    #         subdomains, = list(sd.values())
-    #         domain, = list(sd.keys())
-    #         for k in list(subdomains.keys()):
-    #             if subdomains[k] is None:
-    #                     del subdomains[k]
-    #
-    #         # Finally extract the original subdomain_data (which needs to be present!)
-    #         assert self.measure in subdomains
-    #         subdomain_data = subdomains[self.measure]
-    #
-    #         # Determine the name of the parameter function
-    #         name = get_global_context_value("data").object_names[id(subdomain_data)]
-    #
-    #         # Trigger the generation of code for this thing in the parameter class
-    #         from ufl.checks import is_cellwise_constant
-    #         cellwise_constant = is_cellwise_constant(o)
-    #         from dune.perftool.pdelab.parameter import intersection_parameter_function
-    #         intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int')
-    #
-    #         predicates = predicates.union(['{} == {}'.format(name, self.subdomain_id)])
-    return frozenset([])
+def boundary_predicates(expr, measure, subdomain_id):
+    predicates = frozenset([])
+
+    if subdomain_id not in ['everywhere', 'otherwise']:
+        # We need to reconstruct the subdomain_data parameter of the measure
+        # I am *totally* confused as to why this information is not at hand anyway,
+        # but conversation with Martin pointed me to dolfin.fem.assembly where this
+        # is done in preprocessing with the limitation of only one possible type of
+        # modified measure per integral type.
+
+        # Get the original form and inspect the present measures
+        from dune.perftool.generation import get_global_context_value
+        original_form = get_global_context_value("formdata").original_form
+
+        sd = original_form.subdomain_data()
+        assert len(sd) == 1
+        subdomains, = list(sd.values())
+        domain, = list(sd.keys())
+        for k in list(subdomains.keys()):
+            if subdomains[k] is None:
+                    del subdomains[k]
+
+        # Finally extract the original subdomain_data (which needs to be present!)
+        assert measure in subdomains
+        subdomain_data = subdomains[measure]
+
+        # Determine the name of the parameter function
+        name = get_global_context_value("data").object_names[id(subdomain_data)]
+
+        # Trigger the generation of code for this thing in the parameter class
+        from ufl.checks import is_cellwise_constant
+        cellwise_constant = is_cellwise_constant(expr)
+        from dune.perftool.pdelab.parameter import intersection_parameter_function
+        intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int')
+
+        predicates = predicates.union(['{} == {}'.format(name, subdomain_id)])
+
+    return predicates
 
 
 @iname
@@ -355,7 +361,7 @@ def grad_iname(ma):
     return name
 
 
-def generate_accumulation_instruction(visitor, accterm):
+def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
     # First we do the tree traversal to get a pymbolic expression representing this expression
     pymbolic_expr = visitor(accterm.term)
 
@@ -376,14 +382,14 @@ def generate_accumulation_instruction(visitor, accterm):
     pymbolic_expr = Product((pymbolic_expr, test_expr))
 
     # Collect the lfs and lfs indices for the accumulate call
-    test_lfs = determine_accumulation_space(accterm.argument.expr, 0)
+    test_lfs = determine_accumulation_space(accterm.argument.expr, 0, measure)
     # In the jacobian case, also determine the space for the ansatz space
-    ansatz_lfs = determine_accumulation_space(accterm.term, 1)
+    ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure)
 
     from dune.perftool.pdelab.argument import name_accumulation_variable
     accumvar = name_accumulation_variable((ansatz_lfs.get_restriction() + test_lfs.get_restriction()))
 
-    predicates = boundary_predicates()
+    predicates = boundary_predicates(accterm.term, measure, subdomain_id)
 
     rank = 1 if ansatz_lfs.lfs is None else 2
 
@@ -447,7 +453,7 @@ def generate_kernel(integrals):
 
         # Iterate over the terms and generate a kernel
         for term in accterms:
-            generate_accumulation_instruction(visitor, term)
+            generate_accumulation_instruction(visitor, term, measure, subdomain_id)
 
     # Extract the information, which is needed to create a loopy kernel.
     # First extracting it, might be useful to alter it before kernel generation.

python/dune/perftool/ufl/transformations/extract_accumulation_terms.py

@@ -58,7 +58,8 @@ def split_into_accumulation_terms(expr):
         # 4) Further split according to trial function in jacobian terms
         if all_jacobian_args:
             for jac_arg in all_jacobian_args:
-                pass
+                if not isinstance(replace_expr, Zero):
+                    ret.append(AccumulationTerm(replace_expr, test_arg))
         else:
             if not isinstance(replace_expr, Zero):
                 ret.append(AccumulationTerm(replace_expr, test_arg))