Mesh transformation

Hi there,

Im trying to implement geometric parametrization in fenics. For this I need two domains, where one domain is a transformed version of the other (mesh2 = J*mesh1).

In order to compare the results that I would get from solving some differential equation on both domains, my idea was to use:

u_error.vector()[:] = u1.vector() - u2.vector()

But this does not work for the finite element that I use, which is ‘BDM’. However, it works with ‘CG’.

Here is an mwe, that demonstrates the issue. I would assume, u1 and u2 both to be (0,1) everywhere in their domain.

from dolfin import *
import numpy as np
from matplotlib import pyplot as plt

mesh1 = UnitSquareMesh(1,1)

Element = FiniteElement("BDM", mesh1.ufl_cell(), 1)
#Element = VectorElement("CG", mesh1.ufl_cell(), 1)

W1  = FunctionSpace(mesh1, Element)
g  = Constant((0,1))
u1  = interpolate(g,W1)

J = np.mat(  [[1,1], 
              [0,1]]  )

mesh2 = UnitSquareMesh(1,1)
mesh2.coordinates()[:] = (J*mesh2.coordinates().T).T

W2  = FunctionSpace(mesh2, Element)
u2  = Function(W2)
u2.vector()[:] = u1.vector()[:]

cb = plot(u1)

cb = plot(u2)

As BDM is an integral moment, using plot to visualize the function is not really a good idea. plot computes the vertex values of a function, which is a CG interpolation.
To visualize function with integral moments, you have to make your own plot routine.
To get started with this I suggest you have a look at the Brezzi-Douglas-Marini element in DefFem.