Gradient of a "CG" order one element.

Question

Hi all,

I am running an iterative scheme which includes the gradient of a function. Therefore, I use the gradient operator in Fenics, but this not returning a Fenics function, which makes sense because the gradient will not always live in some finite element space. However, in the case of piece-wise linear continuous functions the derivative would be just piece-wise constant.
I know that I can just use the projection into the "DG" 0 (originally DG1, that was a typo) space which should give me the result and to some extend I can improve the speed of this operation a lot. Now to my question, is it possible to compute the derivative in this case without taking the workaround of projecting the gradient into the space of piece-wise constants? This is performance wise and in terms of memory usage quite inconvenient for me, since the scales of my problems will be quite big. Basically I have to solve a second huge linear system which gives me no new information, computing the gradient of something piece-wise linear is a lot easier problem than that. So far I am not really using this fact. Does anyone have any idea to get around this problem?

Best,
Sebastian

Answer 1

I asked a similar question, some years ago

https://fenicsproject.org/qa/1899/fast-projection-in-case-of-explicit-relations

Comments

Hm. I am familiar with this speed up option. However, it does not really solve the problem because I still have storage of the assembled matrices which produces a huge memory overhead. But I will read through it maybe I can get some additional improvement. Actually, I was hoping for something more elegant around this issue.

Edit:
I actually tried the code and there is no attribute 'compress()' for matrices anymore, was it removed?

Answer 2

Hi, the following uses the fact that for DG0 space(a natural space for a derivative of CG1) the mass matrix is diagonal with the entries being volumes of the cells. Thus you can actually assemble the vector representing the gradient directly and then use it to define function.

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(20, 20, 20)
V = FunctionSpace(mesh, 'CG', 1)

f = Function(V)
f_vec = f.vector()
f_vec.set_local(np.random.rand(f_vec.local_size()))
f_vec.apply('insert')

Q = VectorFunctionSpace(mesh, 'DG', 0)
# Grad via projection
p, q = TrialFunction(Q), TestFunction(Q)
L = inner(grad(f), q)*dx
M = assemble(inner(p, q)*dx)
b = assemble(L)

grad_f0 = Function(Q)
x0 = grad_f0.vector()
solve(M, x0, b)

# Assemble the action of Minv onto the rhs - no mass matrix needed
MinvL = (1/CellVolume(mesh))*inner(grad(f), q)*dx
x = assemble(MinvL)

print (x - x0).norm('linf')

# You can define a function now
grad_f = Function(Q, x) 

Why do you want to represent the gradient in DG1?

Comments

Oh sorry, this is a typo :/ Clearly I want to be in DG0. I will try this approach, thank you very much.

---

Okay that is definitely what I was looking for, but I have another question related to this and I am not sure if it is worth opening a new topic.
Will the gradient be exact with this method up to machine precision or is there some numerical round off? Additionally, is there a way to influence the assemble command in performance and accuracy wise?

---

For precision consider the following

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(20, 20, 20)
V = FunctionSpace(mesh, 'CG', 1)

f = interpolate(Expression('2*x[0]+3*x[1]+x[2]', degree=1), V)

Q = VectorFunctionSpace(mesh, 'DG', 0)
q = TestFunction(Q)

MinvL = (1/CellVolume(mesh))*inner(grad(f), q)*dx
x = assemble(MinvL)
grad_f = Function(Q, x)

grad_f0 = interpolate(Constant((2, 3, 1)), Q)

print (grad_f.vector() - grad_f0.vector()).norm('linf')