Fix constraints assembly for systems + lfs index bug in stokes operator

python/dune/perftool/compile.py

@@ -81,14 +81,15 @@ def compile_form():
     from dune.perftool.options import get_option
     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(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(formdata, namedata)
+    from dune.perftool.generation import cache_context, global_context
+    with global_context(formdata=formdata, namedata=namedata):
+        if get_option("driver_file"):
+            with cache_context('driver', delete=True):
+                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(formdata, namedata)
 
         # TODO insert sophisticated analysis/feedback loops here
         if get_option("interactive"):

python/dune/perftool/loopy/transformer.py

@@ -177,6 +177,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
                 from dune.perftool.pdelab.basis import lfs_child
                 lfs = lfs_child(lfs, dimension_iname(count=count, context='arg'))
                 residual_shape[icount] = name_lfs_bound(name_leaf_lfs(subel.sub_elements()[0], ma.restriction))
+                subel = subel.sub_elements()[0]
             else:
                 residual_shape[icount] = name_lfs_bound(lfs)
 
@@ -251,9 +252,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
         element = select_subelement(o.ufl_element(), self.component)
         leaf_element = element
 
-        # Have the issued instruction depend on the iname for this localfunction space
-        self.inames.append(LFSIndexPlaceholder(element, restriction, o.number()))
-
         # Now treat the case of this being a vector finite element
         if element.num_sub_elements() > 0:
             # I cannot handle general mixed elements here...
@@ -272,12 +270,15 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
             # For the purpose of basis evaluation, we need to take the leaf element
             leaf_element = element.sub_elements()[0]
 
+        # Have the issued instruction depend on the iname for this localfunction space
+        self.inames.append(LFSIndexPlaceholder(leaf_element, restriction, o.number()))
+
         if self.grad:
             from dune.perftool.pdelab.argument import name_testfunction_gradient
-            return Subscript(Variable(name_testfunction_gradient(leaf_element, restriction)), (LFSIndexPlaceholder(element, restriction, o.number()),))
+            return Subscript(Variable(name_testfunction_gradient(leaf_element, restriction)), (LFSIndexPlaceholder(leaf_element, restriction, o.number()),))
         else:
             from dune.perftool.pdelab.argument import name_testfunction
-            return Subscript(Variable(name_testfunction(leaf_element, restriction)), (LFSIndexPlaceholder(element, restriction, o.number()),))
+            return Subscript(Variable(name_testfunction(leaf_element, restriction)), (LFSIndexPlaceholder(leaf_element, restriction, o.number()),))
 
     def coefficient(self, o):
         # If this is a trialfunction, we do something entirely different

python/dune/perftool/pdelab/driver.py

@@ -13,6 +13,13 @@ from dune.perftool.generation import (generator_factory,
                                       )
 
 
+def has_constraints(element):
+    from ufl import MixedElement, VectorElement
+    if isinstance(element, MixedElement) and not isinstance(element, VectorElement):
+        return any(has_constraints(el) for el in element.sub_elements())
+    return hasattr(element, "dirichlet_constraints")
+
+
 # Those symbol generators that depend on the meta data we hacked into
 # the ufl.FiniteElement should use fem_metadata_dependent_symbol instead of symbol.
 def get_constraints_metadata(element):
@@ -353,15 +360,18 @@ def define_constraintscontainer(expr, name):
     return ["{} {};".format(cctype, name), "{}.clear();".format(name)]
 
 
-@symbol
-def name_constraintscontainer(element, dirichlet):
-    name = "{}{}_cc".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").lower()
-    define_constraintscontainer(element, name)
+# TODO this does not yet fully support having same elements with different constraints!
+@fem_metadata_dependent_symbol
+def name_constraintscontainer(expr):
+    name = "{}{}_cc".format(FEM_name_mangling(expr), "_dirichlet" if has_constraints(expr) else "").lower()
+    define_constraintscontainer(expr, name)
     return name
 
 
 @preamble
 def define_intersection_lambda(expression, name):
+    if expression is None:
+        return "auto {} = [&](const auto& x){{ return 0; }};".format(name)
     if expression.is_global:
         return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr)
     else:
@@ -369,16 +379,16 @@ def define_intersection_lambda(expression, name):
 
 
 @symbol
-def name_bctype_lambda(dirichlet):
-    define_intersection_lambda(dirichlet, "glambda")
-    # TODO get access to the name data dict here!
-    return "glambda"
+def name_bctype_lambda(name, dirichlet):
+    name = name + "_lambda"
+    define_intersection_lambda(dirichlet, name)
+    return name
 
 
 @preamble
 def define_bctype_function(dirichlet, name):
     gv = name_leafview()
-    bctype_lambda = name_bctype_lambda(dirichlet)
+    bctype_lambda = name_bctype_lambda(name, dirichlet)
     include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
     return "auto {} = Dune::PDELab::makeBoundaryConditionFromCallable({}, {});".format(name,
                                                                                        gv,
@@ -386,30 +396,65 @@ def define_bctype_function(dirichlet, name):
                                                                                        )
 
 
+@preamble
+def define_powergfs_constraints(name, child, dim):
+    include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
+    return "Dune::PDELab::PowerConstraintsParameters<decltype({}), {}> {}({});".format(child, dim, name, child)
+
+
+@preamble
+def define_compositegfs_constraints(name, *children):
+    include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
+    return "Dune::PDELab::CompositeConstraintsParameters<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
+                                                                             name,
+                                                                             ', '.join(c for c in children)
+                                                                             )
+
+
 @symbol
-def name_bctype_function(dirichlet):
-    define_bctype_function(dirichlet, "g")
-    # TODO get access to the name data dict here!
-    return "g"
+def name_bctype_function(expr):
+    from ufl import MixedElement, VectorElement
+    if isinstance(expr, VectorElement):
+        element, (dirichlet, _) = get_constraints_metadata(expr)
+        # get the correct name from the corresponding UFL file
+        from dune.perftool.generation import get_global_context_value
+        name = get_global_context_value("namedata").get(id(dirichlet), "everywhere")
+        define_bctype_function(dirichlet, name)
+        pgfs_name = '{}_{}'.format(name, expr.num_sub_elements())
+        define_powergfs_constraints(pgfs_name, name, expr.num_sub_elements())
+        return pgfs_name
+    if isinstance(expr, MixedElement):
+        children = tuple(name_bctype_function(el) for el in expr.sub_elements())
+        name = '_'.join(children)
+        define_compositegfs_constraints(name, *children)
+        return name
+
+    # So, this seems to be a leaf!
+    element, (dirichlet, _) = get_constraints_metadata(expr)
+
+    # get the correct name from the corresponding UFL file
+    from dune.perftool.generation import get_global_context_value
+    name = get_global_context_value("namedata").get(id(dirichlet), "everywhere")
+
+    define_bctype_function(dirichlet, name)
+    return name
 
 
 @preamble
-def assemble_constraints(name, element, dirichlet):
-    assert dirichlet
-    bctype_function = name_bctype_function(dirichlet)
-    gfs = name_gfs(element)
+def assemble_constraints(name, expr):
+    bctype_function = name_bctype_function(expr)
+    gfs = name_gfs(expr)
     return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
                                                            gfs,
                                                            name,
                                                            )
 
 
-@fem_metadata_dependent_symbol
+@symbol
 def name_assembled_constraints(expr):
-    element, (dirichlet, _) = get_constraints_metadata(expr)
-    name = name_constraintscontainer(element, dirichlet)
-    if dirichlet:
-        assemble_constraints(name, element, dirichlet)
+    name = name_constraintscontainer(expr)
+    if has_constraints(expr):
+        assemble_constraints(name, expr)
     return name
 
 
@@ -790,10 +835,18 @@ def print_matrix():
 
 
 @preamble
-def define_gfs_name(element):
+def define_gfs_name(element, prefix=()):
+    from ufl import MixedElement, VectorElement
+    if isinstance(element, VectorElement):
+        gfs = name_gfs(element)
+        return "{}.name(\"data_{}\");".format(gfs, "_".join(str(p) for p in prefix))
+    if isinstance(element, MixedElement):
+        for (i, el) in enumerate(element.sub_elements()):
+            define_gfs_name(el, prefix + (i,))
+        return ""
+    # This is a scalar leaf
     gfs = name_gfs(element)
-    # TODO make something sensible here
-    return "{}.name(\"bla\");".format(gfs)
+    return "{}.name(\"data_{}\");".format(gfs, "_".join(str(p) for p in prefix))
 
 
 @preamble

python/dune/perftool/pdelab/localoperator.py

@@ -276,62 +276,60 @@ def generate_localoperator_kernels(formdata, namedata):
     # Have a data structure collect the generated kernels
     operator_kernels = {}
 
-    from dune.perftool.generation import global_context
-    with global_context(formdata=formdata, namedata=namedata):
-        with global_context(form_type='residual'):
-            # Generate the necessary residual methods
-            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)
-
-                    enum_pattern()
-                    pattern_baseclass()
-                    enum_alpha()
-                    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(measure)
-                        base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
-
-                        # Add the initializer list for that base class
-                        ini = name_initree_member()
-                        ini_constructor = name_initree_constructor_param()
-                        initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
-                                         ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
-                                         classtag="operator",
-                                         )
+    with global_context(form_type='residual'):
+        # Generate the necessary residual methods
+        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)
+
+                enum_pattern()
+                pattern_baseclass()
+                enum_alpha()
+                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(measure)
+                    base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
+
+                    # Add the initializer list for that base class
+                    ini = name_initree_member()
+                    ini_constructor = name_initree_constructor_param()
+                    initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
+                                     ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
+                                     classtag="operator",
+                                     )
 
                 operator_kernels[(measure, 'residual')] = kernel
 
-        # Generate the necessary jacobian methods
-        if not get_option("numerical_jacobian"):
-            from ufl import derivative
-            from ufl.algorithms import expand_derivatives
-            jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
-
-            with global_context(form_type="jacobian"):
-                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=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
-                # we use empty jacobian assembly methods instead
-                alpha_measures = set(i.integral_type() for i in form.integrals())
-                jacobian_measures = set(i.integral_type() for i in jacform.integrals())
-                for it in alpha_measures - jacobian_measures:
-                    operator_kernels[(it, 'jacobian')] = None
+    # Generate the necessary jacobian methods
+    if not get_option("numerical_jacobian"):
+        from ufl import derivative
+        from ufl.algorithms import expand_derivatives
+        jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+
+        with global_context(form_type="jacobian"):
+            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=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
+            # we use empty jacobian assembly methods instead
+            alpha_measures = set(i.integral_type() for i in form.integrals())
+            jacobian_measures = set(i.integral_type() for i in jacform.integrals())
+            for it in alpha_measures - jacobian_measures:
+                operator_kernels[(it, 'jacobian')] = None
 
         # TODO: JacobianApply for matrix-free computations.
 

python/dune/perftool/ufl/execution.py

@@ -87,3 +87,22 @@ class FiniteElement(ufl.FiniteElement):
 
         # Initialize the original finite element from ufl
         ufl.FiniteElement.__init__(self, *args, **kwargs)
+
+
+class VectorElement(ufl.VectorElement):
+    def __init__(self, *args, **kwargs):
+        if ('dirichlet_constraints' in kwargs) or ('dirichlet_expression' in kwargs):
+            # Get dirichlet_constraints and convert it to Expression if necessary!
+            self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'return true;')
+            if isinstance(self.dirichlet_constraints, str):
+                self.dirichlet_constraints = Expression(self.dirichlet_constraints)
+            assert isinstance(self.dirichlet_constraints, Expression)
+
+            # Get dirichlet_constraints and convert it to Expression if necessary!
+            self.dirichlet_expression = kwargs.pop('dirichlet_expression', 'return 0.0;')
+            if isinstance(self.dirichlet_expression, str):
+                self.dirichlet_expression = Expression(self.dirichlet_expression)
+            assert isinstance(self.dirichlet_expression, Expression)
+
+        # Initialize the original finite element from ufl
+        ufl.VectorElement.__init__(self, *args, **kwargs)

test/stokes/stokes.ufl

+v_bctype = Expression("if (x[0] < 1e-8) return 1; else return 0;", on_intersection=True)
+v_dirichlet = Expression("Dune::FieldVector<double, 2> y(0.0); y[0] = 4*x[1]*(1.-x[1]); return y;")
+
 cell = triangle
-P2 = VectorElement("Lagrange", cell, 2)
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=v_dirichlet)
 P1 = FiniteElement("Lagrange", cell, 1)
 TH = P2 * P1