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()