Many thanks for your answer. I may not have all the terms of the static
condensation in your solution, as i made some mistakes. Sorry.
My email is way too long! So the short question is:
In order to have nodal values of the solution vector, i've got to work
with the type "std::vector<Tensor<1,dim>>" which is difficult to handle
with the "FullMatrix<double>" of the extracted submatrix of the system
matrix at the element level, as i have to compute [Kpu]{u}. Is it
possible to transform the std::vector<Tensor<1,dim>> to the
std::vector<double> type ?
The detailed but i hope comprehensive explanation of why i want that:
The way i compute the static condensation of the pressure at the element
level is as follow:
*) the weak form is 2*a(grad v, grad u)-a(div v,p)-a(phi,div
u)+a(phi,\lambda p) = a(v,f), lambda being small
*) I create the system matrix [K] at the element level.
*) I extract [Kuu], [Kpu], [Kup], and [Kpp] as FullMatrix, and condense
the pressure (I create A=-[Kup][Kpp]^-1[Kpu]) and scatter the values to
the velocity dof)
Then, i would like to compute the elementary residual at each cell
<R(u),v>, that is to say compute
For each cell, {R(u)} = 2*a(grad v, grad u) - a(v,f) - [Kup][Kpp]^-1[Kpu]{u}
So i've got to, {u} given:
*) compute {R(u)}=2*a(grad v, grad u) - a(v,f) at gauss points
*) compute [Kup][Kpp]^-1[Kpu]{u} at solution node values
*) subtract the latter to take into account the condensed terms, then
distribute_local_to_global
*) compute the L2 error of the resulting vector.
The problem is that i can have the nodal values of {u} with the
following code:
const std::vector<Point<dim>> support_points =
fe.base_element(0).get_unit_support_points();
std::vector<double> weights(support_points.size(),1);
const Quadrature<dim> node_quadrature_formula(support_points,weights);
FEValues<dim> fe_node_values (fe, node_quadrature_formula,
date_values | update_gradients |
date_quadrature_points | update_JxW_values);
But then i will work with the type "std::vector<Tensor<1,dim>>" which is
difficult to handle with the "FullMatrix" of the extracted submatrix, as
i have to compute [Kpu]{u}. Then how can i compute [Kpu]{u} ? Is it
possible to transform the std::vector<Tensor<1,dim>> to the
std::vector<double> type ? Is it the way to do the static condensation
at the element level ? (step-44 does it but not at the element level)
Do that make more sense ? It could still be confusing.
Many thanks,
Best regards,
Yann
Le 15/10/2024 à 19:04, Wolfgang Bangerth a écrit :
On 10/15/24 01:19, Yann Jobic wrote:
My next step is to compute that way (static condensation of the
pressure) the navier-stokes equations with a Newton non linear solver,
as in step-57, but i do it on a simpler problem first (stokes).
Thus if i define
{R(u,p)} = {-2*div u + \nabla p -f},
I want to compute ||R(u,p)||.
Yann,
so the residual for you is a function of x, and what you're looking to
compute is an integral (for the norm). Integrals in deal.II are always
computed using FEValues, and what you want here is going to look
something like the following:
double integral = 0;
for (cell=...)
{
fe_values.reinit(cell);
fe_values[velocities].get_function_divergences(solution, div_u);
fe_values[pressure].get_function_gradients(solution, grad_p);
rhs.value_list(fe_values.quadrature_points(), f_values);
for(q=...)
integral += std::pow(-2*div_u[q] + grad_p[q] - f_values[q], 2) *
fe_values.JxW(q);
}
double norm = std::sqrt(integral);
(Your formula for the residual isn't correct because div u is a scalar
and you're adding grad p to it, which is a vector. This can't be
right. I'll let you fix this, and just wrote the code matching your
formula.)
You might want to take a look at step-15 to see some of the parts I
left out above.
Best
W.
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/6caa1836-256a-45c0-a39c-196c877e97c7%40gmail.com.
