### How to compute the gradient of a scalar P1 function in C++ without solving a variational problem?

273

views

0

Hi, I need to compute the gradient of a scalar P1 function, say u. At the moment I compute it by solving the following linear variational problem:

$find\quad\textbf{g }\in\left[P_0\right]^3\quad s.t.\quad\int_{\Omega}\textbf{g}\cdot\nabla\psi=\int_{\Omega}\nabla u\cdot\nabla\psi\quad\quad\quad\quad\quad\forall\psi\in P_1$

where V is a discontinuous Lagrange space P0. This approach is used also in: https://fenicsproject.org/qa/10947/computing-derivatives-in-c .

The fact is that solving this problem is a lot more expensive than solving the original scalar problem in which we compute u. So that, when the mesh is thick, the program goes out of memory when solving the above problem.

I tried to use finite difference (first order since the function is P1), but I don't know how to size the increment h without going out of the cell. In order to avoid this, one should use h = Cell::inradius(), and centering the finite difference on the incenter of the cell (and I don't know how to compute the coordinates of the incenter of a tetrahedron).

Do you have any suggestions?

$find\quad\textbf{g }\in\left[P_0\right]^3\quad s.t.\quad\int_{\Omega}\textbf{g}\cdot\nabla\psi=\int_{\Omega}\nabla u\cdot\nabla\psi\quad\quad\quad\quad\quad\forall\psi\in P_1$

`ƒ``i``n``d`**g**∈[`P`_{0}]^{3}`s`.`t`. ∫_{Ω}**g**·∇`ψ`=∫_{Ω}∇`u`·∇`ψ`∀`ψ`∈`P`_{1}where V is a discontinuous Lagrange space P0. This approach is used also in: https://fenicsproject.org/qa/10947/computing-derivatives-in-c .

The fact is that solving this problem is a lot more expensive than solving the original scalar problem in which we compute u. So that, when the mesh is thick, the program goes out of memory when solving the above problem.

I tried to use finite difference (first order since the function is P1), but I don't know how to size the increment h without going out of the cell. In order to avoid this, one should use h = Cell::inradius(), and centering the finite difference on the incenter of the cell (and I don't know how to compute the coordinates of the incenter of a tetrahedron).

Do you have any suggestions?

Community: FEniCS Project

You might find this question useful.

written
8 months ago by
Miguel de Benito

### 1 Answer

1

What do you want to do with the gradient ? Your gradient problem does not make sense to me. I would have done an L2 projection. This however needs solving a matrix problem. If $g=\left(g_1,g_2,g_3\right)$

$\int_{\Omega}g_i\phi=\int_{\Omega}\frac{\partial u}{\partial x_i}\phi,\qquad\forall\phi\in P_1$∫

You have to solve mass matrix for which you can use CG.

Alternately, you can search for gradient recovery techniques. For P1, there is a simple procedure which involves averaging gradient from all elements that meet at a vertex. This does not need any matrix solves.
Thank you for the answer. What I want to compute is the following quantity:

$\nabla u\left(x\right)=\sum u_i\nabla\psi_i\left(x\right)$∇

Hypothetically, I need a function defined in the following way:

I heard about the gradient recovery techniques, but I don't need a reconstruction, I need the L

`g`=(`g`_{1},`g`_{2},`g`_{3}) , then solve three problems$\int_{\Omega}g_i\phi=\int_{\Omega}\frac{\partial u}{\partial x_i}\phi,\qquad\forall\phi\in P_1$∫

_{Ω}`g`_{i}`ϕ`=∫_{Ω}∂`u`∂`x`_{i}`ϕ`, ∀`ϕ`∈`P`_{1}You have to solve mass matrix for which you can use CG.

Alternately, you can search for gradient recovery techniques. For P1, there is a simple procedure which involves averaging gradient from all elements that meet at a vertex. This does not need any matrix solves.

$\nabla u\left(x\right)=\sum u_i\nabla\psi_i\left(x\right)$∇

`u`(

`x`)=∑

`u`

_{i}∇

`ψ`

_{i}(

`x`)

Hypothetically, I need a function defined in the following way:

`double compute_gradient_from_component(Function u, int cell, int component)`

which takes the function to be derived, the cell index where I need the value and the component (x, y or z) and return the value of the gradient on a cell. Consider the mesh given.I heard about the gradient recovery techniques, but I don't need a reconstruction, I need the L

^{2}projection.

written
8 months ago by
Francesco Clerici

Please login to add an answer/comment or follow this question.