diff --git a/.gitmodules b/.gitmodules
index 0de93ddeb2191afbd96edcb149447785602fbb95..982bfd26f3903235440402a9fedd5777b2a21785 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -11,3 +11,6 @@
path = python/ufl
url = https://bitbucket.org/fenics-project/ufl.git
ignore = untracked
+[submodule "python/pymbolic"]
+ path = python/pymbolic
+ url = https://github.com/inducer/pymbolic.git
diff --git a/python/dune/perftool/__init__.py b/python/dune/perftool/__init__.py
index eb279cf11e48b26b210e9e010b4371635c939f61..02ef874ff476ba7fc9bca5ecf261f0fcf7073edb 100644
--- a/python/dune/perftool/__init__.py
+++ b/python/dune/perftool/__init__.py
@@ -1,4 +1,4 @@
class Restriction:
NONE = 0
- INSIDE = 1
- OUTSIDE = 2
+ NEGATIVE = 1
+ POSITIVE = 2
diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 83265ae709953d9731f41ee069acb1768f2b1b04..231792238a7e4e7c60494a241775e178881c3ec2 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -49,6 +49,8 @@ def read_ufl(uflfile):
form = transform_form(form, reindexing)
# form = transform_form(form, split_arguments)
+ formdata.preprocessed_form = form
+
return formdata, data.object_names
diff --git a/python/dune/perftool/interactive.py b/python/dune/perftool/interactive.py
index a492bc12d0bd564c65efdb837a8fd5fd3ee09647..b84edff8377f366bb051a4ef91db5942c8bf53be 100644
--- a/python/dune/perftool/interactive.py
+++ b/python/dune/perftool/interactive.py
@@ -73,9 +73,11 @@ def show_code(which, kernel):
clear()
print("Showing the generated dune-pdelab code for {}:\n".format(kernel_name(which)))
- from dune.perftool.pdelab.localoperator import measure_specific_details, AssemblyMethod
- signature = measure_specific_details(which[0])["{}_signature".format(which[1])]
- print("".join(AssemblyMethod(signature, kernel).generate()))
+ from dune.perftool.generation import global_context
+ with global_context(integral_type=which[0], form_type=which[1]):
+ from dune.perftool.pdelab.localoperator import assembly_routine_signature, AssemblyMethod
+ signature = assembly_routine_signature()
+ print("".join(AssemblyMethod(signature, kernel).generate()))
print("Press Return to return to the previous menu")
input()
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 21c2265d0ebac22dbe7ca454ae20fdd330021203..271abc12123f1e15eb79fc120c545016f9a5e800 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -8,8 +8,9 @@ from __future__ import absolute_import
from dune.perftool import Restriction
from dune.perftool.ufl.modified_terminals import ModifiedTerminalTracker
from dune.perftool.pymbolic.uflmapper import UFL2PymbolicMapper
-
+from dune.perftool.pdelab.geometry import GeometryMapper
from dune.perftool.generation import (domain,
+ global_context,
globalarg,
iname,
instruction,
@@ -31,7 +32,7 @@ def index_sum_iname(i):
return name_index(i)
-class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper):
+class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapper):
def __init__(self):
super(UFL2LoopyVisitor, self).__init__()
@@ -73,12 +74,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper):
# Now continue processing the expression
return self.call(o.ufl_operands[0])
- # TODO use multiple inheritance and have a geometry transformer in the pdelab subpackage
- def facet_area(self, o):
- from pymbolic.primitives import Variable
- from dune.perftool.pdelab.geometry import name_facetarea
- return Variable(name_facetarea())
-
class _Counter:
counter = 0
@@ -106,19 +101,21 @@ def transform_accumulation_term(term, measure, subdomain_id):
rmap = {}
for ma in test_ma:
- # Set up the local function space structure
- traverse_lfs_tree(ma)
+ with global_context(restriction=ma.restriction):
+ # Set up the local function space structure
+ traverse_lfs_tree(ma)
- # Get the expression for the modified argument representing the test function
- from dune.perftool.pdelab.argument import pymbolic_testfunction
- rmap[ma.expr] = pymbolic_testfunction(ma)
+ # Get the expression for the modified argument representing the test function
+ from dune.perftool.pdelab.argument import pymbolic_testfunction
+ rmap[ma.expr] = pymbolic_testfunction(ma)
for ma in trial_ma:
- # Set up the local function space structure
- traverse_lfs_tree(ma)
+ with global_context(restriction=ma.restriction):
+ # Set up the local function space structure
+ traverse_lfs_tree(ma)
- # Get the expression for the modified argument representing the trial function
- from dune.perftool.pdelab.argument import pymbolic_trialfunction
- rmap[ma.expr] = pymbolic_trialfunction(ma)
+ # Get the expression for the modified argument representing the trial function
+ from dune.perftool.pdelab.argument import pymbolic_trialfunction
+ rmap[ma.expr] = pymbolic_trialfunction(ma)
# Get the transformer!
ufl2l_mf = UFL2LoopyVisitor()
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 4d5762dd4687a2e7fa8cd7b8e4d2a4865484c3e6..801c50f8ed32307a8018b62bc2d5c91e3e170643 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -10,7 +10,10 @@ from dune.perftool.generation import (preamble,
def name_index(index):
from ufl.classes import MultiIndex, Index
if isinstance(index, Index):
- return str(index)
+ # This failed for index > 9 because ufl placed curly brackets around
+ # return str(index)
+ return "i_{}".format(index.count())
if isinstance(index, MultiIndex):
assert len(index) == 1
- return str(index._indices[0])
+ # return str(index._indices[0])
+ return "i_{}".format(index._indices[0].count())
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 2015ac7322debd8ddd2d477cac755ef8cc9e3140..b8d5a07c6b4e2a24e4ffd7f34d903f4afc26cf43 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -269,7 +269,7 @@ def evaluate_trialfunction(element, name):
temporary_variable(name, shape=())
lfs = name_lfs(element)
index = lfs_iname(element)
- basis = name_basis()
+ basis = name_basis(element)
instruction(inames=(quadrature_iname(),
index,
),
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 556fbfb818ed12ab3efb2f226afb2b448c9971e3..56dfb86092a6c2cae9b151cbb68e94e227528c50 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -63,6 +63,10 @@ def isQk(fem):
return isLagrange(fem) and isQuadrilateral(fem)
+def isDG(fem):
+ return fem._short_name is 'DG'
+
+
def FEM_name_mangling(fem):
from ufl import MixedElement, VectorElement, FiniteElement
if isinstance(fem, MixedElement):
@@ -79,6 +83,9 @@ def FEM_name_mangling(fem):
return "P" + str(fem._degree)
if isQk(fem):
return "Q" + str(fem._degree)
+ if isDG(fem):
+ return "DG" + str(fem._degree)
+
raise NotImplementedError("FEM NAME MANGLING")
@@ -242,12 +249,23 @@ def typedef_fem(expr, name):
gv = type_leafview()
df = type_domainfield()
r = type_range()
+ dim = name_dimension()
if isPk(expr):
include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
return "typedef Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}> {};".format(gv, df, r, expr._degree, name)
- if isQk(generator._kwargs['expr']):
+ if isQk(expr):
include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
return "typedef Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}> {};".format(gv, df, r, expr._degree, name)
+ if isDG(expr):
+ if isQuadrilateral(expr):
+ include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
+ # TODO allow switching the basis here!
+ return "typedef Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}> {}".format(df, r, expr._degree, dim, name)
+ if isSimplical(expr):
+ include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
+ return "typedef Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::cube> {}".format(df, r, expr._degree, dim, name)
+ raise NotImplementedError("Geometry type not known in code generation")
+
raise NotImplementedError("FEM not implemented in dune-perftool")
@@ -261,8 +279,11 @@ def type_fem(expr):
@preamble
def define_fem(expr, name):
femtype = type_fem(expr)
- gv = name_leafview()
- return "{} {}({});".format(femtype, name, gv)
+ if isDG(expr):
+ return "{} {};".format(femtype, name)
+ else:
+ gv = name_leafview()
+ return "{} {}({});".format(femtype, name, gv)
@symbol
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 4e0f5ab6e75b70f644bbfa44e6c8993ce63700ab..4bbfa6d54d5d25a28e139da8ffe46383f91b389a 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -6,6 +6,37 @@ from dune.perftool.generation import (preamble,
from dune.perftool.pdelab.quadrature import (name_quadrature_position,
quadrature_preamble,
)
+from ufl.algorithms import MultiFunction
+
+
+class GeometryMapper(MultiFunction):
+ """
+ A collection of visitors for geometry related UFL nodes
+ NB: This is kind of 'abstract' as it needs to be combined
+ with a ModifiedTerminalTracker through multi inheritance.
+ """
+ def __init__(self):
+ super(GeometryMapper, self).__init__()
+
+ def facet_normal(self, o):
+ # The normal must be restricted to be well-defined
+ assert self.restriction is not Restriction.NONE
+
+ from pymbolic.primitives import Variable
+ if self.restriction == Restriction.POSITIVE:
+ return Variable(name_unit_outer_normal())
+ if self.restriction == Restriction.NEGATIVE:
+ # It is highly unnatural to have this generator function,
+ # but I do run into subtle trouble with return -1*outer
+ # as the indexing into the normal happens only later.
+ # Not investing more time into this cornercase right now.
+ return Variable(name_unit_inner_normal())
+
+ # TODO This one was just copied over so, it might need some tweaking
+ def facet_area(self, o):
+ from pymbolic.primitives import Variable
+ from dune.perftool.pdelab.geometry import name_facetarea
+ return Variable(name_facetarea())
@symbol
@@ -55,7 +86,7 @@ def type_geometry_wrapper():
@preamble
def define_restricted_cell(name, restriction):
ig = name_intersection_geometry_wrapper()
- which = "inside" if restriction == Restriction.INSIDE else "outside"
+ which = "inside" if restriction == Restriction.NEGATIVE else "outside"
return "const auto& {} = {}.{}();".format(name,
ig,
which,
@@ -68,7 +99,7 @@ def _name_cell(restriction):
eg = name_element_geometry_wrapper()
return "{}.entity()".format(eg)
- which = "inside" if restriction == Restriction.INSIDE else "outside"
+ which = "inside" if restriction == Restriction.NEGATIVE else "outside"
name = "cell_{}".format(which)
define_restricted_cell(name, restriction)
return name
@@ -81,7 +112,7 @@ def name_cell():
if it == 'cell':
r = Restriction.NONE
if it == 'exterior_facet':
- r = Restriction.INSIDE
+ r = Restriction.NEGATIVE
if it == 'interior_facet':
raise NotImplementedError
@@ -138,7 +169,7 @@ def name_geometry():
def define_in_cell_geometry(restriction, name):
cell = _name_cell(restriction)
ig = name_intersection_geometry_wrapper()
- which = "In" if restriction == Restriction.INSIDE else "Out"
+ which = "In" if restriction == Restriction.NEGATIVE else "Out"
return "auto {} = {}.geometryIn{}side();".format(name,
ig,
which
@@ -149,7 +180,7 @@ def define_in_cell_geometry(restriction, name):
def name_in_cell_geometry(restriction):
assert restriction is not Restriction.NONE
- name = "geo_in_{}side".format("in" if restriction is Restriction.INSIDE else "out")
+ name = "geo_in_{}side".format("in" if restriction is Restriction.NEGATIVE else "out")
define_in_cell_geometry(restriction, name)
return name
@@ -168,7 +199,7 @@ def apply_in_cell_transformation(name, local, restriction):
def name_in_cell_coordinates(local, basename, restriction):
name = "{}_in_inside".format(basename)
temporary_variable(name, shape=(name_dimension(),), shape_impl=("fv",))
- apply_in_cell_transformation(name, local, restriction=Restriction.INSIDE)
+ apply_in_cell_transformation(name, local, restriction=Restriction.NEGATIVE)
return name
@@ -179,9 +210,10 @@ def to_cell_coordinates(local, basename):
if it == 'cell':
return local
if it == 'exterior_facet':
- return name_in_cell_coordinates(local, basename, Restriction.INSIDE)
+ return name_in_cell_coordinates(local, basename, Restriction.NEGATIVE)
if it == 'interior_facet':
- raise NotImplementedError
+ restriction = get_global_context_value("restriction")
+ return name_in_cell_coordinates(local, basename, restriction)
@preamble
@@ -196,6 +228,44 @@ def name_dimension():
return "dim"
+def evaluate_unit_outer_normal(name):
+ ig = name_intersection_geometry_wrapper()
+ qp = name_quadrature_position()
+ return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ assignees=frozenset({name}),
+ )
+
+
+@preamble
+def declare_normal(name, shape, shape_impl):
+ ig = name_intersection_geometry_wrapper()
+ return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)
+
+
+@symbol
+def name_unit_outer_normal():
+ name = "outer_normal"
+ temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
+ evaluate_unit_outer_normal(name)
+ return "outer_normal"
+
+
+def evaluate_unit_inner_normal(name):
+ outer = name_unit_outer_normal()
+ return quadrature_preamble("auto {} = -1. * {};".format(name, outer),
+ assignees=frozenset({name}),
+ read_variables=frozenset({outer}),
+ )
+
+
+@symbol
+def name_unit_inner_normal():
+ name = "inner_normal"
+ temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
+ evaluate_unit_inner_normal(name)
+ return "inner_normal"
+
+
@symbol
def type_jacobian_inverse_transposed():
geo = type_element_geometry_wrapper()
diff --git a/python/dune/perftool/pymbolic/uflmapper.py b/python/dune/perftool/pymbolic/uflmapper.py
index 18b03cc1c82b95ed4d53f07d9b1a3420eb2e3d47..854c476235165ff0f617a96bf9bb3f1b7da692af 100644
--- a/python/dune/perftool/pymbolic/uflmapper.py
+++ b/python/dune/perftool/pymbolic/uflmapper.py
@@ -24,12 +24,12 @@ class UFL2PymbolicMapper(MultiFunction):
return Product(tuple(self.call(op) for op in get_operands(o)))
def multi_index(self, o):
- return tuple(self.call(op) for op in o.ufl_operands)
+ from dune.perftool.pdelab import name_index
+ return tuple(Variable(name_index(op)) for op in o.indices())
def index(self, o):
# One might as well take the uflname as string here, but I apply this function
from dune.perftool.pdelab import name_index
-
return Variable(name_index(o))
def fixed_index(self, o):
diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py
index 6beed959b0c0cd372114be603a8eeaae2190f122..b992dc05f2fe76fb68c76b336bd46cf580c96277 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -18,14 +18,22 @@ class ModifiedTerminalTracker(MultiFunction):
def positive_restricted(self, o):
assert self.restriction == Restriction.NONE
self.restriction = Restriction.POSITIVE
- ret = self.call(o.ufl_operands[0])
+
+ from dune.perftool.generation import global_context
+ with global_context(restriction=Restriction.POSITIVE):
+ ret = self.call(o.ufl_operands[0])
+
self.restriction = Restriction.NONE
return ret
def negative_restricted(self, o):
assert self.restriction == Restriction.NONE
self.restriction = Restriction.NEGATIVE
- ret = self.call(o.ufl_operands[0])
+
+ from dune.perftool.generation import global_context
+ with global_context(restriction=Restriction.NEGATIVE):
+ ret = self.call(o.ufl_operands[0])
+
self.restriction = Restriction.NONE
return ret
diff --git a/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py b/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py
index 2b8505c56d77e790ad875bc454bd995c1423f2b8..c783d88229a995eafe741d9c1d0443c822dd08dc 100644
--- a/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py
+++ b/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py
@@ -19,7 +19,7 @@ class _ReplacementDict(dict):
def __init__(self, good=[], bad=[]):
dict.__init__(self)
for a in bad:
- self[a] = Zero()
+ self[a] = Zero(shape=a.ufl_shape, free_indices=a.ufl_free_indices, index_dimensions=a.ufl_index_dimensions)
for a in good:
self[a] = a
diff --git a/python/pymbolic b/python/pymbolic
new file mode 160000
index 0000000000000000000000000000000000000000..80ef7a7745829f53ae949f68ad458b88c06d66a0
--- /dev/null
+++ b/python/pymbolic
@@ -0,0 +1 @@
+Subproject commit 80ef7a7745829f53ae949f68ad458b88c06d66a0
diff --git a/test/laplace/laplace_dg.ufl b/test/laplace/laplace_dg.ufl
index f1bffa46e0cf2945d06f29cafb83b89b5cb06297..3c89fbb7648410f3117a781bc41ca12333deff57 100644
--- a/test/laplace/laplace_dg.ufl
+++ b/test/laplace/laplace_dg.ufl
@@ -4,14 +4,17 @@ V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell)
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
- - inner(n, avg(grad(u)))*jump(v)*(dS+ds) \
- + gamma*jump(u)*jump(v)*(dS+ds) \
- - theta*jump(u)*inner(grad(v), n)*(dS+ds)
+ - inner(n, avg(grad(u)))*jump(v)*dS \
+ + gamma*jump(u)*jump(v)*dS \
+ - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+ - inner(n, grad(u))*v*ds \
+ + gamma*u*v*ds \
+ - theta*u*inner(grad(v), n)*ds
forms = [r]