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

10 months ago by
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$ƒ ind g [P0]3 s.t. Ωg·ψ=Ωu·ψψP1    

where V is a discontinuous Lagrange space P0. This approach is used also in: .
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 10 months ago by Miguel de Benito  

1 Answer

10 months ago by
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)$g=(g1,g2,g3) , then solve three problems

 $\int_{\Omega}g_i\phi=\int_{\Omega}\frac{\partial u}{\partial x_i}\phi,\qquad\forall\phi\in P_1$Ωgiϕ=Ωuxi ϕ, ∀ϕP1 
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)$u(x)=uiψ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 L2 projection.
written 10 months ago by Francesco Clerici  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »