diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 27a0b2450e8c8a58a16322e69518999397263055..872f6bd809ade33a923ab59907b3c58474ecd2d5 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -6,7 +6,30 @@ Currently, these are hardcoded as strings. It would be possible
to switch these to cgen expression. OTOH, there is not much to be
gained there.
"""
-from dune.perftool.generation import include_file, symbol, preamble
+from dune.perftool.generation import (generator_factory,
+ include_file,
+ symbol,
+ preamble,
+ )
+
+
+# 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):
+ """
+ Normally, we suppress the constraints information on a finite element.
+ However, we might want to use it explicitly in those case where it matters
+ This function does the job!
+ """
+ dc = getattr(element, 'dirichlet_constraints', None)
+ de = getattr(element, 'dirichlet_expression', None)
+ return (element, (dc, de))
+
+
+fem_metadata_dependent_symbol = generator_factory(item_tags=("symbol",),
+ cache_key_generator=get_constraints_metadata,
+ )
+
# Have a global variable with the entire form data. This allows functions that depend
# deterministically on the entire data set to directly access it instead of passing it
@@ -267,28 +290,37 @@ def type_orderingtag():
@preamble
-def typedef_constraintsassembler(name):
+def typedef_dirichlet_constraintsassembler(name):
include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
return "typedef Dune::PDELab::ConformingDirichletConstraints {};".format(name)
+@preamble
+def typedef_no_constraintsassembler(name):
+ return "typedef Dune::PDELab::NoConstraints {};".format(name)
+
+
@symbol
-def type_constraintsassembler():
- typedef_constraintsassembler("ConstraintsAssembler")
- return "ConstraintsAssembler"
+def type_constraintsassembler(dirichlet):
+ if dirichlet:
+ typedef_dirichlet_constraintsassembler("DirichletConstraintsAssember")
+ return "DirichletConstraintsAssember"
+ else:
+ typedef_no_constraintsassembler("NoConstraintsAssembler")
+ return "NoConstraintsAssembler"
@preamble
-def typedef_constraintscontainer(expr, name):
- gfs = type_gfs(expr)
+def typedef_constraintscontainer(element, name):
+ gfs = type_gfs(element)
r = type_range()
return "typedef {}::ConstraintsContainer<{}>::Type {};".format(gfs, r, name)
@symbol
-def type_constraintscontainer(expr):
- name = "{}_cc".format(FEM_name_mangling(expr)).upper()
- typedef_constraintscontainer(expr, name)
+def type_constraintscontainer(element):
+ name = "{}_cc".format(FEM_name_mangling(element)).upper()
+ typedef_constraintscontainer(element, name)
return name
@@ -299,20 +331,76 @@ def define_constraintscontainer(expr, name):
@symbol
-def name_constraintscontainer(expr):
- name = "{}_cc".format(FEM_name_mangling(expr)).lower()
- define_constraintscontainer(expr, name)
+def name_constraintscontainer(element, dirichlet):
+ name = "{}{}_cc".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").lower()
+ define_constraintscontainer(element, name)
return name
@preamble
-def typedef_gfs(expr, name):
+def define_cell_lambda(expression, name):
+ param = name_parameters()
+ return "auto {} = [&](const auto& e, const auto& x){{ return {}; }};".format(name, expression)
+
+
+@preamble
+def define_intersection_lambda(expression, name):
+ return "auto {} = [&](const auto& is, const auto& x){{ return {}; }};".format(name, expression)
+
+
+@symbol
+def name_bctype_lambda(dirichlet):
+ define_intersection_lambda(dirichlet, "glambda")
+ # TODO get access to the name data dict here!
+ return "glambda"
+
+
+@preamble
+def define_bctype_function(dirichlet, name):
+ gv = name_leafview()
+ bctype_lambda = name_bctype_lambda(dirichlet)
+ include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
+ return "auto {} = Dune::PDELab::makeBoundaryConditionFromCallable({}, {});".format(name,
+ gv,
+ bctype_lambda,
+ )
+
+
+@symbol
+def name_bctype_function(dirichlet):
+ define_bctype_function(dirichlet, "g")
+ # TODO get access to the name data dict here!
+ return "g"
+
+
+@preamble
+def assemble_constraints(name, element, dirichlet):
+ assert dirichlet
+ bctype_function = name_bctype_function(dirichlet)
+ gfs = name_gfs(element)
+ return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
+ gfs,
+ name,
+ )
+
+
+@fem_metadata_dependent_symbol
+def name_assembled_constraints(expr):
+ element, (dirichlet, _) = get_constraints_metadata(expr)
+ name = name_constraintscontainer(element, dirichlet)
+ if dirichlet:
+ assemble_constraints(name, element, dirichlet)
+ return name
+
+
+@preamble
+def typedef_gfs(element, dirichlet, name):
vb = type_vectorbackend()
from ufl import FiniteElement, MixedElement, VectorElement, EnrichedElement, RestrictedElement
- if isinstance(expr, FiniteElement):
+ if isinstance(element, FiniteElement):
gv = type_leafview()
- fem = type_fem(expr)
- cass = type_constraintsassembler()
+ fem = type_fem(element)
+ cass = type_constraintsassembler(dirichlet)
return "typedef Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}> {};".format(gv, fem, cass, vb, name)
if isinstance(expr, MixedElement):
ot = type_orderingtag()
@@ -323,7 +411,7 @@ def typedef_gfs(expr, name):
gv = type_leafview()
fem = type_fem(expr._sub_elements[0])
dim = name_dimension()
- cass = type_constraintsassembler()
+ cass = type_constraintsassembler(dirichlet)
return "typedef Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}, {}, {}> {};".format(gv, fem, dim, vb, vb, cass, name)
if isinstance(expr, EnrichedElement):
raise NotImplementedError("Dune does not support enriched elements!")
@@ -331,38 +419,37 @@ def typedef_gfs(expr, name):
raise NotImplementedError("Dune does not support restricted elements!")
-@symbol
+@fem_metadata_dependent_symbol
def type_gfs(expr):
- name = "{}_GFS".format(FEM_name_mangling(expr).upper())
- typedef_gfs(expr, name)
+ element, dirichlet = get_constraints_metadata(expr)
+ name = "{}{}_gfs".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").upper()
+ typedef_gfs(element, dirichlet, name)
return name
@preamble
-def define_gfs(expr, name):
- gfstype = type_gfs(expr)
+def define_gfs(element, dirichlet, name):
+ gfstype = type_gfs(element)
from ufl import FiniteElement, MixedElement, VectorElement, EnrichedElement, RestrictedElement
- if isinstance(expr, FiniteElement):
+ if isinstance(element, FiniteElement):
gv = name_leafview()
- fem = name_fem(expr)
+ fem = name_fem(element)
return "{} {}({}, {});".format(gfstype, name, gv, fem)
if isinstance(expr, MixedElement):
- args = ", ".join(name_gfs(childgfs) for childgfs in expr._sub_elements)
+ args = ", ".join(name_gfs(childgfs) for childgfs in element._sub_elements)
return "{} {}({});".format(gfstype, name, args)
if isinstance(expr, VectorElement):
gv = name_leafview()
- fem = name_fem(expr._sub_elements[0])
+ fem = name_fem(element._sub_elements[0])
return "{} {}({}, {});".format(gfstype, name, gv, fem)
- if isinstance(expr, EnrichedElement):
- raise NotImplementedError("Dune does not support enriched elements!")
- if isinstance(expr, RestrictedElement):
- raise NotImplementedError("Dune does not support restricted elements!")
+ raise NotImplementedError("Dune does not support this type of finite element!")
-@symbol
+@fem_metadata_dependent_symbol
def name_gfs(expr):
- name = "{}_gfs".format(FEM_name_mangling(expr)).lower()
- define_gfs(expr, name)
+ element, dirichlet = get_constraints_metadata(expr)
+ name = "{}{}_gfs".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").lower()
+ define_gfs(element, dirichlet, name)
return name
@@ -483,9 +570,9 @@ def type_gridoperator():
def define_gridoperator(name):
gotype = type_gridoperator()
ugfs = name_gfs(_form.coefficients()[0].element())
- ucc = name_constraintscontainer(_form.coefficients()[0].element())
+ ucc = name_assembled_constraints(_form.coefficients()[0].element())
vgfs = name_gfs(_form.arguments()[0].element())
- vcc = name_constraintscontainer(_form.arguments()[0].element())
+ vcc = name_assembled_constraints(_form.arguments()[0].element())
lop = name_localoperator()
mb = name_matrixbackend()
return ["{} {}({}, {}, {}, {}, {}, {});".format(gotype, name, ugfs, ucc, vgfs, vcc, lop, mb),
@@ -518,9 +605,49 @@ def define_vector(name):
return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)]
+@preamble
+def define_boundary_lambda(boundary, name):
+ return "auto {} = [&](const auto& x){{ return {}; }};".format(name, boundary)
+
+
+@symbol
+def name_boundary_lambda(boundary, name):
+ define_boundary_lambda(boundary, name + "lambda")
+ return name + "lambda"
+
+
+@preamble
+def define_boundary_function(boundary, name):
+ gv = name_leafview()
+ lambdaname = name_boundary_lambda(boundary, name)
+ return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
+ gv,
+ lambdaname,
+ )
+
+
+@symbol
+def name_boundary_function(boundary):
+ define_boundary_function(boundary, "b")
+ # TODO add namespace data here
+ return "b"
+
+
+@preamble
+def interpolate_vector(name):
+ define_vector(name)
+ element, (_, boundary) = get_constraints_metadata(_form.coefficients()[0].element())
+ b = name_boundary_function(boundary)
+ gfs = name_gfs(element)
+ return "Dune::PDELab::interpolate({}, {}, {});".format(b,
+ gfs,
+ name,
+ )
+
+
@symbol
def name_vector():
- define_vector("x")
+ interpolate_vector("x")
return "x"
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index a2cd7a143290f0ff461056ef24122e3d0de21ad0..59a8d05f70f1aeb2a0aaf33dc04741fc6c730413 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -4,9 +4,19 @@ UFL.
"""
import ufl
+
from ufl import *
from ufl.split_functions import split
+# Monkey patch to the FiniteElementBase
+# We want to subclass finite elements in order to allow meta data
+# connected to constraints. This is an intermediate solution until
+# we have a proper notion of (global) function spaces in UFL. As is,
+# two finite elements with different meta data would not compare equally,
+# but they should in order to allow hashing generation magic based on
+# them (as done everywhere)
+ufl.FiniteElementBase.__eq__ = lambda s, o: repr(s) == repr(o)
+
class TrialFunction(ufl.Coefficient):
""" A coefficient that always takes the reserved index 0 """
@@ -42,3 +52,11 @@ class Expression(ufl.Coefficient):
# Initialize a coefficient with a dummy finite element map.
ufl.Coefficient.__init__(self, FiniteElement("DG", "triangle", 0))
+
+
+class FiniteElement(ufl.FiniteElement):
+ def __init__(self, *args, **kwargs):
+ if 'dirichlet_constraints' or 'dirichlet_expression' in kwargs:
+ self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'true')
+ self.dirichlet_expression = kwargs.pop('dirichlet_expression', '0.0')
+ ufl.FiniteElement.__init__(self, *args, **kwargs)