Composing ufl block matrices from tensors

1

I would like to compose blocks of arbitrary ufl tensors as in the attached example, such that the resulting ufl object is a 2D matrix rather than a higher-order tensor. Is there a simpler way to get the desired result or do I have do everything using indices as in the MWE?

from dolfin import *
mesh = UnitSquareMesh(2,2)
I = Identity(2)
test = as_matrix([[I, -I],
                  [-I, I]])
# This results in a tensor of shape (2,2,2,2)
print(project(test, TensorFunctionSpace(mesh, 'CG', 1, (2,2,2,2)))(0.5,.5).reshape(2,2,2,2))

# want a matrix of shape (4,4)
desired = as_matrix([[I[0,0], I[0,1], -I[0,0], -I[0,1]],
                     [I[1,0], I[1,1], -I[1,0], -I[1,1]],
                     [-I[0,0], -I[0,1], I[0,0], I[0,1]],
                     [-I[1,0], -I[1,1], I[0,0], I[1,1]]])

print(project(desired, TensorFunctionSpace(mesh, 'CG', 1, (4,4)))(0.5,.5).reshape(4,4))
2

I don’t know of a built-in way to do this, but you could hack something together yourself, like the following:

from dolfin import *
mesh = UnitSquareMesh(2,2)

# Test with non-symmetric matrix:
#I = Identity(2)
I = as_matrix([[1,2],
               [3,4]])

from ufl import shape
def as_block_matrix(mlist):
    rows = []
    # Block rows:
    for i in range(0,len(mlist)):
        block_row_height = shape(mlist[i][0])[0]
        # Matrix rows in block row i:
        for j in range(0,block_row_height):
            row = []
            # Block columns:
            for k in range(0,len(mlist[i])):
                block_column_width = shape(mlist[i][k])[1]
                # Matrix columns in block column k:
                for l in range(0,block_column_width):
                    row += [mlist[i][k][j,l],]
            rows += [row,]
    return as_matrix(rows)

# Check using MWE:
test = as_block_matrix([[I, -I],
                        [-I, I]])
# (Contained an index typo in original MWE)
desired = as_matrix([[I[0,0], I[0,1], -I[0,0], -I[0,1]],
                     [I[1,0], I[1,1], -I[1,0], -I[1,1]],
                     [-I[0,0], -I[0,1], I[0,0], I[0,1]],
                     [-I[1,0], -I[1,1], I[1,0], I[1,1]]])
err = test - desired
print(sqrt(assemble(inner(err,err)*dx(domain=mesh))))

Testing this actually caught an index typo in desired[3,2], which is presumably what you’re trying to avoid with block matrix functionality!

3 Likes
3

Thank you, this is what I was looking for. I ended coming up with a similar (but sloppy) solution. Indeed, doing this manually is a typo minefield