Towards neumann boundaries: subdomain_data implemented

python/dune/perftool/compile.py

@@ -33,7 +33,8 @@ def read_ufl(uflfile):
     # Extract the form from the given data
     data = interpret_ufl_namespace(namespace)
     form = data.forms[0]
-    form = compute_form_data(form).preprocessed_form
+    formdata = compute_form_data(form)
+    form = formdata.preprocessed_form
 
     # We do not expect more than one form
     assert len(data.forms) == 1
@@ -48,7 +49,7 @@ def read_ufl(uflfile):
     form = transform_form(form, reindexing)
 #    form = transform_form(form, split_arguments)
 
-    return form, data.object_names
+    return formdata, data.object_names
 
 
 def generate_driver(form, filename):
@@ -78,16 +79,16 @@ def generate_driver(form, filename):
 # This function is the entrypoint of the ufl2pdelab executable
 def compile_form():
     from dune.perftool.options import get_option
-    form, namedata = read_ufl(get_option("uflfile"))
+    formdata, namedata = read_ufl(get_option("uflfile"))
 
     from dune.perftool.generation import cache_context
     if get_option("driver_file"):
         with cache_context('driver', delete=True):
-            generate_driver(form, get_option("driver_file"))
+            generate_driver(formdata.preprocessed_form, get_option("driver_file"))
 
     if get_option("operator_file"):
         from dune.perftool.pdelab.localoperator import generate_localoperator_kernels
-        kernels = generate_localoperator_kernels(form, namedata)
+        kernels = generate_localoperator_kernels(formdata, namedata)
 
         # TODO insert sophisticated analysis/feedback loops here
         if get_option("interactive"):

python/dune/perftool/loopy/transformer.py

@@ -85,7 +85,7 @@ def get_count():
     return c
 
 
-def transform_accumulation_term(term):
+def transform_accumulation_term(term, measure, subdomain_id):
     from dune.perftool.ufl.transformations.replace import ReplaceExpression
     from pymbolic.primitives import Variable
 
@@ -180,11 +180,47 @@ def transform_accumulation_term(term):
 
     from dune.perftool.pdelab.quadrature import name_factor
     factor = name_factor()
-    instruction(code="{}.accumulate({}, {}*{});".format(accumvar,
-                                                        ", ".join(accumargs),
-                                                        expr_tv_name,
-                                                        factor,
-                                                        ),
+
+    # Generate the code snippet for this accumulation instruction
+    code = "{}.accumulate({}, {}*{});".format(accumvar,
+                                              ", ".join(accumargs),
+                                              expr_tv_name,
+                                              factor,
+                                              )
+
+    # Maybe wrap this instruction into a condiditional. This mostly happens with mixed boundary conditions
+    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_generic_context_value
+        original_form = get_generic_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!)
+        subdomain_data = subdomains[measure]
+
+        # Determine the name of the parameter function
+        name = get_generic_context_value("namedata")[id(subdomain_data)]
+
+        # Trigger the generation of code for this thing in the parameter class
+        from dune.perftool.pdelab.parameter import parameter_function
+        parameter_function(name, subdomain_data, t='int')
+
+        code = "if ({} == {})\n  {}".format(name, subdomain_id, code)
+
+    # Finally, issue the instruction
+    instruction(code=code,
                 assignees=frozenset({accumvar}),
                 read_variables=frozenset({accumvar, factor, expr_tv_name}),
                 forced_iname_deps=acc_inames,

python/dune/perftool/pdelab/localoperator.py

@@ -135,8 +135,11 @@ def measure_specific_details(measure):
     return ret
 
 
-def generate_kernel(integrand=None, measure=None):
-    assert integrand and measure
+def generate_kernel(integral):
+    integrand = integral.integrand()
+    measure = integral.integral_type()
+    subdomain_id = integral.subdomain_id()
+    subdomain_data = integral.subdomain_data()
 
     # Get the measure specifics
     specifics = measure_specific_details(measure)
@@ -148,7 +151,7 @@ def generate_kernel(integrand=None, measure=None):
     # 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)
+        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.
@@ -204,7 +207,10 @@ def cgen_class_from_cache(tag, members=[]):
     return Class(basename, base_classes=base_classes, members=members + pm, constructors=[constructor], tparam_decls=tparams)
 
 
-def generate_localoperator_kernels(form, namedata):
+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()))
@@ -240,28 +246,30 @@ def generate_localoperator_kernels(form, namedata):
     import functools
     from dune.perftool.generation import generic_context
     namedata_context = functools.partial(generic_context, "namedata")
+    formdata_context = functools.partial(generic_context, "formdata")
 
-    # Generate the necessary residual methods
-    for integral in form.integrals():
-        with namedata_context(namedata):
-            kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
-        operator_kernels[(integral.integral_type(), 'residual')] = kernel
-
-    # Generate the necessary jacobian methods
-    from dune.perftool.options import get_option
-    if get_option("numerical_jacobian"):
-        include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
-    else:
-        from ufl import derivative
-        from ufl.algorithms import expand_derivatives
-        jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
-
-        for integral in jacform.integrals():
+    with formdata_context(formdata):
+        # Generate the necessary residual methods
+        for integral in form.integrals():
             with namedata_context(namedata):
-                kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
-            operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
+                kernel = generate_kernel(integral)
+            operator_kernels[(integral.integral_type(), 'residual')] = kernel
 
-    # TODO: JacobianApply for matrix-free computations.
+        # Generate the necessary jacobian methods
+        from dune.perftool.options import get_option
+        if get_option("numerical_jacobian"):
+            include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
+        else:
+            from ufl import derivative
+            from ufl.algorithms import expand_derivatives
+            jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+
+            for integral in jacform.integrals():
+                with namedata_context(namedata):
+                    kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
+                operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
+
+        # TODO: JacobianApply for matrix-free computations.
 
     # Return the set of generated kernels
     return operator_kernels

python/dune/perftool/pdelab/parameter.py

@@ -39,9 +39,9 @@ def name_paramclass():
 
 
 @class_member("parameterclass", access=AccessModifier.PUBLIC)
-def define_parameter_function_class_member(name, expr):
+def define_parameter_function_class_member(name, expr, t):
     result = ["template<typename E, typename X>",
-              "double {}(const E& e, const X& local) const".format(name),
+              "{} {}(const E& e, const X& local) const".format(t, name),
               "{",
               ]
 
@@ -70,7 +70,7 @@ def evaluate_parameter_function(name):
                                )
 
 
-def parameter_function(name, expr):
+def parameter_function(name, expr, t='double'):
     temporary_variable(name, shape=())
-    define_parameter_function_class_member(name, expr)
+    define_parameter_function_class_member(name, expr, t)
     evaluate_parameter_function(name)

python/dune/perftool/ufl/execution.py

@@ -44,6 +44,8 @@ def TrialFunctions(element):
     return split(TrialFunction(element))
 
 
+_dg0 = FiniteElement("DG", "triangle", 0)
+
 class Expression(Coefficient):
     def __init__(self, expr, is_global=True):
         assert isinstance(expr, str)
@@ -51,7 +53,7 @@ class Expression(Coefficient):
         self.is_global = is_global
 
         # Initialize a coefficient with a dummy finite element map.
-        Coefficient.__init__(self, FiniteElement("DG", "triangle", 0))
+        Coefficient.__init__(self, _dg0)
 
 
 class FiniteElement(ufl.FiniteElement):

test/poisson/CMakeLists.txt

+# 1. Poisson Test Case: source f, bc g + numerical differentiation
 add_generated_executable(UFLFILE poisson.ufl
                          TARGET poisson_numdiff
                          FORM_COMPILER_ARGS --numerical-jacobian
@@ -12,6 +13,7 @@ dune_add_system_test(TARGET poisson_numdiff
 dune_symlink_to_source_files(FILES poisson_ref.vtu)
 
 
+# 2. Poisson Test Case: source f, bc g + symbolic differentiation
 add_generated_executable(UFLFILE poisson.ufl
                          TARGET poisson_symdiff
                          )
@@ -20,3 +22,13 @@ dune_add_system_test(TARGET poisson_symdiff
                      INIFILE poisson.mini
                      SCRIPT dune_vtkcompare.py
                      )
+
+# 3. Poisson Test Case: source f, mixed neumann/dirichlet boundary + numerical differentiation
+add_generated_executable(UFLFILE poisson_neumann.ufl
+                         TARGET poisson_neumann_numdiff
+                         FORM_COMPILER_ARGS --numerical-jacobian
+                         )
+
+dune_add_system_test(TARGET poisson_neumann_numdiff
+                     INIFILE poisson.mini
+                     )

test/poisson/poisson_neumann.ufl

+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
+j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return s*4*x[1]*c.two_norm2()*std::exp(-1.*c.two_norm2());") 
+bctype = Expression("if ((x[1]<1e-8) || (x[1]>1-1e-8)) return 0; else return 1;")
+
+V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g, dirichlet_constraints=bctype)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define the boundary measure that knows where we are...
+ds = ds(subdomain_data=bctype)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx - j*v*ds(0)]