From state discretization to control discretization

General
1

Dear FEniCS Community,

I am currently solving optimal control problems (OCP) with PDEs as constraints, in particular the Poisson equation.
I initially wanted to use the adjoint to get some information for the control but there I encountered a weird phenomenon.
The adjoint lies in the same discretization as the state (on nodes), which is not the same as the one from the control (on cells).
Now I wanted to bring the adjoint to the same discretization as the control.

For example if we have 4x4 unit square mesh, then we have 32 cells (triangles) and 25 nodes, so 16 values for the controls and 25 values for the adjoint and state.
So my first idea (Variant A) was to just calculate the mean of the 3 nodes of each triangle to get a constant value for each triangle.
The second idea (Variant B) was to “cut” out 2 arbitrary boarders to get from 25 values to 16 values, which can easily be doubled to 32 values. This actually results in squares instead of triangles, but doesn’t make a difference for the sake of this phenomenon.

I tried this out and got some weird pictures for the calclulated controls.
I found out that the same happens for just the state, so I made a little example showing this phenomenon without the OCP or adjoint and using the state instead of the adjoint for this phenomenon. We look at the problem:

\Delta y = f \text{ in } \Omega
y = 0 \text{ on } \Omega .

Here y is the state and f is the (initial) control from the OCP.
In the code I created functionspaces Y for states and U for controls.
The Poisson equation is then solved and afterwards the state is transformed via variant A and B to the control space.
The expected result should look like the solution y, but in both cases the result look different if compared to the solution y.
And that’s my question and phenomenon: why does this happen and (how) can it be fixed?

Code 
from dolfin import *
import numpy as np

# Create Mesh
n = 64
mesh = UnitSquareMesh(n, n)

# Functionspaces, Functions, Boundary
Y = FunctionSpace(mesh, 'CG', 1)
U = FunctionSpace(mesh, 'DG', 0)
y_ = TrialFunction(Y)
v_ = TestFunction(Y)
f = interpolate(Expression("1", degree=1), U)
bc = DirichletBC(Y, 0.0, "on_boundary")

# Solve Weak Formulation
A = dot(grad(y_), grad(v_)) * dx
L = f * v_ * dx
y = Function(Y)
solve(A == L, y, bc)
File("y.pvd") << y

# Transformation Variant A (Triangles)
u_A = Function(U)
u_A_vec = np.zeros(len(u_A.vector()[:]))
for cell in cells(mesh):
    p_value = 0
    for vertex in vertices(cell):
        p_value += y.vector()[vertex.index()]
    p_value /= 1/3
    u_A_vec[cell.index()] = p_value

u_A.vector().set_local(u_A_vec)
File("u_A.pvd") << u_A

# Transfomration Varaint B (Cut off 2 boarders)
y_cut = y.vector()[:].reshape(n+1, n+1)
y_cut = y_cut[0:-1, 0:-1]
u_B = Function(U)
u_B.vector().set_local(np.repeat(y_cut, 2))
File("u_B.pvd") << u_B

My guess would be that it has something to do with the shape function and it has somehow to be inverted(?). But I am not quite sure if it’s true and my understanding is not too deep in that area.
If it’s necessary I can show the Code for the OCP as well, which uses dolfin_adjoint.
Here are the pictures I get from the Code:

Solution y Variant A Variant B
Solution y Variant A Variant B

If compared, the solution y looks different from Variant A and B, which raised my question. It’s not surprising that Variant A and B look the same.
I hope my explanation of the question is clear.

Thanks for any help or suggestion!

2

First of all, I do not understand what information you are trying to gain by putting the adjoint into the control space.

Secondly, you need to use vertex_to_dofmap to map the vertex index to the correct dof index, i.e.

from dolfin import *
import numpy as np

# Create Mesh
n = 64
mesh = UnitSquareMesh(n, n)

# Functionspaces, Functions, Boundary
Y = FunctionSpace(mesh, 'CG', 1)
U = FunctionSpace(mesh, 'DG', 0)
y_ = TrialFunction(Y)
v_ = TestFunction(Y)
f = interpolate(Expression("1", degree=1), U)
bc = DirichletBC(Y, 0.0, "on_boundary")

# Solve Weak Formulation
A = dot(grad(y_), grad(v_)) * dx
L = f * v_ * dx
y = Function(Y)
solve(A == L, y, bc)
File("y.pvd") << y

# Transformation Variant A (Triangles)
u_A = Function(U)
u_A_vec = np.zeros(len(u_A.vector()[:]))
v_to_d = vertex_to_dof_map(Y)
for cell in cells(mesh):
    p_value = 0
    for vertex in vertices(cell):
        p_value += y.vector()[v_to_d[vertex.index()]]
    p_value /= 1/3
    u_A_vec[cell.index()] = p_value

u_A.vector().set_local(u_A_vec)
File("u_A.pvd") << u_A
3

Thanks dokken!
That’s exactly what I needed.