|Subject:||[Getfem-users] [ERRATA CORRIGE] Mixed Darcy problem: how assemble the "gradient term" and to compute the divergence of a vector field|
|Date:||Wed, 4 Mar 2015 23:17:02 +0000|
I am sorry, I missed something.
Da: Domenico Notaro
Inviato: giovedì 5 marzo 2015 00.13
Oggetto: Mixed Darcy problem: how assemble the "gradient term" and to compute the divergence of a vector field
Dear GetFem Users,
I am completely new of GetFem++ I am trying to implement a mixed formulation for the 3D Darcy problem.
By simplifying terms useless for this issue, the weak formulation of the problem isFind (u, p) in VxQ s.t.
(1) (1/k u, v) + (GRAD(p), v) = 0 in \Omega
(2) (u, GRAD(q)) - ((\alpha p, q)) + (f(p), q) = 0 in \Omega
where (.,.) and ((.,.)) indicate the L2 product on \Omega and \Gamma, respectively.
[I integrate by part the divergence term because of the Robin BC: u.n = \alpha p on \Gamma]
I have one question for each of the following step:
a) First of all, I tried to implement the assembly procedure for (GRAD(p), v).
b) Then I tried to evaluate the satisfaction of the divergence constraint (||DIV(u)-f||), i.e. the strong equivalent of (2), for which I need to compute the
divergence of a vector field.
a) The implementation is the following, It seems to work properly - I have just a small doubt below - but I would like a double check from you expert users because this is my first implementation:
/// Build the mixed pressure term
(?) The output of this asm procedure is G^T not G, isn't it? (because A(i,j)=a(\phi_j,\phi_i) is a(.,.) is a non symmetric bilinear form)
I mean, in this way I am assembling a mf_p.nb_dof() x mf_u.nb_dof() matrix while I need the transpose (to be multiplied then for the pressure vector P).
b) I tried to address this issue in two ways.
b.1) By using the function getfem::compute_gradient - that seems to be the only way to compute derivatives - I computed the gradient tensor of the vector velocity and then extracted the divergence:
// Compute GRAD(U)
//mf_U is at this level 'FEM_PK(3,1)'
// Compute DIV(U)
// Compute ||DIV(U)-F|| by using asm_L2_norm
(?) Here I assumed - but I am not sure at all - the function compute_gradient stores derivatives in the following order: GRAD(U) = [DxUx, DxUy, DxUz, DyUx, DyUy, ...]
b.2) In order to feel more confident about the previous implementation I tried also to compute ||DIV(u)-f|| with an assembly approach:
/// Compute the L2 norm of the residual of the divergence constraint
(?) The results are quite different from those of b.1 (also if the order is the same) so I don't know what to trust - if at least one is correct!
That's all! I am very sorry if it was too long and boring.
|[Prev in Thread]||Current Thread||[Next in Thread]|