firedrakeproject/firedrake — pull request #23 by wence- (2013-10-17 15:59:52+00:00)
Title: Vector function spaces

Assembling into vector function spaces is currently unsupported. There are a few reasons why:

1. The `ufl_element` on a `VectorFunctionSpace` needs to be a `VectorElement`.

2. If this is done, the extruded code breaks, I think because ufl doesn't know how to build a `MixedElement` (of which a `VectorElement` is a special case) on top of an `OuterProductCell`. This is a problem for a combination of Andrew and Doru I think.

3. Ignoring the extrusion failures, our ffc fork currently generates the wrong code for assembling into vector spaces. I think this is a known problem (Florian?)

Corollaries:

1. Vector `Expression`s don't work, that is: u.interpolate(Expression(("0", "1"))) doesn't work. Given all the other bits and pieces being in place, I think that is fixed by the following:

````patch
@@ -1008,6 +1012,9 @@ the :class:`FunctionSpace`.

to_pts = []

+ if expression.rank() != self.function_space().rank:
+ raise RuntimeError('Rank mismatch between Expression and FunctionSpace')
+
for dual in to_element.dual_basis():
if not isinstance(dual, FIAT.functional.PointEvaluation):
raise NotImplementedError("Can only interpolate onto point evaluation operators. Try projecting instead")
@@ -1019,8 +1026,9 @@ the :class:`FunctionSpace`.
# Produce C array notation of X.
X_str = "{{"+"},\n{".join([ ",".join(map(str,x)) for x in X.T])+"}}"

+ assign_expression = ";\n".join(["A[%(i)d] = %(code)s" % { 'i': i, 'code': code } for i, code in enumerate(expression.code)])
_expression_template = """
-void expression_kernel(double A[1], double **x_, int k)
+void expression_kernel(double A[%(rank)d], double **x_, int k)
{
const double X[%(ndof)d][%(xndof)d] = %(x_array)s;

@@ -1034,15 +1042,15 @@ void expression_kernel(double A[1], double **x_, int k)
};
};

- // Note this currently only does the scalar case.
- A[0] = %(expression)s;
+ %(assign_expression)s;
}
"""
kernel = op2.Kernel(_expression_template % { "x_array" : X_str,
"dim" : coords_space.dim,
"xndof" : coords_element.space_dimension(),
"ndof" : to_element.space_dimension(),
- "expression" : expression.code },
+ "assign_expression" : assign_expression,
+ "rank" : expression.rank() },
"expression_kernel")

op2.par_loop(kernel, self.cell_set,
diff --git a/python/firedrake/expression.py b/python/firedrake/expression.py
index ea8e2e6..56b558a 100644
--- a/python/firedrake/expression.py
+++ b/python/firedrake/expression.py
@@ -1,11 +1,20 @@
import ufl
-
+import numpy as np

class Expression(ufl.Coefficient):

def __init__(self, code=None, element=None, cell=None, degree=None, **kwargs):

- self.code = code
+ self.code = np.array(code)
+ shape = self.code.shape
+ self._rank = len(shape)
+ self._shape = shape
self.element = element
self.cell = cell
self.degree = degree
+
+ def rank(self):
+ return self._rank
+
+ def shape(self):
+ return self._shape
````