Fix error calculation for systems

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)

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()

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))

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)
+    return _define_gfs_name(element, is_dirichlet, 0)
\ No newline at end of file