From 067f6783e0f1928227d4b69e34c429f00409778b Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Thu, 28 Apr 2016 13:30:31 +0200
Subject: [PATCH] Big chunk of work piping restrictions through the
infrastructure
---
python/dune/perftool/loopy/transformer.py | 78 +++---
python/dune/perftool/pdelab/__init__.py | 10 +
python/dune/perftool/pdelab/argument.py | 64 +++--
python/dune/perftool/pdelab/basis.py | 119 +++++-----
python/dune/perftool/pdelab/driver.py | 2 +-
python/dune/perftool/pdelab/geometry.py | 101 ++++----
python/dune/perftool/pdelab/localoperator.py | 83 +++----
python/dune/perftool/pdelab/parameter.py | 10 +-
python/dune/perftool/pdelab/quadrature.py | 4 +-
python/dune/perftool/pdelab/signatures.py | 238 +++++++++++++++++++
python/dune/perftool/pymbolic/uflmapper.py | 3 +
test/laplace/CMakeLists.txt | 7 +
test/laplace/laplace_dg_numdiff.mini | 2 +-
test/laplace/laplace_dg_symdiff.mini | 6 +
14 files changed, 499 insertions(+), 228 deletions(-)
create mode 100644 python/dune/perftool/pdelab/signatures.py
create mode 100644 test/laplace/laplace_dg_symdiff.mini
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 49a3f5d1..5ec24a16 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -55,7 +55,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
if o.on_intersection:
intersection_parameter_function(name, o)
else:
- cell_parameter_function(name, o)
+ cell_parameter_function(name, o, self.restriction)
# And return a symbol
from pymbolic.primitives import Variable
@@ -101,21 +101,28 @@ def transform_accumulation_term(term, measure, subdomain_id):
rmap = {}
for ma in test_ma:
- with global_context(restriction=ma.restriction):
- # Set up the local function space structure
- traverse_lfs_tree(ma)
+ # If this is a boundary integral, all terms are implicitly restricted to the inside cell
+ if measure == 'exterior_facet':
+ ma.restriction = Restriction.NEGATIVE
+
+ # 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:
- with global_context(restriction=ma.restriction):
- # Set up the local function space structure
- traverse_lfs_tree(ma)
+ # If this is a boundary integral, all terms are implicitly restricted to the inside cell
+ if measure == 'exterior_facet':
+ ma.restriction = Restriction.NEGATIVE
+
+ # 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()
@@ -135,7 +142,6 @@ def transform_accumulation_term(term, measure, subdomain_id):
# Define a temporary variable for this expression
expr_tv_name = "expr_" + str(get_count()).zfill(4)
- expr_tv = temporary_variable(expr_tv_name)
from pymbolic.primitives import Variable
# This is a bit hacky now. To correctly determine the iname dependencies of
@@ -151,40 +157,54 @@ def transform_accumulation_term(term, measure, subdomain_id):
forced_iname_deps_is_final=True,
)
+ # Actually, if we have a cache hit, we need to change our temporary
+ from dune.perftool.generation import retrieve_cache_items
+ expr_tv_name = filter(lambda i: i.id == insn_id, retrieve_cache_items("instruction"))[0].assignee_name
+
+ # Now register the temporary variable with loopy
+ expr_tv = temporary_variable(expr_tv_name)
+
# The data that is used to collect the arguments for the accumulate function
accumargs = []
residual_shape = {}
+ arg_restr = [None] * len(test_ma)
+
+ # Determine whether the inames in the jacobian case clash. This is the case, when
+ # we have a galerkin method.
+ galerkin = False
+ if len(test_ma) == 2:
+ galerkin = (test_ma[0].argexpr.element() == test_ma[1].argexpr.element())
# Generate the code for the modified arguments:
for arg in test_ma:
from dune.perftool.pdelab.argument import pymbolic_argument
from dune.perftool.pdelab.basis import name_lfs
accumargs.append(name_lfs(arg.argexpr.element()))
- accumargs.append(lfs_iname(arg.argexpr.element(), argcount=arg.argexpr.count()))
+
+ if galerkin and arg.argexpr.count() == 1:
+ iname = lfs_iname(arg.argexpr.element(), context="arg")
+ else:
+ iname = lfs_iname(arg.argexpr.element())
+
+ accumargs.append(iname)
+
+ # Add this iname to the set of used inames for dependencies later on
+ acc_inames = acc_inames.union(frozenset({iname}))
# Determine the shape
residual_shape[arg.argexpr.number()] = name_lfs_bound(name_lfs(arg.argexpr.element()))
- from dune.perftool.pdelab.argument import name_jacobian, name_residual
- if len(test_ma) == 1:
- accumvar = name_residual()
- else:
- accumvar = name_jacobian()
+ # Determine the restriction for later
+ arg_restr[arg.argexpr.count()] = arg.restriction
+
+ from dune.perftool.pdelab.argument import name_accumulation_variable
+ accumvar = name_accumulation_variable(arg_restr)
# The residual/the jacobian should be represented through a loopy global argument
# TODO this seems still a bit hacky, esp. w.r.t. systems
shape = tuple(v for k, v in sorted(residual_shape.items(), key=lambda (k, v): k))
globalarg(accumvar, shape=shape)
- # In the jacobian case we might need to enforce the second loop!
- if len(test_ma) == 2:
- def _get_element(test_ma):
- for ma in test_ma:
- if ma.argexpr.count() == 1:
- return ma.argexpr.element()
-
- acc_inames = acc_inames.union(frozenset({lfs_iname(_get_element(test_ma), argcount=1)}))
-
from dune.perftool.pdelab.quadrature import name_factor
factor = name_factor()
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 801c50f8..0f576010 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -17,3 +17,13 @@ def name_index(index):
assert len(index) == 1
# return str(index._indices[0])
return "i_{}".format(index._indices[0].count())
+
+
+def restricted_name(name, restriction):
+ from dune.perftool import Restriction
+ if restriction == Restriction.NONE:
+ return name
+ if restriction == Restriction.POSITIVE:
+ return name + '_n'
+ if restriction == Restriction.NEGATIVE:
+ return name + '_s'
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 3ea5c1d8..f70f3fc9 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -2,7 +2,9 @@
from dune.perftool.generation import domain, iname, pymbolic_expr, symbol, globalarg
from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor
-from dune.perftool.pdelab import name_index
+from dune.perftool.pdelab import (name_index,
+ restricted_name,
+ )
from dune.perftool.pdelab.basis import (evaluate_trialfunction,
evaluate_trialfunction_gradient,
lfs_iname,
@@ -10,7 +12,7 @@ from dune.perftool.pdelab.basis import (evaluate_trialfunction,
name_basis_gradient,
name_lfs_bound,
)
-
+from dune.perftool import Restriction
from pymbolic.primitives import Subscript, Variable
@@ -18,14 +20,14 @@ from pymbolic.primitives import Subscript, Variable
def name_testfunction_gradient(ma):
assert ma.grad
- return name_basis_gradient(ma.argexpr.element())
+ return name_basis_gradient(ma.argexpr.element(), ma.restriction)
@symbol
def name_testfunction(ma):
assert not ma.grad
- return name_basis(ma.argexpr.element())
+ return name_basis(ma.argexpr.element(), ma.restriction)
@pymbolic_expr
@@ -43,17 +45,18 @@ def pymbolic_testfunction(ma):
def name_trialfunction_gradient(ma):
assert ma.grad
- evaluate_trialfunction_gradient(ma.argexpr.element(), "gradu")
- return "gradu"
+ name = restricted_name("gradu", ma.restriction)
+ evaluate_trialfunction_gradient(ma.argexpr.element(), name, ma.restriction)
+ return name
@symbol
def name_trialfunction(ma):
assert not ma.grad
- # TODO trigger the evaluation of the trial function
- evaluate_trialfunction(ma.argexpr.element(), "u")
- return "u"
+ name = restricted_name("u", ma.restriction)
+ evaluate_trialfunction(ma.argexpr.element(), name, ma.restriction)
+ return name
@pymbolic_expr
@@ -68,13 +71,13 @@ def pymbolic_trialfunction(ma):
@symbol
-def name_testfunctionspace():
- return "lfsv"
+def name_testfunctionspace(restriction):
+ return restricted_name("lfsv", restriction)
@symbol
-def name_trialfunctionspace():
- return "lfsu"
+def name_trialfunctionspace(restriction):
+ return restricted_name("lfsu", restriction)
@symbol
@@ -88,8 +91,8 @@ def type_trialfunctionspace():
@symbol
-def name_coefficientcontainer():
- return "x"
+def name_coefficientcontainer(restriction):
+ return restricted_name("x", restriction)
@symbol
@@ -98,10 +101,15 @@ def type_coefficientcontainer():
def name_argumentspace(ma):
- if ma.argexpr.number() == 0:
- return name_testfunctionspace(ma)
- if ma.argexpr.number() == 1:
- return name_trialfunctionspace(ma)
+ from ufl.classes import Argument, Coefficient
+ if isinstance(ma.argexpr, Argument):
+ if ma.argexpr.number() == 0:
+ return name_testfunctionspace(ma.restriction)
+ if ma.argexpr.number() == 1:
+ return name_trialfunctionspace(ma.restriction)
+ if isinstance(ma.argexpr, Coefficient):
+ assert ma.argexpr.count() == 0
+ return name_trialfunctionspace(ma.restriction)
# We should never encounter an argument other than 0 or 1
assert False
@@ -124,8 +132,12 @@ def pymbolic_argument(ma):
@symbol
-def name_jacobian():
- return "jac"
+def name_jacobian(restriction1, restriction2):
+ # Restrictions may only differ if NONE
+ if (restriction1 == Restriction.NONE) or (restriction2 == Restriction.NONE):
+ assert restriction1 == restriction2
+
+ return restricted_name(restricted_name("jac", restriction1), restriction2)
@symbol
@@ -134,8 +146,8 @@ def type_jacobian():
@symbol
-def name_residual():
- return "r"
+def name_residual(restriction):
+ return restricted_name("r", restriction)
@symbol
@@ -143,13 +155,13 @@ def type_residual():
return "R"
-def name_accumulation_variable():
+def name_accumulation_variable(restrictions=(Restriction.NONE,)):
from dune.perftool.generation import get_global_context_value
ft = get_global_context_value("form_type")
if ft == 'residual':
- return name_residual()
+ return name_residual(*restrictions)
if ft == 'jacobian':
- return name_jacobian()
+ return name_jacobian(*restrictions)
assert False
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index b8d5a07c..38cb06a6 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -22,6 +22,7 @@ from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
lop_template_test_gfs,
)
from dune.perftool.pdelab.driver import FEM_name_mangling
+from dune.perftool.pdelab import restricted_name
@symbol
@@ -114,27 +115,10 @@ def traverse_lfs_tree(arg):
# First we need to determine the basename as given in the signature of
# this kernel method!
- lfs_basename = None
- gfs_basename = None
- from ufl.classes import Argument, Coefficient
- if isinstance(arg.argexpr, Argument):
- if arg.argexpr.count() == 0:
- lfs_basename = 'lfsv'
- gfs_basename = 'GFSV'
- if arg.argexpr.count() == 1:
- lfs_basename = 'lfsu'
- gfs_basename = 'GFSU'
- # TODO add restrictions here.
-
- if isinstance(arg.argexpr, Coefficient):
- # We should only ever call this for a trialfunction, which in our case
- # is the coefficient of reserved index 0.
- assert arg.argexpr.count() == 0
-
- lfs_basename = 'lfsu'
- gfs_basename = 'GFSU'
-
- assert lfs_basename and gfs_basename
+ from dune.perftool.pdelab.argument import name_argumentspace
+ lfs_basename = name_argumentspace(arg)
+ from dune.perftool.pdelab.localoperator import lop_template_gfs
+ gfs_basename = lop_template_gfs(arg)[1]
# Now start recursively extracting local function spaces and fill the cache with
# all those values. That way we can later get a correct local function space with
@@ -144,13 +128,10 @@ def traverse_lfs_tree(arg):
@iname
-def _lfs_iname(element, argcount, context):
+def _lfs_iname(element, context):
name = name_lfs(element)
bound = name_lfs_bound(name)
- if argcount != 0:
- name = 'lfsu'
-
if context:
context = '_' + context
@@ -160,7 +141,7 @@ def _lfs_iname(element, argcount, context):
return name
-def lfs_iname(element, argcount=0, context=''):
+def lfs_iname(element, context=''):
""" Get the iname to iterate over the local function space given by element
Arguments:
@@ -174,7 +155,7 @@ def lfs_iname(element, argcount=0, context=''):
a given purpose, see the 'Loops and dependencies' of the loopy docs:
https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies
"""
- return _lfs_iname(element, argcount, context)
+ return _lfs_iname(element, context)
@preamble
@@ -188,12 +169,12 @@ def declare_cache_temporary(element, name, which):
@cached
-def evaluate_basis(element, name):
+def evaluate_basis(element, name, restriction):
temporary_variable(name, shape=(name_lfs_bound(element),), decl_method=None)
declare_cache_temporary(element, name, which='Function')
cache = name_localbasis_cache(element)
lfs = name_lfs(element)
- qp = name_quadrature_position_in_cell()
+ qp = name_quadrature_position_in_cell(restriction)
instruction(inames=(quadrature_iname(),
),
code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
@@ -207,18 +188,19 @@ def evaluate_basis(element, name):
@symbol
-def name_basis(element):
- evaluate_basis(element, "phi")
- return "phi"
+def name_basis(element, restriction):
+ name = restricted_name("phi", restriction)
+ evaluate_basis(element, name, restriction)
+ return name
@cached
-def evaluate_reference_gradient(element, name):
+def evaluate_reference_gradient(element, name, restriction):
temporary_variable(name, shape=(name_lfs_bound(element), name_dimension()), decl_method=None)
declare_cache_temporary(element, name, which='Jacobian')
cache = name_localbasis_cache(element)
lfs = name_lfs(element)
- qp = name_quadrature_position_in_cell()
+ qp = name_quadrature_position_in_cell(restriction)
instruction(inames=(quadrature_iname(),
),
code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
@@ -232,18 +214,18 @@ def evaluate_reference_gradient(element, name):
@symbol
-def name_reference_gradient(element):
- evaluate_reference_gradient(element, "js")
- return "js"
+def name_reference_gradient(element, restriction):
+ name = restricted_name("js", restriction)
+ evaluate_reference_gradient(element, name, restriction)
+ return name
@cached
-def evaluate_basis_gradient(element, name):
- # TODO this is of course not yet correct
+def evaluate_basis_gradient(element, name, restriction):
temporary_variable(name, shape=(name_lfs_bound(element), name_dimension()), shape_impl=('vec', 'fv'))
- jac = name_jacobian_inverse_transposed()
+ jac = name_jacobian_inverse_transposed(restriction)
index = lfs_iname(element, context='transformgrads')
- reference_gradients = name_reference_gradient(element)
+ reference_gradients = name_reference_gradient(element, restriction)
instruction(inames=(index,
quadrature_iname(),
),
@@ -259,28 +241,42 @@ def evaluate_basis_gradient(element, name):
@symbol
-def name_basis_gradient(element):
- evaluate_basis_gradient(element, "gradphi")
- return "gradphi"
+def name_basis_gradient(element, restriction):
+ name = restricted_name("gradphi", restriction)
+ evaluate_basis_gradient(element, name, restriction)
+ return name
+
+
+def reset_trialfunction(name):
+ return instruction(inames=(quadrature_iname(),
+ ),
+ code="{} = 0.0;".format(name),
+ assignees=frozenset({name}),
+ )
@cached
-def evaluate_trialfunction(element, name):
+def evaluate_trialfunction(element, name, restriction):
temporary_variable(name, shape=())
+ reset = reset_trialfunction(name)
lfs = name_lfs(element)
index = lfs_iname(element)
- basis = name_basis(element)
+ basis = name_basis(element, restriction)
+ from dune.perftool.pdelab.argument import name_coefficientcontainer
+ coeffs = name_coefficientcontainer(restriction)
instruction(inames=(quadrature_iname(),
index,
),
- code='{} += x({}, {})*{}[{}];'.format(name,
- lfs,
- index,
- basis,
- index
- ),
+ code='{} += {}({}, {})*{}[{}];'.format(name,
+ coeffs,
+ lfs,
+ index,
+ basis,
+ index
+ ),
assignees=frozenset({name}),
read_variables=frozenset({basis}),
+ depends_on=frozenset({reset}),
)
@@ -293,22 +289,25 @@ def reset_trialfunction_gradient(name):
@cached
-def evaluate_trialfunction_gradient(element, name):
+def evaluate_trialfunction_gradient(element, name, restriction):
# TODO this is of course not yet correct
temporary_variable(name, shape=(name_dimension(),), shape_impl=('fv',))
reset = reset_trialfunction_gradient(name)
lfs = name_lfs(element)
index = lfs_iname(element, context='trialgrad')
- basis = name_basis_gradient(element)
+ basis = name_basis_gradient(element, restriction)
+ from dune.perftool.pdelab.argument import name_coefficientcontainer
+ coeffs = name_coefficientcontainer(restriction)
instruction(inames=(quadrature_iname(),
index,
),
- code='{}.axpy(x({}, {}), {}[{}]);'.format(name,
- lfs,
- index,
- basis,
- index,
- ),
+ code='{}.axpy({}({}, {}), {}[{}]);'.format(name,
+ coeffs,
+ lfs,
+ index,
+ basis,
+ index,
+ ),
assignees=frozenset({name}),
read_variables=frozenset({basis}),
forced_iname_deps=frozenset({quadrature_iname(), index}),
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 56dfb860..d8d69a50 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -263,7 +263,7 @@ def typedef_fem(expr, name):
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)
+ 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")
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 4bbfa6d5..1cdd8d62 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -1,9 +1,11 @@
from dune.perftool import Restriction
+from dune.perftool.pdelab import restricted_name
from dune.perftool.generation import (preamble,
symbol,
temporary_variable,
)
from dune.perftool.pdelab.quadrature import (name_quadrature_position,
+ name_quadrature_position_in_cell,
quadrature_preamble,
)
from ufl.algorithms import MultiFunction
@@ -94,43 +96,40 @@ def define_restricted_cell(name, restriction):
@symbol
-def _name_cell(restriction):
+def name_cell(restriction):
if restriction == Restriction.NONE:
eg = name_element_geometry_wrapper()
return "{}.entity()".format(eg)
-
- which = "inside" if restriction == Restriction.NEGATIVE else "outside"
- name = "cell_{}".format(which)
- define_restricted_cell(name, restriction)
- return name
-
-
-def name_cell():
- from dune.perftool.generation import get_global_context_value
- it = get_global_context_value("integral_type")
-
- if it == 'cell':
- r = Restriction.NONE
- if it == 'exterior_facet':
- r = Restriction.NEGATIVE
- if it == 'interior_facet':
- raise NotImplementedError
-
- return _name_cell(r)
+ else:
+ which = "inside" if restriction == Restriction.NEGATIVE else "outside"
+ name = "{}_cell".format(which)
+ define_restricted_cell(name, restriction)
+ return name
@preamble
-def define_cell_geometry(name):
- eg = name_element_geometry_wrapper()
+def define_cell_geometry(name, restriction):
+ cell = name_cell(restriction)
return "auto {} = {}.geometry();".format(name,
- eg
+ cell
)
@symbol
-def name_cell_geometry():
- define_cell_geometry("cell_geo")
- return "cell_geo"
+def name_cell_geometry(restriction):
+ name = restricted_name("cell_geo", restriction)
+ define_cell_geometry(name, restriction)
+ return name
+
+
+@symbol
+def type_cell_geometry(restriction):
+ if restriction == Restriction.NONE:
+ eg = type_element_geometry_wrapper()
+ return "{}::Geometry".format(eg)
+ else:
+ ig = type_intersection_geometry_wrapper()
+ return "{}::Entity::Geometry".format(ig)
@preamble
@@ -159,7 +158,7 @@ def name_geometry():
it = get_global_context_value("integral_type")
if it == 'cell':
- return name_cell_geometry()
+ return name_cell_geometry(Restriction.NONE)
if it == 'exterior_facet' or it == 'interior_facet':
return name_intersection_geometry()
assert False
@@ -167,7 +166,6 @@ def name_geometry():
@preamble
def define_in_cell_geometry(restriction, name):
- cell = _name_cell(restriction)
ig = name_intersection_geometry_wrapper()
which = "In" if restriction == Restriction.NEGATIVE else "Out"
return "auto {} = {}.geometryIn{}side();".format(name,
@@ -203,16 +201,10 @@ def name_in_cell_coordinates(local, basename, restriction):
return name
-def to_cell_coordinates(local, basename):
- from dune.perftool.generation import get_global_context_value
- it = get_global_context_value("integral_type")
-
- if it == 'cell':
+def to_cell_coordinates(local, basename, restriction):
+ if restriction == Restriction.NONE:
return local
- if it == 'exterior_facet':
- return name_in_cell_coordinates(local, basename, Restriction.NEGATIVE)
- if it == 'interior_facet':
- restriction = get_global_context_value("restriction")
+ else:
return name_in_cell_coordinates(local, basename, restriction)
@@ -267,24 +259,26 @@ def name_unit_inner_normal():
@symbol
-def type_jacobian_inverse_transposed():
- geo = type_element_geometry_wrapper()
- return "typename {}::Geometry::JacobianInverseTransposed".format(geo)
+def type_jacobian_inverse_transposed(restriction):
+ geo = type_cell_geometry(restriction)
+ return "typename {}::JacobianInverseTransposed".format(geo)
-@preamble
-def define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
- t = type_jacobian_inverse_transposed()
- return "{} {};".format(t,
- name,
- )
+def define_jacobian_inverse_transposed_temporary(restriction):
+ @preamble
+ def _define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
+ t = type_jacobian_inverse_transposed(restriction)
+ return "{} {};".format(t,
+ name,
+ )
+ return _define_jacobian_inverse_transposed_temporary
-def define_jacobian_inverse_transposed(name):
+def define_jacobian_inverse_transposed(name, restriction):
dim = name_dimension()
- temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary, shape=(dim, dim))
- geo = name_cell_geometry()
- pos = name_quadrature_position()
+ temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
+ geo = name_cell_geometry(restriction)
+ pos = name_quadrature_position_in_cell(restriction)
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
pos,
@@ -294,6 +288,7 @@ def define_jacobian_inverse_transposed(name):
@symbol
-def name_jacobian_inverse_transposed():
- define_jacobian_inverse_transposed("jit")
- return "jit"
+def name_jacobian_inverse_transposed(restriction):
+ name = restricted_name("jit", restriction)
+ define_jacobian_inverse_transposed(name, restriction)
+ return name
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 90555510..97e1f7e1 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -15,7 +15,7 @@ from dune.perftool.cgen.clazz import (AccessModifier,
BaseClass,
ClassMember,
)
-
+from dune.perftool import Restriction
from pytools import memoize
@@ -29,6 +29,14 @@ def lop_template_test_gfs():
return "GFSV"
+def lop_template_gfs(ma):
+ if ma.argexpr.count() == 0:
+ return lop_template_test_gfs()
+ if ma.argexpr.count() == 1:
+ return lop_template_ansatz_gfs()
+ assert False
+
+
@symbol
def name_initree_constructor_param():
return "iniParams"
@@ -116,56 +124,29 @@ def assembly_routine_signature():
integral_type = get_global_context_value("integral_type")
form_type = get_global_context_value("form_type")
- aj = "alpha" if form_type == 'residual' else "jacobian"
- which = ufl_measure_to_pdelab_measure(integral_type).lower()
-
- if integral_type == 'interior_facet':
- raise NotImplementedError
-
- from dune.perftool.pdelab.geometry import (name_geometry_wrapper,
- type_geometry_wrapper,
- )
- geot = type_geometry_wrapper()
- geo = name_geometry_wrapper()
-
- from dune.perftool.pdelab.argument import (name_accumulation_variable,
- type_accumulation_variable,
- name_coefficientcontainer,
- type_coefficientcontainer,
- name_testfunctionspace,
- type_testfunctionspace,
- name_trialfunctionspace,
- type_trialfunctionspace,
- )
- lfsut = type_trialfunctionspace()
- lfsu = name_trialfunctionspace()
- lfsvt = type_testfunctionspace()
- lfsv = name_testfunctionspace()
- cct = type_coefficientcontainer()
- cc = name_coefficientcontainer()
- avt = type_accumulation_variable()
- av = name_accumulation_variable()
-
- return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
- lfsut,
- cct,
- lfsvt,
- avt,
- ),
- 'void {}_{}(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(aj,
- which,
- geot,
- geo,
- lfsut,
- lfsu,
- cct,
- cc,
- lfsvt,
- lfsv,
- avt,
- av,
- )
- ]
+ if form_type == 'residual':
+ if integral_type == 'cell':
+ from dune.perftool.pdelab.signatures import alpha_volume_signature
+ return alpha_volume_signature()
+ if integral_type == 'exterior_facet':
+ from dune.perftool.pdelab.signatures import alpha_boundary_signature
+ return alpha_boundary_signature()
+ if integral_type == 'interior_facet':
+ from dune.perftool.pdelab.signatures import alpha_skeleton_signature
+ return alpha_skeleton_signature()
+
+ if form_type == 'jacobian':
+ if integral_type == 'cell':
+ from dune.perftool.pdelab.signatures import jacobian_volume_signature
+ return jacobian_volume_signature()
+ if integral_type == 'exterior_facet':
+ from dune.perftool.pdelab.signatures import jacobian_boundary_signature
+ return jacobian_boundary_signature()
+ if integral_type == 'interior_facet':
+ from dune.perftool.pdelab.signatures import jacobian_skeleton_signature
+ return jacobian_skeleton_signature()
+
+ assert False
def generate_kernel(integral):
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index dc717063..abffb46d 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -62,10 +62,10 @@ def define_parameter_function_class_member(name, expr, t, cell):
return result
-def evaluate_cell_parameter_function(name):
+def evaluate_cell_parameter_function(name, restriction):
param = name_paramclass()
- entity = name_cell()
- pos = name_quadrature_position_in_cell()
+ entity = name_cell(restriction)
+ pos = name_quadrature_position_in_cell(restriction)
return quadrature_preamble('{} = {}.{}({}, {});'.format(name,
name_paramclass(),
name,
@@ -105,10 +105,10 @@ def define_parameter_temporary(t):
return _define
-def cell_parameter_function(name, expr, t='double'):
+def cell_parameter_function(name, expr, restriction, t='double'):
temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
define_parameter_function_class_member(name, expr, t, True)
- evaluate_cell_parameter_function(name)
+ evaluate_cell_parameter_function(name, restriction)
def intersection_parameter_function(name, expr, t='double'):
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 2e717f76..a874c923 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -30,9 +30,9 @@ def name_quadrature_position():
@symbol
-def name_quadrature_position_in_cell():
+def name_quadrature_position_in_cell(restriction):
from dune.perftool.pdelab.geometry import to_cell_coordinates
- return to_cell_coordinates(name_quadrature_position(), name_quadrature_point())
+ return to_cell_coordinates(name_quadrature_position(), name_quadrature_point(), restriction)
@symbol
diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py
new file mode 100644
index 00000000..b0e02593
--- /dev/null
+++ b/python/dune/perftool/pdelab/signatures.py
@@ -0,0 +1,238 @@
+""" Signatures for PDELab local opreator assembly functions """
+
+from dune.perftool import Restriction
+from dune.perftool.generation import symbol
+from dune.perftool.pdelab.geometry import (name_geometry_wrapper,
+ type_geometry_wrapper,
+ )
+from dune.perftool.pdelab.argument import (name_accumulation_variable,
+ type_accumulation_variable,
+ name_coefficientcontainer,
+ type_coefficientcontainer,
+ name_testfunctionspace,
+ type_testfunctionspace,
+ name_trialfunctionspace,
+ type_trialfunctionspace,
+ )
+
+
+@symbol
+def alpha_volume_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NONE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NONE)
+ cct = type_coefficientcontainer()
+ cc = name_coefficientcontainer(Restriction.NONE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NONE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void alpha_volume(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ cc,
+ lfsvt,
+ lfsv,
+ avt,
+ av,
+ )
+ ]
+
+
+@symbol
+def alpha_boundary_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NEGATIVE)
+ cct = type_coefficientcontainer()
+ cc = name_coefficientcontainer(Restriction.NEGATIVE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NEGATIVE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void alpha_boundary(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ cc,
+ lfsvt,
+ lfsv,
+ avt,
+ av,
+ )
+ ]
+
+
+@symbol
+def alpha_skeleton_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
+ lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
+ cct = type_coefficientcontainer()
+ cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
+ cc_n = name_coefficientcontainer(Restriction.POSITIVE)
+ avt = type_accumulation_variable()
+ av_s = name_accumulation_variable((Restriction.NEGATIVE,))
+ av_n = name_accumulation_variable((Restriction.POSITIVE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void alpha_skeleton(const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu_s,
+ cct,
+ cc_s,
+ lfsvt,
+ lfsv_s,
+ lfsut,
+ lfsu_n,
+ cct,
+ cc_n,
+ lfsvt,
+ lfsv_n,
+ avt,
+ av_s,
+ avt,
+ av_n,
+ )
+ ]
+
+
+@symbol
+def jacobian_volume_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NONE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NONE)
+ cct = type_coefficientcontainer()
+ cc = name_coefficientcontainer(Restriction.NONE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NONE, Restriction.NONE))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_volume(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ cc,
+ lfsvt,
+ lfsv,
+ avt,
+ av,
+ )
+ ]
+
+
+@symbol
+def jacobian_boundary_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NEGATIVE)
+ cct = type_coefficientcontainer()
+ cc = name_coefficientcontainer(Restriction.NEGATIVE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NEGATIVE, Restriction.NEGATIVE))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_boundary(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ cc,
+ lfsvt,
+ lfsv,
+ avt,
+ av,
+ )
+ ]
+
+
+@symbol
+def jacobian_skeleton_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
+ lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
+ cct = type_coefficientcontainer()
+ cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
+ cc_n = name_coefficientcontainer(Restriction.POSITIVE)
+ avt = type_accumulation_variable()
+ av_ss = name_accumulation_variable((Restriction.NEGATIVE, Restriction.NEGATIVE))
+ av_sn = name_accumulation_variable((Restriction.NEGATIVE, Restriction.POSITIVE))
+ av_ns = name_accumulation_variable((Restriction.POSITIVE, Restriction.NEGATIVE))
+ av_nn = name_accumulation_variable((Restriction.POSITIVE, Restriction.POSITIVE))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_skeleton(const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}, {}& {}, {}& {}, {}& {}) const'.format(geot,
+ geo,
+ lfsut,
+ lfsu_s,
+ cct,
+ cc_s,
+ lfsvt,
+ lfsv_s,
+ lfsut,
+ lfsu_n,
+ cct,
+ cc_n,
+ lfsvt,
+ lfsv_n,
+ avt,
+ av_ss,
+ avt,
+ av_sn,
+ avt,
+ av_ns,
+ avt,
+ av_nn,
+ )
+ ]
diff --git a/python/dune/perftool/pymbolic/uflmapper.py b/python/dune/perftool/pymbolic/uflmapper.py
index 854c4762..19d6f250 100644
--- a/python/dune/perftool/pymbolic/uflmapper.py
+++ b/python/dune/perftool/pymbolic/uflmapper.py
@@ -58,3 +58,6 @@ class UFL2PymbolicMapper(MultiFunction):
def sum(self, o):
return Sum(tuple(self.call(op) for op in get_operands(o)))
+
+ def zero(self, o):
+ return 0
diff --git a/test/laplace/CMakeLists.txt b/test/laplace/CMakeLists.txt
index e268854d..ae7f7d55 100644
--- a/test/laplace/CMakeLists.txt
+++ b/test/laplace/CMakeLists.txt
@@ -22,3 +22,10 @@ add_generated_executable(UFLFILE laplace_dg.ufl
dune_add_system_test(TARGET laplace_dg_numdiff
INIFILE laplace_dg_numdiff.mini)
+
+add_generated_executable(UFLFILE laplace_dg.ufl
+ TARGET laplace_dg_symdiff
+ )
+
+dune_add_system_test(TARGET laplace_dg_symdiff
+ INIFILE laplace_dg_symdiff.mini)
diff --git a/test/laplace/laplace_dg_numdiff.mini b/test/laplace/laplace_dg_numdiff.mini
index bf8ddc7a..ad527212 100644
--- a/test/laplace/laplace_dg_numdiff.mini
+++ b/test/laplace/laplace_dg_numdiff.mini
@@ -1,4 +1,4 @@
-__name = laplacg_dg_numdiff
+__name = laplace_dg_numdiff
lowerleft = 0.0 0.0
upperright = 1.0 1.0
diff --git a/test/laplace/laplace_dg_symdiff.mini b/test/laplace/laplace_dg_symdiff.mini
new file mode 100644
index 00000000..edefa9db
--- /dev/null
+++ b/test/laplace/laplace_dg_symdiff.mini
@@ -0,0 +1,6 @@
+__name = laplace_dg_symdiff
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 4 4
+elementType = simplical
--
GitLab