Also handle boundary type lambdas

Need to use handler from pymbolic as loopy codegen handlers
are too involved for standalone usage.

python/dune/perftool/pdelab/driver.py

@@ -430,12 +430,27 @@ def name_constraintscontainer(expr):
 
 @preamble
 def define_intersection_lambda(expression, name):
+    from dune.perftool.ufl.execution import Expression
+    from ufl.classes import Expr
     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[0])
-    else:
-        return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
+    elif isinstance(expression, Expression):
+        if expression.is_global:
+            return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr[0])
+        else:
+            return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
+    elif isinstance(expression, Expr):
+        # Set up a visitor
+        with global_context(integral_type="exterior_facet", formdata=_driver_data["formdata"], driver=True):
+            from dune.perftool.ufl.visitor import UFL2LoopyVisitor
+            from dune.perftool.pdelab import PDELabInterface
+            visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
+            from pymbolic.mapper.c_code import CCodeMapper
+            ccm = CCodeMapper()
+            expr = visitor(expression)
+            return "auto {} = [&](const auto& x){{ return {}; }};".format(name, ccm(expr))
+
+    raise ValueError("Expression not understood")
 
 
 def name_bctype_lambda(name, dirichlet):

python/dune/perftool/pdelab/localoperator.py

@@ -30,6 +30,7 @@ from dune.perftool.ufl.modified_terminals import Restriction
 import dune.perftool.loopy.mangler
 
 from pymbolic.primitives import Variable
+import pymbolic.primitives as prim
 from pytools import Record
 
 import loopy as lp
@@ -250,7 +251,7 @@ def determine_accumulation_space(expr, number, measure):
                              )
 
 
-def boundary_predicates(expr, measure, subdomain_id):
+def boundary_predicates(expr, measure, subdomain_id, visitor):
     predicates = frozenset([])
 
     if subdomain_id not in ['everywhere', 'otherwise']:
@@ -276,16 +277,22 @@ def boundary_predicates(expr, measure, subdomain_id):
         assert measure in subdomains
         subdomain_data = subdomains[measure]
 
-        # Determine the name of the parameter function
-        name = get_global_context_value("data").object_names[id(subdomain_data)]
+        from ufl.classes import Expr
+        if isinstance(subdomain_data, Expr):
+            cond = visitor(subdomain_data)
+        else:
+            # Determine the name of the parameter function
+            cond = get_global_context_value("data").object_names[id(subdomain_data)]
+
+            # Trigger the generation of code for this thing in the parameter class
+            from ufl.checks import is_cellwise_constant
+            cellwise_constant = is_cellwise_constant(expr)
+            from dune.perftool.pdelab.parameter import intersection_parameter_function
+            intersection_parameter_function(cond, subdomain_data, cellwise_constant, t='int32')
 
-        # Trigger the generation of code for this thing in the parameter class
-        from ufl.checks import is_cellwise_constant
-        cellwise_constant = is_cellwise_constant(expr)
-        from dune.perftool.pdelab.parameter import intersection_parameter_function
-        intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int32')
+            cond = prim.Variable(cond)
 
-        predicates = predicates.union(['{} == {}'.format(name, subdomain_id)])
+        predicates = predicates.union([prim.Comparison(cond, '==', subdomain_id)])
 
     return predicates
 
@@ -331,7 +338,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
     from dune.perftool.pdelab.argument import name_accumulation_variable
     accumvar = name_accumulation_variable((ansatz_lfs.get_restriction() + test_lfs.get_restriction()))
 
-    predicates = boundary_predicates(accterm.term, measure, subdomain_id)
+    predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
 
     rank = 1 if ansatz_lfs.lfs is None else 2
 

test/poisson/poisson_neumann.ufl

 cell = triangle
+x = SpatialCoordinate(cell)
 
-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 -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 0; else return 1;", on_intersection=True)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
+sgn = conditional(x[1] > 0.5, 1., -1.)
+j = -2.*sgn*(x[1]-0.5)*g
+
+bctype = conditional(Or(x[1]<1e-8, x[1]>1.-1e-8), 0, 1)
 
 V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g, dirichlet_constraints=bctype)
 u = TrialFunction(V)
 v = TestFunction(V)
 
-x = SpatialCoordinate(cell)
-
-c = (0.5-x[0])**2 + (0.5-x[1])**2
-f = 4*(1.-c)*exp(-1.*c)
-sgn = conditional(x[1] > 0.5, 1., -1.)
-j = -2.*sgn*(x[1]-0.5)*exp(-1.*c)
-
 # Define the boundary measure that knows where we are...
 ds = ds(subdomain_data=bctype)