diff --git a/cmake/modules/StandardMain.cmake b/cmake/modules/StandardMain.cmake
index 8c8dd056d9916538d2a216ffff6e1a8db88803d8..028c2efc0b208705845326fc9a621697c72ec408 100644
--- a/cmake/modules/StandardMain.cmake
+++ b/cmake/modules/StandardMain.cmake
@@ -18,9 +18,10 @@ int main(int argc, char** argv)
std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
<<" processes!"<<std::endl;
- driver(argc, argv);
-
- return 0;
+ if (driver(argc, argv))
+ return 1;
+ else
+ return 0;
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 78e872f6b3410c0323f0b91edeca496c7c7132cc..0129f4d70888be3db7f4ef1dd418956252754fba 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -236,9 +236,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
# Check if this is a parameter function
else:
- # We expect all coefficients to be of type Expression or VectorExpression!
- from dune.perftool.ufl.execution import Expression, VectorExpression
- assert isinstance(o, (Expression, VectorExpression))
+ # We expect all coefficients to be of type Expression!
+ from dune.perftool.ufl.execution import Expression
+ assert isinstance(o, Expression)
# Determine the name of the parameter function
from dune.perftool.generation import get_global_context_value
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 57ae37c02406953ba799be1322cbc24f1ca42ac5..1fca4361b09ee9758a05924a01a5a6af46ba04c3 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -423,9 +423,9 @@ 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)
+ 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)
+ return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
@symbol
@@ -741,9 +741,9 @@ def define_boundary_lambda(boundary, name):
if boundary is None:
return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
if boundary.is_global:
- return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr)
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr[0])
else:
- return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr)
+ return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr[0])
@symbol
@@ -777,41 +777,61 @@ def define_compositegfs_parameterfunction(name, *children):
', '.join(children))
+def expression_splitter(expr, length):
+ if expr is None:
+ return (None,) * length
+ else:
+ from dune.perftool.ufl.execution import split_expression
+ return split_expression(expr)
+
+
@symbol
-def name_boundary_function(expr):
+def _name_boundary_function(element, boundary, name=None):
from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(expr, (VectorElement, TensorElement)):
- element, (_, boundary) = get_constraints_metadata(expr)
+ if isinstance(element, (VectorElement, TensorElement)):
from dune.perftool.generation import get_global_context_value
- name = get_global_context_value("data").object_names.get(id(boundary), "zero")
- define_boundary_function(boundary, name)
- return name
- if isinstance(expr, MixedElement):
- children = tuple(name_boundary_function(el) for el in expr.sub_elements())
+ basename = get_global_context_value("data").object_names.get(id(boundary), "veczero")
+ children = tuple(_name_boundary_function(el, exp, basename + '_' + str(i)) for el, exp, i in zip(element.sub_elements(), expression_splitter(boundary, element.num_sub_elements()), range(element.num_sub_elements())))
+ define_compositegfs_parameterfunction(basename, *children)
+ return basename
+ if isinstance(element, MixedElement):
+ children = tuple(name_boundary_function(el) for el in element.sub_elements())
name = '_'.join(children)
define_compositegfs_parameterfunction(name, *children)
return name
- # This is a scalar leaf of the ansatz function tree
- element, (_, boundary) = get_constraints_metadata(expr)
-
from dune.perftool.generation import get_global_context_value
- name = get_global_context_value("data").object_names.get(id(boundary), "zero")
+ if name is None:
+ name = get_global_context_value("data").object_names.get(id(boundary), "zero")
define_boundary_function(boundary, name)
return name
+def name_boundary_function(expr):
+ # This is a scalar leaf of the ansatz function tree
+ element, (_, boundary) = get_constraints_metadata(expr)
+
+ return _name_boundary_function(element, boundary)
+
+
@symbol
-def name_solution_function(expr):
- # Get expression representing solution
+def name_solution_function(tree_path=()):
+ from dune.perftool.generation import get_global_context_value
+ expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
+
+ from dune.perftool.ufl.execution import split_expression
+ for i in tree_path:
+ expr = split_expression(expr)[int(i)]
+
from dune.perftool.generation import get_global_context_value
- solution_expression_name = get_option("exact_solution_expression")
- solution_expression = get_global_context_value("data").object_by_name[solution_expression_name]
+ try:
+ name = get_global_context_value("data").object_names[id(expr)]
+ except KeyError:
+ from dune.perftool.generation.cache import _GlobalCounter
+ name = 'generic_{}'.format(_GlobalCounter().get())
- # Create function just like it is done for boundary_function
- name = "solution_function"
- define_boundary_function(solution_expression, name)
+ define_boundary_function(expr, name)
return name
@@ -833,7 +853,7 @@ def interpolate_solution_expression(name):
formdata = _driver_data['formdata']
define_vector(name, formdata)
element = _driver_data['form'].coefficients()[0].ufl_element()
- sol = name_solution_function(element)
+ sol = name_solution_function()
gfs = name_gfs(element)
return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
gfs,
@@ -1156,6 +1176,7 @@ def name_vtkfile():
def compare_dofs():
v = name_vector(_driver_data['formdata'])
solution = name_vector_from_solution_expression()
+ fail = name_test_fail_variable()
return ["",
"// Maximum norm of difference between dof vectors of the",
"// numerical solution and the interpolation of the exact solution",
@@ -1166,44 +1187,161 @@ def compare_dofs():
" maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
"std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
"if (maxerror>{})".format(get_option("compare_dofs")),
- " throw;"]
+ " {} = true;".format(fail)]
+
+
+@symbol
+def type_discrete_grid_function(gfs):
+ return "DGF_{}".format(gfs.upper())
@preamble
def define_discrete_grid_function(gfs, vector_name, dgf_name):
- return ["typedef Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})> DGF;".format(gfs, vector_name),
- "DGF {}({},{});".format(dgf_name, gfs, vector_name)]
+ dgf_type = type_discrete_grid_function(gfs)
+ return ["typedef Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})> {};".format(gfs, vector_name, dgf_type),
+ "{} {}({},{});".format(dgf_type, dgf_name, gfs, vector_name)]
@symbol
def name_discrete_grid_function(gfs, vector_name):
- dgf_name = vector_name + "_dgf"
+ dgf_name = gfs + "_dgf"
define_discrete_grid_function(gfs, vector_name, dgf_name)
return dgf_name
+@symbol
+def type_subgfs(element, tree_path):
+ include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
+ gfs = type_gfs(element)
+ return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(tree_path))
+
+
@preamble
-def compare_L2_squared():
+def define_subgfs(name, element, tree_path):
+ t = type_subgfs(element, tree_path)
+ gfs = name_gfs(element)
+
+ return '{} {}({});'.format(t, name, gfs)
+
+
+@symbol
+def name_subgfs(element, tree_path):
+ gfs = name_gfs(element)
+ if len(tree_path) == 0:
+ return gfs
+ else:
+ name = "subgfs_{}_{}".format(gfs, '_'.join(tree_path))
+ define_subgfs(name, element, tree_path)
+ return name
+
+
+@preamble
+def typedef_difference_squared_adapter(name, tree_path):
element = _driver_data['form'].coefficients()[0].ufl_element()
formdata = _driver_data['formdata']
- v = name_vector(formdata)
- gfs = name_gfs(element)
- vdgf = name_discrete_grid_function(gfs, v)
- solution_function = name_solution_function(element)
+
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return 'typedef DifferenceSquaredAdapter<decltype({}), decltype({})> {};'.format(solution_function, dgf, name)
+
+
+@symbol
+def type_difference_squared_adapter(tree_path):
+ t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
+ typedef_difference_squared_adapter(t, tree_path)
+ return t
+
+
+@preamble
+def define_difference_squared_adapter(name, tree_path):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ formdata = _driver_data['formdata']
+
+ t = type_difference_squared_adapter(tree_path)
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return '{} {}({}, {});'.format(t, name, solution_function, dgf)
+
+
+def name_difference_squared_adapter(tree_path):
+ name = 'differencesquared_{}'.format('_'.join(tree_path))
+ define_difference_squared_adapter(name, tree_path)
+ return name
+
+
+@preamble
+def _accumulate_L2_squared(tree_path):
+ dsa = name_difference_squared_adapter(tree_path)
+ t_dsa = type_difference_squared_adapter(tree_path)
+ accum_error = name_accumulated_L2_error()
include_file("dune/pdelab/gridfunctionspace/gridfunctionadapter.hh", filetag="driver")
include_file("dune/pdelab/common/functionutilities.hh", filetag="driver")
- return ["",
- "// L2 error squared of difference between numerical",
- "// solution and the interpolation of exact solution",
- "typedef DifferenceSquaredAdapter<decltype({}),decltype({})> DifferenceSquared;".format(solution_function, vdgf),
- "DifferenceSquared differencesquared({},{});".format(solution_function, vdgf),
- "typename DifferenceSquared::Traits::RangeType l2errorsquared(0.0);",
- "Dune::PDELab::integrateGridFunction(differencesquared,l2errorsquared,10);",
- "std::cout << \"l2errorsquared: \" << l2errorsquared << std::endl;",
- "if (l2errorsquared>{})".format(get_option("compare_l2errorsquared")),
- " throw;"]
+ return ["{",
+ " // L2 error squared of difference between numerical",
+ " // solution and the interpolation of exact solution",
+ " // only subgfs with treepath: {}".format(', '.join(tree_path)),
+ " typename {}::Traits::RangeType err(0.0);".format(t_dsa),
+ " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+ " {} += err;".format(accum_error),
+ " std::cout << \"Accumlated L2 Error for tree path {}: \" << err << std::endl;".format(tree_path),
+ "}"]
+
+
+def accumulate_L2_squared(tree_path=()):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ from ufl.functionview import select_subelement
+ from ufl.classes import MultiIndex, FixedIndex
+ element = select_subelement(element, MultiIndex(tuple(FixedIndex(int(i)) for i in tree_path)))
+
+ from ufl import MixedElement
+ if isinstance(element, MixedElement):
+ for i, subel in enumerate(element.sub_elements()):
+ accumulate_L2_squared(tree_path + (str(i),))
+ else:
+ _accumulate_L2_squared(tree_path)
+
+
+@preamble
+def define_accumulated_L2_error(name):
+ return "double {}(0.0);".format(name)
+
+
+@symbol
+def name_accumulated_L2_error():
+ name = 'l2error'
+ define_accumulated_L2_error(name)
+ return name
+
+
+@preamble
+def compare_L2_squared():
+ accumulate_L2_squared()
+
+ accum_error = name_accumulated_L2_error()
+ fail = name_test_fail_variable()
+ return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
+ "if ({}>{})".format(accum_error, get_option("compare_l2errorsquared")),
+ " {} = true;".format(fail)]
+
+
+@preamble
+def define_test_fail_variable(name):
+ return 'bool {}(false);'.format(name)
+
+
+@symbol
+def name_test_fail_variable():
+ name = "testfail"
+ define_test_fail_variable(name)
+ return name
@preamble
@@ -1289,10 +1427,6 @@ def vtkoutput():
print_matrix()
if get_option("exact_solution_expression"):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
- raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
-
if get_option("compare_dofs"):
compare_dofs()
if get_option("compare_l2errorsquared"):
@@ -1424,16 +1558,18 @@ def solve_instationary():
print_residual()
print_matrix()
if get_option("exact_solution_expression"):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
- raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
-
if get_option("compare_dofs"):
compare_dofs()
if get_option("compare_l2errorsquared"):
compare_L2_squared()
+@preamble
+def return_statement():
+ fail = name_test_fail_variable()
+ return "return {};".format(fail)
+
+
def generate_driver(formdatas, data):
# The driver module uses a global dictionary for storing necessary data
set_driver_data(formdatas, data)
@@ -1445,9 +1581,11 @@ def generate_driver(formdatas, data):
else:
solve_instationary()
+ return_statement()
+
from dune.perftool.generation import retrieve_cache_items
from cgen import FunctionDeclaration, FunctionBody, Block, Value
- driver_signature = FunctionDeclaration(Value('void', 'driver'), [Value('int', 'argc'), Value('char**', 'argv')])
+ driver_signature = FunctionDeclaration(Value('bool', 'driver'), [Value('int', 'argc'), Value('char**', 'argv')])
driver_body = Block(contents=[i for i in retrieve_cache_items("preamble", make_generable=True)])
driver = FunctionBody(driver_signature, driver_body)
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index 79fede9c1eb72e20613b2db92cd6df8b1463050c..624f82dbaa97c4e39bc9ad0223e38613874b4ded 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -69,8 +69,22 @@ def define_set_time_method():
return result
+def component_to_tree_path(element, component):
+ subel = element.extract_subelement_component(component)
+
+ def _flatten(x):
+ if isinstance(x, tuple):
+ return '_'.join(_flatten(i) for i in x if i != ())
+ else:
+ return str(x)
+
+ return _flatten(subel)
+
+
@class_member("parameterclass", access=AccessModifier.PUBLIC)
-def define_parameter_function_class_member(name, expr, t, cell):
+def define_parameter_function_class_member(name, expr, baset, shape, cell):
+ t = construct_nested_fieldvector(baset, shape)
+
geot = "E" if cell else "I"
geo = geot.lower()
result = ["template<typename {}, typename X>".format(geot),
@@ -78,12 +92,32 @@ def define_parameter_function_class_member(name, expr, t, cell):
"{",
]
- if expr.is_global:
- result.append(" auto x = {}.geometry().global(local);".format(geo))
+ # In the case of a non-scalar parameter function, recurse into leafs
+ if expr.element.value_shape():
+ # Check that this is a VectorElement, as I have no idea how a parameter function
+ # over a non-vector mixed element should be well-defined in PDELab.
+ from ufl import VectorElement
+ assert isinstance(expr.element, VectorElement)
+
+ result.append(" {} result(0.0);".format(t))
+
+ from dune.perftool.ufl.execution import split_expression
+ for i, subexpr in enumerate(split_expression(expr)):
+ child_name = "{}_{}".format(name, component_to_tree_path(expr.element, i))
+ result.append(" result[{}] = {}({}, local);".format(i, child_name, geo))
+ define_parameter_function_class_member(child_name, subexpr, baset, shape[1:], cell)
+
+ result.append(" return result;")
+
else:
- result.append(" auto x = local;")
+ # Evaluate a scalar parameter function
+ if expr.is_global:
+ result.append(" auto x = {}.geometry().global(local);".format(geo))
+ else:
+ result.append(" auto x = local;")
+
+ result.append(" " + expr.c_expr[0])
- result.append(" " + expr.c_expr)
result.append("}")
return result
@@ -175,8 +209,7 @@ def construct_nested_fieldvector(t, shape):
def cell_parameter_function(name, expr, restriction, cellwise_constant, t='double'):
shape = expr.ufl_element().value_shape()
shape_impl = ('fv',) * len(shape)
- t = construct_nested_fieldvector(t, shape)
- define_parameter_function_class_member(name, expr, t, True)
+ define_parameter_function_class_member(name, expr, t, shape, True)
if cellwise_constant:
from dune.perftool.generation.loopy import default_declaration
default_declaration(name, shape, shape_impl)
@@ -190,8 +223,7 @@ def cell_parameter_function(name, expr, restriction, cellwise_constant, t='doubl
def intersection_parameter_function(name, expr, cellwise_constant, t='double'):
shape = expr.ufl_element().value_shape()
shape_impl = ('fv',) * len(shape)
- t = construct_nested_fieldvector(t, shape)
- define_parameter_function_class_member(name, expr, t, False)
+ define_parameter_function_class_member(name, expr, t, shape, False)
if cellwise_constant:
from dune.perftool.generation.loopy import default_declaration
default_declaration(name, shape, shape_impl)
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index fb0aa6b5d008dd028f3d7903c3b57a43ed728842..32717a9ed2b703eaf2fe44e9081e61fdc98f82a2 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -38,7 +38,10 @@ class Coefficient(ufl.Coefficient):
ufl.Coefficient.__init__(self, element, count)
-split = ufl.split_functions.split2
+def split(obj):
+ if isinstance(obj, Expression):
+ return split_expression(obj)
+ return ufl.split_functions.split2(obj)
def Coefficients(element):
@@ -54,42 +57,66 @@ def TrialFunctions(element):
class Expression(Coefficient):
- def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle", cellwise_constant=False):
- assert isinstance(expr, str)
- self.c_expr = expr
- self.is_global = is_global
- self.on_intersection = on_intersection
-
- # Avoid ufl complaining about not matching dimension/cells
- if cellwise_constant:
- _dg = FiniteElement("DG", cell_type, 0)
- else:
- _dg = FiniteElement("DG", cell_type, 1)
+ def __init__(self, cppcode=None, element=None, cell="triangle", degree=None, is_global=True, on_intersection=False, cellwise_constant=False):
+ assert cppcode
- # Initialize a coefficient with a dummy finite element map.
- Coefficient.__init__(self, _dg)
+ if isinstance(cppcode, str):
+ cppcode = (cppcode,)
- # TODO the subdomain_data code needs a uflid, not idea how to get it here
- # The standard way through class decorator fails here...
- def ufl_id(self):
- return 0
+ def wrap_return(e):
+ if "return" not in e:
+ return "return {};".format(e)
+ else:
+ return e
+ cppcode = tuple(wrap_return(e) for e in cppcode)
-class VectorExpression(Coefficient):
- def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle", cellwise_constant=False):
- assert isinstance(expr, str)
- self.c_expr = expr
+ if cellwise_constant:
+ if element:
+ assert element.degree() == 0
+ elif degree is not None:
+ assert degree == 0
+ else:
+ element = FiniteElement("DG", cell, 0)
+
+ if degree is None:
+ degree = 1
+
+ if element is None:
+ if len(cppcode) == 1:
+ element = FiniteElement("DG", cell, degree)
+ else:
+ element = VectorElement("DG", cell, degree, len(cppcode))
+
+ def element_length(elem):
+ if isinstance(elem, ufl.FiniteElement):
+ return 1
+ else:
+ return elem.value_shape()[0]
+
+ assert element_length(element) == len(cppcode)
+
+ self.c_expr = cppcode
self.is_global = is_global
self.on_intersection = on_intersection
-
- # Avoid ufl complaining about not matching dimension/cells
- if cellwise_constant:
- _dgvec = VectorElement("DG", cell_type, 0)
- else:
- _dgvec = VectorElement("DG", cell_type, 1)
+ self.element = element
# Initialize a coefficient with a dummy finite element map.
- Coefficient.__init__(self, _dgvec)
+ Coefficient.__init__(self, element)
+
+ def __mul__(self, other):
+ # Allow the combination of Expressions
+ if isinstance(other, Expression):
+ def get_sub(elem):
+ if isinstance(elem, MixedElement) and not isinstance(elem, (VectorElement, TensorElement)):
+ return elem.sub_elements()
+ else:
+ return [elem]
+
+ newelem = MixedElement(*(get_sub(self.element) + get_sub(other.element)))
+ return Expression(cppcode=self.c_expr + other.c_expr, element=newelem)
+ else:
+ return Coefficient.__mul__(self, other)
# TODO the subdomain_data code needs a uflid, not idea how to get it here
# The standard way through class decorator fails here...
@@ -97,6 +124,19 @@ class VectorExpression(Coefficient):
return 0
+def split_expression(expr):
+ assert isinstance(expr, Expression)
+
+ def element_slice(expression, sub):
+ offset = sum(subel.value_size() for subel in expression.element.sub_elements()[:sub])
+ return expression.c_expr[offset:offset + expr.element.sub_elements()[sub].value_size()]
+
+ return tuple(Expression(cppcode=element_slice(expr, i),
+ element=expr.element.sub_elements()[i])
+ for i in range(expr.element.num_sub_elements())
+ )
+
+
class FiniteElement(ufl.FiniteElement):
def __init__(self, *args, **kwargs):
if ('dirichlet_constraints' in kwargs) or ('dirichlet_expression' in kwargs):
diff --git a/test/heatequation/heatequation.mini b/test/heatequation/heatequation.mini
index 6087d4160a3084387be6fa5af703745893d6dfde..b5897d254180c71726bf0b181e6ae72e6afa8714 100644
--- a/test/heatequation/heatequation.mini
+++ b/test/heatequation/heatequation.mini
@@ -15,3 +15,6 @@ extension = vtu
explicit_time_stepping = 0, 1 | expand scheme
exact_solution_expression = g
compare_l2errorsquared = 1e-7
+
+# Disable explicit tests for now
+{__exec_suffix} == explicit | exclude
diff --git a/test/heatequation/heatequation_dg.mini b/test/heatequation/heatequation_dg.mini
index e2a455372a077c34f9a46f239da702363a8cdde1..47559b75a86880bddb1134be30a31dfa4494b8e0 100644
--- a/test/heatequation/heatequation_dg.mini
+++ b/test/heatequation/heatequation_dg.mini
@@ -15,3 +15,6 @@ extension = vtu
explicit_time_stepping = 0, 1 | expand scheme
exact_solution_expression = g
compare_l2errorsquared = 1e-5
+
+# Disable explicit tests for now
+{__exec_suffix} == explicit | exclude
diff --git a/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl b/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
index 864810e8a0dd87330992901f0badaf53c9eedbe3..6bda534f16236d66b351c0ac1c0671a41f03f9de 100644
--- a/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
@@ -1,9 +1,9 @@
-cell_type = "interval"
+cell = "interval"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", cell=cell)
-V = FiniteElement("CG", cell_type, 1, dirichlet_expression=g)
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
v = TestFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
index f7775baf11b5225addc511a2016ddaac8ad0d9fe..907f19841d51695f3944617ff162b33fce7b5053 100644
--- a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
@@ -1,14 +1,14 @@
-cell_type = "interval"
+cell = "interval"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", on_intersection=True, cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
-V = FiniteElement("DG", cell_type, 1)
+V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell_type)('+')
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
index 7d43c7eb3ee86e180ad51443d2d4bb370e79e4c2..5a2c570583a2fe200e91a1dbf20b0dde699ec156 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
@@ -1,9 +1,9 @@
-cell_type = "quadrilateral"
+cell = "quadrilateral"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", cell=cell)
-V = FiniteElement("CG", cell_type, 1, dirichlet_expression=g)
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
v = TestFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
index 24dbccae56b8b93ae305b561e1ac5025d3a53bd7..1ca65e1eaf16d4ff71444bd71626c277b82328c9 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
@@ -1,14 +1,14 @@
-cell_type = "quadrilateral"
+cell = "quadrilateral"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", on_intersection=True, cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
-V = FiniteElement("DG", cell_type, 1)
+V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell_type)('+')
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
index 5c0091d74b1e9349fd6631b18b2ffdfbd42e990d..4cad67e31f354a7eb0c4c6d25b95713f31067612 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
@@ -1,9 +1,9 @@
-cell_type = "hexahedron"
+cell = "hexahedron"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", cell=cell)
-V = FiniteElement("CG", cell_type, 1, dirichlet_expression=g)
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
v = TestFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
index 1456671d43ce0de2d3bfbf8d1ca13fda9f8e0733..d8f43c0941267480a07eba857b93d635571af43d 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
@@ -1,7 +1,7 @@
-cell_type = "tetrahedron"
+cell = "tetrahedron"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", cell=cell)
V = FiniteElement("CG", "tetrahedron", 1, dirichlet_expression=g)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
index 3a410e3011596ab6372028e1e705d26cd5747216..bf5cc370a5651c1a408d2f1b6687d4c3255c0172 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
@@ -1,14 +1,14 @@
-cell_type = "hexahedron"
+cell = "hexahedron"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", on_intersection=True, cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
-V = FiniteElement("DG", cell_type, 1)
+V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell_type)('+')
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
index a2d8135baa6d8f6556dcff607989af38a9ffa826..4c73c7debf2ac3a5f81941d2e6d12cb5ee4ff4fa 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
@@ -1,14 +1,14 @@
-cell_type = "tetrahedron"
+cell = "tetrahedron"
-f = Expression("return -2.0*x.size();", cell_type=cell_type)
-g = Expression("return x.two_norm2();", on_intersection=True, cell_type=cell_type)
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
-V = FiniteElement("DG", cell_type, 1)
+V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell_type)('+')
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index 8d6fa96d619fc308cba1a75296dc91c84fe43413..0d3751b8483822901ea723d51bd14e71a68657fb 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -5,17 +5,14 @@ dune_symlink_to_source_files(FILES hagenpoiseuille_ref.vtu
dune_add_formcompiler_system_test(UFLFILE stokes.ufl
BASENAME stokes
INIFILE stokes.mini
- SCRIPT dune_vtkcompare.py
)
dune_add_formcompiler_system_test(UFLFILE stokes_dg.ufl
BASENAME stokes_dg
INIFILE stokes_dg.mini
- SCRIPT dune_vtkcompare.py
)
dune_add_formcompiler_system_test(UFLFILE stokes_stress.ufl
BASENAME stokes_stress
INIFILE stokes_stress.mini
- SCRIPT dune_vtkcompare.py
)
diff --git a/test/stokes/stokes.mini b/test/stokes/stokes.mini
index 23a7932929a7783a92ee032727b9749022bdfe47..6bb834ff2fa259e99a2fa835cd4564765d63cd9b 100644
--- a/test/stokes/stokes.mini
+++ b/test/stokes/stokes.mini
@@ -14,3 +14,5 @@ extension = vtu
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes.ufl b/test/stokes/stokes.ufl
index 9252834820f67c72a6b28f08b6b1bd1b164ad04c..da632caf9ecbf0f93eb5aa1c70cd5df17aba01b6 100644
--- a/test/stokes/stokes.ufl
+++ b/test/stokes/stokes.ufl
@@ -1,8 +1,11 @@
v_bctype = Expression("if (x[0] < 1. - 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;")
+g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
+g_p = Expression("8*(1.-x[0])")
+
+g = g_v * g_p
cell = triangle
-P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=v_dirichlet)
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
P1 = FiniteElement("Lagrange", cell, 1)
TH = P2 * P1
diff --git a/test/stokes/stokes_dg.mini b/test/stokes/stokes_dg.mini
index 6448f4f5b36dfad891d76e069c3e5573546e6f85..20d92508ac79daf523858c1c99b4476c2cc0b65e 100644
--- a/test/stokes/stokes_dg.mini
+++ b/test/stokes/stokes_dg.mini
@@ -16,3 +16,5 @@ zeroThreshold.data_1 = 1e-6
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-5
diff --git a/test/stokes/stokes_dg.ufl b/test/stokes/stokes_dg.ufl
index 7b15e24ce48c5aac90f274581d3e13f44ae80fb5..b63893f0fe98622e1b24019e930026c615074010 100644
--- a/test/stokes/stokes_dg.ufl
+++ b/test/stokes/stokes_dg.ufl
@@ -1,4 +1,6 @@
-g = VectorExpression("Dune::FieldVector<double, 2> y(0.0); y[0]=4*x[1]*(1.-x[1]); return y;", on_intersection=True)
+g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"), on_intersection=True)
+g_p = Expression("8*(1.-x[0])")
+g = g_v * g_p
bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
cell = triangle
@@ -19,15 +21,15 @@ r = inner(grad(u), grad(v))*dx \
+ inner(avg(grad(u))*n, jump(v))*dS \
- eps * inner(avg(grad(v))*n, jump(u))*dS \
- inner(grad(u)*n, v)*ds \
- + eps * inner(grad(v)*n, u-g)*ds \
+ + eps * inner(grad(v)*n, u-g_v)*ds \
+ sigma * inner(jump(u), jump(v))*dS \
- + sigma * inner(u-g, v)*ds \
+ + sigma * inner(u-g_v, v)*ds \
- p*div(v)*dx \
- avg(p)*inner(jump(v), n)*dS \
+ p*inner(v, n)*ds \
- q*div(u)*dx \
- avg(q)*inner(jump(u), n)*dS \
+ q*inner(u, n)*ds \
- - q*inner(g, n)*ds
+ - q*inner(g_v, n)*ds
forms = [r]
diff --git a/test/stokes/stokes_stress.mini b/test/stokes/stokes_stress.mini
index 642ef6d2ea8768f1fad13f29c6773e2ee173b830..4481e91d21e683dc8452c84c35ca03eca7a2182c 100644
--- a/test/stokes/stokes_stress.mini
+++ b/test/stokes/stokes_stress.mini
@@ -14,3 +14,5 @@ extension = vtu
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes_stress.ufl b/test/stokes/stokes_stress.ufl
index 9eb4f3b750e26bb930acda17387da3644773d044..d5c41e301ce7919fbd94fb86f9c6cab250ae9083 100644
--- a/test/stokes/stokes_stress.ufl
+++ b/test/stokes/stokes_stress.ufl
@@ -1,16 +1,19 @@
v_bctype = Expression("if (x[0] < 1. - 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;")
+g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
+g_p = Expression("8*(1.-x[0])")
+g_S = Expression(("0.0", "0.0", "-8*x[1] + 4.", "0.0"))
+g = g_v * g_p * g_S
cell = triangle
-P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=v_dirichlet)
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
P1 = FiniteElement("Lagrange", cell, 1)
-P2_stress = TensorElement("Lagrange", cell, 2)
+P2_stress = TensorElement("DG", cell, 1)
TH = MixedElement(P2, P1, P2_stress)
v, q, T = TestFunctions(TH)
u, p, S = TrialFunctions(TH)
-r = (inner(grad(v), S) + inner(S - grad(u), T)- div(v)*p - q*div(u))*dx
+r = (inner(grad(v), S) + inner(grad(u) - S, T) - div(v)*p - q*div(u))*dx
forms = [r]