diff --git a/python/dune/perftool/pdelab/driver/constraints.py b/python/dune/perftool/pdelab/driver/constraints.py
index 97b8bf9dbf41157d4c03ae7dc5493ce5d96aa65b..5c9319726fefbebdaeb57fdcc28f90f05b9a806d 100644
--- a/python/dune/perftool/pdelab/driver/constraints.py
+++ b/python/dune/perftool/pdelab/driver/constraints.py
@@ -41,19 +41,19 @@ def assemble_constraints(name):
def name_bctype_function(element, is_dirichlet):
if isinstance(element, (VectorElement, TensorElement)):
- subel = element.sub_element()[0]
- subgfs = name_bctype_function(subel, is_dirichlet[:subel.value_size()])
- name = "{}_pow{}bctype".format(subgfs, element.num_sub_elements())
- define_power_bctype_function(element, name, subgfs)
+ subel = element.sub_elements()[0]
+ child = name_bctype_function(subel, is_dirichlet[:subel.value_size()])
+ name = "{}_pow{}bctype".format(child, element.num_sub_elements())
+ define_power_bctype_function(element, name, child)
return name
if isinstance(element, MixedElement):
k = 0
- subgfs = []
+ childs = []
for subel in element.sub_elements():
- subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
+ childs.append(name_bctype_function(subel, is_dirichlet[k:k + subel.value_size()]))
k = k + subel.value_size()
- name = "{}_bctype".format("_".join(subgfs))
- define_composite_bctype_function(element, is_dirichlet, name, tuple(subgfs))
+ name = "{}_bctype".format("_".join(childs))
+ define_composite_bctype_function(element, is_dirichlet, name, tuple(childs))
return name
else:
assert isinstance(element, FiniteElement)
diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
index 3dc50671c2a040d0b090532bc22e45d571d9dabf..c598f2df1c7d41f1bca46da7bc57342ead935155 100644
--- a/python/dune/perftool/pdelab/driver/error.py
+++ b/python/dune/perftool/pdelab/driver/error.py
@@ -11,6 +11,7 @@ from dune.perftool.pdelab.driver import (get_formdata,
)
from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
name_trial_gfs,
+ name_trial_subgfs,
type_range,
)
from dune.perftool.pdelab.driver.interpolate import (interpolate_vector,
@@ -20,6 +21,7 @@ from dune.perftool.pdelab.driver.solve import (define_vector,
dune_solve,
name_vector,
)
+from ufl import MixedElement
@preamble
@@ -33,9 +35,13 @@ def name_test_fail_variable():
return name
-def name_exact_solution_gridfunction():
+def name_exact_solution_gridfunction(treepath):
element = get_trial_element()
func = preprocess_leaf_data(element, "exact_solution")
+ if isinstance(element, MixedElement):
+ index = treepath_to_index(element, treepath)
+ func = (func[index],)
+ element = element.extract_component(index)[1]
return name_boundary_function(element, func)
@@ -57,54 +63,86 @@ def name_discrete_grid_function(gfs, vector_name):
@preamble
-def typedef_difference_squared_adapter(name):
- sol = name_exact_solution_gridfunction()
+def typedef_difference_squared_adapter(name, treepath):
+ sol = name_exact_solution_gridfunction(treepath)
vector = name_vector(get_formdata())
- gfs = name_trial_gfs()
+ gfs = name_trial_subgfs(treepath)
dgf = name_discrete_grid_function(gfs, vector)
return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, sol, dgf)
-def type_difference_squared_adapter():
- name = 'DifferenceSquaredAdapter'
- typedef_difference_squared_adapter(name)
+def type_difference_squared_adapter(treepath):
+ name = 'DifferenceSquaredAdapter_{}'.format("_".join(str(t) for t in treepath))
+ typedef_difference_squared_adapter(name, treepath)
return name
@preamble
-def define_difference_squared_adapter(name):
- t = type_difference_squared_adapter()
- sol = name_exact_solution_gridfunction()
+def define_difference_squared_adapter(name, treepath):
+ t = type_difference_squared_adapter(treepath)
+ sol = name_exact_solution_gridfunction(treepath)
vector = name_vector(get_formdata())
- gfs = name_trial_gfs()
+ gfs = name_trial_subgfs(treepath)
dgf = name_discrete_grid_function(gfs, vector)
return '{} {}({}, {});'.format(t, name, sol, dgf)
-def name_difference_squared_adapter():
- name = 'differencesquared_adapter'
- define_difference_squared_adapter(name)
+def name_difference_squared_adapter(treepath):
+ name = 'dsa_{}'.format("_".join(str(t) for t in treepath))
+ define_difference_squared_adapter(name, treepath)
return name
@preamble
-def accumulate_L2_squared():
- dsa = name_difference_squared_adapter()
+def _accumulate_L2_squared(treepath):
+ dsa = name_difference_squared_adapter(treepath)
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",
- "typename decltype({})::Traits::RangeType err(0.0);".format(dsa),
- "Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
- "{} += err;".format(accum_error)
+ strtp = ", ".join(str(t) for t in treepath)
+
+ return ["{",
+ " // L2 error squared of difference between numerical",
+ " // solution and the interpolation of exact solution",
+ " // for treepath ({})".format(strtp),
+ " typename decltype({})::Traits::RangeType err(0.0);".format(dsa),
+ " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+ " {} += err;".format(accum_error),
+ " std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+ "}",
]
+def get_treepath(element, index):
+ if isinstance(element, MixedElement):
+ i, rest = element.extract_subelement_component(index)
+ return (i,) + rest
+ else:
+ return ()
+
+
+def treepath_to_index(element, treepath, offset=0):
+ if len(treepath) == 0:
+ return offset
+ index = treepath[0]
+ offset = offset + sum(element.sub_elements()[i].value_size() for i in range(index))
+ subel = element.sub_elements()[treepath[0]]
+ return treepath_to_index(subel, treepath[1:], offset)
+
+
+def accumulate_L2_squared():
+ element = get_trial_element()
+ if isinstance(element, MixedElement):
+ for i in range(element.value_size()):
+ _accumulate_L2_squared(get_treepath(element, i))
+ else:
+ _accumulate_L2_squared(())
+
+
@preamble
def define_accumulated_L2_error(name):
t = type_range()
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index e264e2ab46fb2ca4d0c8a4292d2c8a8e6b7d214f..fb0533076b34351411f2486289da9c06fd958f08 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -344,3 +344,30 @@ def type_constraintsassembler(is_dirichlet):
name = "NoConstraintsAssembler"
typedef_no_constraintsassembler(name)
return name
+
+
+def name_trial_subgfs(treepath):
+ if len(treepath) == 0:
+ return name_trial_gfs()
+ else:
+ return name_subgfs(treepath)
+
+
+def name_subgfs(treepath):
+ gfs = name_trial_gfs()
+ name = "{}_{}".format(gfs, "_".join(str(t) for t in treepath))
+ define_subgfs(name, treepath)
+ return name
+
+
+@preamble
+def define_subgfs(name, treepath):
+ t = type_subgfs(treepath)
+ gfs = name_trial_gfs()
+ return "{} {}({});".format(t, name, gfs)
+
+
+def type_subgfs(tree_path):
+ include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
+ gfs = type_trial_gfs()
+ return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(str(t) for t in tree_path))
diff --git a/python/dune/perftool/pdelab/driver/vtk.py b/python/dune/perftool/pdelab/driver/vtk.py
index 15ca6e7f132871aab9284468c369003c34add62a..2d2fa20ad8a157ddbbb2c20ae54dab86923d00ef 100644
--- a/python/dune/perftool/pdelab/driver/vtk.py
+++ b/python/dune/perftool/pdelab/driver/vtk.py
@@ -135,12 +135,15 @@ def visualize_initial_condition():
@preamble
-def _define_gfs_name(element, is_dirichlet, offset=0):
- from ufl import MixedElement
- if isinstance(element, MixedElement):
+def _define_gfs_name(element, is_dirichlet, offset):
+ from ufl import MixedElement, VectorElement
+ if isinstance(element, VectorElement):
+ gfs = name_gfs(element, is_dirichlet)
+ return "{}.name(\"data_{}to{}\");".format(gfs, offset, offset + element.value_size())
+ elif isinstance(element, MixedElement):
k = 0
for subel in element.sub_elements():
- _define_gfs_name(subel, is_dirichlet[k:k + subel.value_size()], offset=offset + k)
+ _define_gfs_name(subel, is_dirichlet[k:k + subel.value_size()], offset + k)
k = k + subel.value_size()
return []
else:
@@ -151,4 +154,4 @@ def _define_gfs_name(element, is_dirichlet, offset=0):
def define_gfs_names(element):
is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- return _define_gfs_name(element, is_dirichlet)
\ No newline at end of file
+ return _define_gfs_name(element, is_dirichlet, 0)
\ No newline at end of file