Wrong solution to trivial variational form in Mixed Function Space

General
1

Hello,

I was attempting to implement a time stepping scheme, but kept getting weird solutions. After some trial and error I reduced the problem to the following:

from dolfin import *
import numpy as np

file = File('data/data.pvd')

mesh = UnitSquareMesh(100, 100)

element = {
            "f": FiniteElement("CG", mesh.ufl_cell(), 2),
            "sigma": VectorElement("CG", mesh.ufl_cell(), 2),
        }

X = FunctionSpace(mesh, MixedElement(*element.values()))

initial_f = Expression('exp(-(pow(x[0]-0.5, 2) + pow(x[1]-0.5, 2))/0.1)', degree = 2)
U0 = Function(X)
f0 = project(initial_f, X.sub(0).collapse())
assign(U0.sub(0), f0)
f0, sigma = split(U0)
file << U0.split(deepcopy=True)[0]


bc = DirichletBC(X.sub(0), Constant(0), "on_boundary")

U = Function(X)
f, sigma = split(U)
phi = TestFunctions(X)
v, tau = phi

# Define the variational form
a = f * v *dx - f0 * v * dx + inner(sigma, tau) * dx
# This should be solved by f = f0 and sigma = 0

solve(a == 0, U, bcs = bc)

f0, sigma0 = split(U)
solutions_f = U.sub(0, deepcopy=True).vector()
solutions_s = U.sub(1, deepcopy=True).vector()

# Update U0
U0.vector().set_local(np.concatenate((solutions_f, solutions_s)))
f0, sigma0 = split(U0)

file << U0.split(deepcopy=True)[0]

As expected a solution is found in a single iteration. I would expect the solution to the variational problem to be f = f0 and sigma = 0. However, plotting the solution in Paraview results in this:

Screenshot from 2024-04-16 15-53-07
Screenshot from 2024-04-16 15-53-07573×572 129 KB

I am using version 2019.2.0.dev0 of dolfin, installed using the Ubuntu ppa. I should also note that this does not happen, when I use a regular Function Space. This only happens when I introduce the mixed function space.

2

Hi.

Welcome to the FEniCS community!

You should not use

U0.vector().set_local(np.concatenate((solutions_f, solutions_s)))

to update functions because the DOF’s ordering might be wrong. To properly update a Function you can use

U0.assign(U)

Cheers.

2 Likes