Possibly wrong dofmap for mixed function space of parametric curves

mesh
1

Assume we have a parametric curve \boldsymbol{r}(u) in 3D

\boldsymbol{r}(u)=\begin{pmatrix}x(u)\\ y(u)\\ z(u) \end{pmatrix},\quad u\in[0,1]

The parameter u is called the material coordinate and it is associated with material points on the curve \boldsymbol{r}(u).

The weak form of the problem we are working on includes different vector-valued functions \boldsymbol{r}(u), \boldsymbol{v}(u), … .

When we create a mixed-function space for the vector-valued functions, there seems to be a problem with how the dofmap is created. To be precise, the values for the coordinates y and z of the first two vertices of \boldsymbol{r}(u) are flipped. Here is a MWE

from fenics import * import dolfin print(f'Dolfin version: {dolfin.__version__}\n')

N = 11
mesh = UnitIntervalMesh(N-1)

P1 = FiniteElement(‘Lagrange’, mesh.ufl_cell(), 1)
P1_3 = MixedElement([P1] * 3)
V3 = FunctionSpace(mesh, P1_3)
W = FunctionSpace(mesh, MixedElement([P1_3, P1_3]))

r = Expression((‘x[0]’, ‘x[0]’, ‘x[0]’), degree = 1)
v = Expression((‘0’, ‘0’, ‘0’), degree = 1)

r = interpolate(r, V3)
v = interpolate(v, V3)

dof_map0 = r.function_space().sub(0).dofmap().dofs()
dof_map1 = r.function_space().sub(1).dofmap().dofs()
dof_map2 = r.function_space().sub(2).dofmap().dofs()

print(‘Vector r from function space V3’)
print(f’x = {r.vector().get_local()[dof_map0]}‘)
print(f’y = {r.vector().get_local()[dof_map1]}’)
print(f’z = {r.vector().get_local()[dof_map2]} \n’)

fa = FunctionAssigner(W, [V3, V3])

u = Function(W)

fa.assign(u, [r, v])

r, v = u.split()

dof_map0 = r.function_space().sub(0).dofmap().dofs()
dof_map1 = r.function_space().sub(1).dofmap().dofs()
dof_map2 = r.function_space().sub(2).dofmap().dofs()

print(‘Vector r from mixed function space W’)
print(f’x = {r.vector().get_local()[dof_map0]}‘)
print(f’y = {r.vector().get_local()[dof_map1]}’)
print(f’z = {r.vector().get_local()[dof_map2]}')

On my machine, this creates the output.

Dolfin version: 2019.1.0

Vector r from function space V3

x = [1. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]
y = [1. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]
z = [1. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]

Vector r from mixed function space W
x = [1. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]
y = [0.9 1. 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]
z = [0.9 1. 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0. ]

Is this a bug or are we doing something wrong here? Thanks for taking the time!

2

Hi,

If you enable deepcopy with r, v = u.split(True), then you get the correct dofmap.

3

Hi Jesus, thank you very much for looking into this that fast. I can confirm that r, v = u.split(True) works on my machine. Can you explain why we need a deepcopy here? Is it in general better to use split(True) ?

4

Hi Lukas.

What happens is that at the beginning, the functions r and v go with the dofmap of the function space V3. But when you create the function u, it goes with the dofmap of the space W. By using u.split() (without the deepcopy), you create shallow copies of subfunctions, populating it with references to the child objects found in the original, this includes the ordering of the function u (which is the order of W). On the other hand, if you use deepcopy=True, you are making recursive copies of the subfunctions, along with their ordering. Since you used assign, this ordering is copied as well. In this case, anything you do on r or v will not affect u.

On the other hand, you can achieve the same ordering by using shallow copies via

u = Function(W)
r1, v1 = u.split()
r1.assign(r)
v1.assign(v)

However, anything that happens over r1 or v1 will affect u. I particularly tend to use deepcopy in most of my simulations.

Everything I tell you now is from experience. Maybe a Fenics developer or someone with more knowledge of this can clarify your doubt much better.

Greetings.