diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index f77f9d308233cd1742eb988e7027e529a3ef479e..dd4b7b2641dbc87d7b91d63ec7df77af00571fc5 100644
--- a/python/dune/perftool/compile.py
+++ b/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"):
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index f043f5f1be99cb81cdc66615d40ffd62bcecd957..4f5ad837457446814048704de4fa808b735510d2 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/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
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 861784c6d5f73509ea2fe3df47c2b3008772eceb..222a67548d3ba0d649f157ea29d4694bd244e6f7 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/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
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 439dcf84dee15bb02a9e06d8c376702ef5c81652..a7c28b06164912e5b1d00aaa9f15e458301e9d0c 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/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.
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 7b5d4a2b289cb8cfc72f2f39485b731b849c8c14..9faaa5e9eec073f621d1075dea29c97636164b7d 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/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)
diff --git a/test/stokes/stokes.ufl b/test/stokes/stokes.ufl
index 0af20a0836bae61ea56ab922e3b767e7f14713d3..e84868a14111d0bd4e3b587ecdc97a32577de9b7 100644
--- a/test/stokes/stokes.ufl
+++ b/test/stokes/stokes.ufl
@@ -1,5 +1,8 @@
+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