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()[:]
plt.figure()
cb = plot(u1)
plt.savefig("u1.png")
plt.close()
plt.figure()
cb = plot(u2)
plt.savefig("u2.png")
plt.close()