Fix a lot of glitches with systems

python/dune/perftool/pdelab/driver/__init__.py

@@ -164,22 +164,37 @@ def _flatten_list(l):
         yield l
 
 
+def _unroll_list_tensors(expr):
+    from ufl.classes import ListTensor
+    if isinstance(expr, ListTensor):
+        for op in expr.ufl_operands:
+            yield op
+    else:
+        yield expr
+
+
+def unroll_list_tensors(data):
+    for expr in data:
+        for e in _unroll_list_tensors(expr):
+            yield e
+
+
 def preprocess_leaf_data(element, data):
     data = get_object(data)
     from ufl import MixedElement
     if isinstance(element, MixedElement):
-        # Dirichlet is None -> no dirichlet boundaries
+        # data is None -> use 0 default
         if data is None:
-            return (0,) * element.value_size()
-        # Dirichlet for MixedElement is not iterable -> Same
-        # constraint on all the leafs.
-        elif not isinstance(data, (tuple, list)):
-            return (data,) * element.value_size()
-        # List sizes do not match -> flatten list
-        elif len(data) != element.value_size():
-            flattened = [i for i in _flatten_list(data)]
-            assert len(flattened) == element.value_size()
-            return flattened
+            data = (0,) * element.value_size()
+
+        # Flatten nested lists
+        data = tuple(i for i in _flatten_list(data))
+
+        # Expand any list tensors
+        data = tuple(i for i in unroll_list_tensors(data))
+
+        assert len(data) == element.value_size()
+        return data
     else:
         # Do not return lists for non-MixedElement
         if not isinstance(data, (tuple, list)):

python/dune/perftool/pdelab/driver/constraints.py

@@ -6,7 +6,8 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
                                          get_formdata,
                                          get_trial_element,
                                          )
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
+                                                           name_leafview,
                                                            name_trial_gfs,
                                                            type_range,
                                                            type_trial_gfs,
@@ -51,8 +52,8 @@ def name_bctype_function(element, is_dirichlet):
         for subel in element.sub_elements():
             subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
             k = k + subel.value_size()
-        name = "_".join(subgfs)
-        define_composite_bctype_function(element, is_dirichlet, name, subgfs)
+        name = "{}_bctype".format("_".join(subgfs))
+        define_composite_bctype_function(element, is_dirichlet, name, tuple(subgfs))
         return name
     else:
         assert isinstance(element, FiniteElement)

python/dune/perftool/pdelab/driver/gridfunctionspace.py

@@ -177,7 +177,7 @@ def name_test_gfs():
 
 def name_gfs(element, is_dirichlet):
     if isinstance(element, (VectorElement, TensorElement)):
-        subel = element.sub_element()[0]
+        subel = element.sub_elements()[0]
         subgfs = name_gfs(subel, is_dirichlet[:subel.value_size()])
         name = "{}_pow{}gfs".format(subgfs, element.num_sub_elements())
         define_power_gfs(element, is_dirichlet, name, subgfs)
@@ -189,7 +189,7 @@ def name_gfs(element, is_dirichlet):
             subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
             k = k + subel.value_size()
         name = "_".join(subgfs)
-        define_composite_gfs(element, is_dirichlet, name, subgfs)
+        define_composite_gfs(element, is_dirichlet, name, tuple(subgfs))
         return name
     else:
         assert isinstance(element, FiniteElement)
@@ -214,7 +214,7 @@ def type_trial_gfs():
 
 def type_gfs(element, is_dirichlet):
     if isinstance(element, (VectorElement, TensorElement)):
-        subel = element.sub_element()[0]
+        subel = element.sub_elements()[0]
         subgfs = type_gfs(subel, is_dirichlet[:subel.value_size()])
         name = "{}_POW{}GFS".format(subgfs, element.num_sub_elements())
         typedef_power_gfs(element, is_dirichlet, name, subgfs)
@@ -226,7 +226,7 @@ def type_gfs(element, is_dirichlet):
             subgfs.append(type_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
             k = k + subel.value_size()
         name = "_".join(subgfs)
-        typedef_composite_gfs(element, name, subgfs)
+        typedef_composite_gfs(element, name, tuple(subgfs))
         return name
     else:
         assert isinstance(element, FiniteElement)
@@ -249,6 +249,7 @@ def define_gfs(element, is_dirichlet, name):
 def define_power_gfs(element, is_dirichlet, name, subgfs):
     gv = name_leafview()
     fem = name_fem(element.sub_elements()[0])
+    gfstype = type_gfs(element, is_dirichlet)
     return "{} {}({}, {});".format(gfstype, name, gv, fem)
 
 

python/dune/perftool/pdelab/driver/interpolate.py

@@ -52,7 +52,7 @@ def name_boundary_function(element, func):
             childs.append(name_boundary_function(subel, func[k:k + subel.value_size()]))
             k = k + subel.value_size()
         name = "_".join(childs)
-        define_composite_boundary_function(name, childs)
+        define_compositegfs_parameterfunction(name, tuple(childs))
         return name
     else:
         assert isinstance(element, FiniteElement)

python/dune/perftool/pdelab/driver/vtk.py

@@ -140,7 +140,7 @@ def _define_gfs_name(element, is_dirichlet, offset=0):
     if 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=offset + k)
             k = k + subel.value_size()
         return []
     else:

test/stokes/stokes.mini

@@ -14,5 +14,4 @@ extension = vtu
 
 [formcompiler]
 numerical_jacobian = 0, 1 | expand num
-exact_solution_expression = g
 compare_l2errorsquared = 1e-11

test/stokes/stokes.ufl

@@ -3,10 +3,8 @@ cell = triangle
 x = SpatialCoordinate(cell)
 v_bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
 g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0))
-g_p = 8.*(1.-x[0])
-g = (g_v, g_p)
 
-P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
+P2 = VectorElement("Lagrange", cell, 2)
 P1 = FiniteElement("Lagrange", cell, 1)
 TH = P2 * P1
 
@@ -16,3 +14,6 @@ u, p = TrialFunctions(TH)
 r = (inner(grad(v), grad(u)) - div(v)*p - q*div(u))*dx
 
 forms = [r]
+is_dirichlet = v_bctype, v_bctype, 0
+dirichlet_expression = g_v, None
+exact_solution = g_v, 8.*(1.-x[0])
\ No newline at end of file