From e27992aa992bcde86f7c01f3704f99c6bb11fbc2 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Fri, 6 May 2016 18:20:49 +0200
Subject: [PATCH] Allow multiple measures of same type -> mixed bc poisson dg

---
 python/dune/perftool/pdelab/localoperator.py | 57 ++++++++++----------
 python/dune/perftool/ufl/execution.py        |  5 ++
 test/poisson/CMakeLists.txt                  |  9 ++++
 3 files changed, 41 insertions(+), 30 deletions(-)

diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 8335f8b0..439dcf84 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -154,20 +154,21 @@ def assembly_routine_signature():
     assert False
 
 
-def generate_kernel(integral):
-    integrand = integral.integrand()
-    measure = integral.integral_type()
-    subdomain_id = integral.subdomain_id()
-    subdomain_data = integral.subdomain_data()
-
-    # Now split the given integrand into accumulation expressions
-    from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
-    accterms = split_into_accumulation_terms(integrand)
-
-    # Iterate over the terms and generate a kernel
-    for term in accterms:
-        from dune.perftool.loopy.transformer import transform_accumulation_term
-        transform_accumulation_term(term, measure, subdomain_id)
+def generate_kernel(integrals):
+    for integral in integrals:
+        integrand = integral.integrand()
+        measure = integral.integral_type()
+        subdomain_id = integral.subdomain_id()
+        subdomain_data = integral.subdomain_data()
+
+        # Now split the given integrand into accumulation expressions
+        from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
+        accterms = split_into_accumulation_terms(integrand)
+
+        # Iterate over the terms and generate a kernel
+        for term in accterms:
+            from dune.perftool.loopy.transformer import transform_accumulation_term
+            transform_accumulation_term(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.
@@ -247,10 +248,6 @@ def generate_localoperator_kernels(formdata, namedata):
     # Extract the relevant attributes of the form data
     form = formdata.preprocessed_form
 
-    # For the moment, I do assume that there is but one integral of each type. This might differ
-    # if you use different quadrature orders for different terms.
-    assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals()))
-
     # Reset the generation cache
     from dune.perftool.generation import delete_cache_items
     delete_cache_items()
@@ -283,23 +280,23 @@ def generate_localoperator_kernels(formdata, namedata):
     with global_context(formdata=formdata, namedata=namedata):
         with global_context(form_type='residual'):
             # Generate the necessary residual methods
-            for integral in form.integrals():
-                # Reset the outer loop
-                from dune.perftool.loopy.transformer import set_outerloop
-                set_outerloop(None)
+            for measure in set(i.integral_type() for i in form.integrals()):
+                with global_context(integral_type=measure):
+                    # Reset the outer loop
+                    from dune.perftool.loopy.transformer import set_outerloop
+                    set_outerloop(None)
 
-                with global_context(integral_type=integral.integral_type()):
                     enum_pattern()
                     pattern_baseclass()
                     enum_alpha()
-                    kernel = generate_kernel(integral)
+                    kernel = generate_kernel(form.integrals_by_type(measure))
 
                     # Maybe add numerical differentiation
                     if get_option("numerical_jacobian"):
                         include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
 
                         _, loptype = class_type_from_cache("operator")
-                        which = ufl_measure_to_pdelab_measure(integral.integral_type())
+                        which = ufl_measure_to_pdelab_measure(measure)
                         base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
 
                         # Add the initializer list for that base class
@@ -310,7 +307,7 @@ def generate_localoperator_kernels(formdata, namedata):
                                          classtag="operator",
                                          )
 
-                operator_kernels[(integral.integral_type(), 'residual')] = kernel
+                operator_kernels[(measure, 'residual')] = kernel
 
         # Generate the necessary jacobian methods
         if not get_option("numerical_jacobian"):
@@ -319,14 +316,14 @@ def generate_localoperator_kernels(formdata, namedata):
             jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
 
             with global_context(form_type="jacobian"):
-                for integral in jacform.integrals():
+                for measure in set(i.integral_type() for i in jacform.integrals()):
                     # Reset the outer loop
                     from dune.perftool.loopy.transformer import set_outerloop
                     set_outerloop(None)
 
-                    with global_context(integral_type=integral.integral_type()):
-                        kernel = generate_kernel(integral)
-                    operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
+                    with global_context(integral_type=measure):
+                        kernel = generate_kernel(jacform.integrals_by_type(measure))
+                    operator_kernels[(measure, 'jacobian')] = kernel
 
                 # Generate dummy functions for those kernels, that vanished in the differentiation process
                 # We *could* solve this problem by using lambda_* terms but we do not really want that, so
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 6002b05f..4f5cf5e2 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -57,6 +57,11 @@ class Expression(Coefficient):
         # Initialize a coefficient with a dummy finite element map.
         Coefficient.__init__(self, _dg0)
 
+    # TODO the subdomain_data code needs a uflid, not idea how to get it here
+    # The standard way through class decorator fails here...
+    def ufl_id(self):
+        return 0
+
 
 class FiniteElement(ufl.FiniteElement):
     def __init__(self, *args, **kwargs):
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index e00a32f4..54b8cfb2 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -77,6 +77,15 @@ dune_add_system_test(TARGET poisson_dg_neumann_numdiff
                      INIFILE poisson_dg_neumann_numdiff.mini
                      )
 
+# 8. Poisson Test Case: DG, mixed bc, symbolic differentiation
+add_generated_executable(UFLFILE poisson_dg_neumann.ufl
+                         TARGET poisson_dg_neumann_symdiff
+                         )
+
+dune_add_system_test(TARGET poisson_dg_neumann_symdiff
+                     INIFILE poisson_dg_neumann_symdiff.mini
+                     )
+
 # Add the reference code with the PDELab localoperator that produced
 # the reference vtk file
 add_executable(poisson_dg_ref reference_main.cc)
-- 
GitLab