getfem-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Mixed Darcy problem: how to assemble the "gradient te


From: Domenico Notaro
Subject: Re: [Getfem-users] Mixed Darcy problem: how to assemble the "gradient term" and to compute the divergence of a vector field
Date: Fri, 6 Mar 2015 20:03:36 +0000

Dear Prof. Renard,



thank you so much for your precious support.

Just a few extra question:

 
Dear Domenico,


Le 05/03/2015 00:17, Domenico Notaro a √©crit :

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 is

  Find (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
/// $ G = \int GRAD(p).v dx $
template<typename MAT>
void
asm_darcy_G(MAT &                G,
            const mesh_im  &    mim,
            const mesh_fem &    mf_p,
            const mesh_fem &    mf_u,
            const mesh_region & rg = mesh_region::all_convexes()
            )
{
    GMM_ASSERT1(mf_p.get_qdim() == 1, "invalid data mesh fem (Qdim=1 required)");
    GMM_ASSERT1(mf_u.get_qdim() > 1,  "invalid data mesh fem (Qdim=2,3 required)");
    generic_assembly assem("M(#1,#2)+=comp(Grad(#1).vBase(#2))(:,i,:,i);");
    
    assem.push_mi(mim);
    assem.push_mf(mf_p);
    assem.push_mf(mf_u);
    assem.push_mat(G);
    assem.assembly(rg);

} /* end of asm_darcy_G */


(?) 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).


Yes, this is correct, and this is indeed a mf_p.nb_dof() x mf_u.nb_dof() matrix. You can also use the high level generic assembly to perform this with an assembly string of the type "Grad_Test_p.Test2_u".

OK, that's great!


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)
    getfem::mesh_fem mf_gradU(mesh);
    bgeot::pgeometric_trans pgt_t = bgeot::geometric_trans_descriptor(MESH_TYPE);
    size_type N = pgt_t->dim();
    mf_gradU.set_qdim(bgeot::dim_type(N), bgeot::dim_type(N)); //3x3
    //mf_gradUt.set_classical_finite_element(0);
    mf_gradUt.set_classical_discontinuous_finite_element(0);
    vector_type gradU(mf_gradU.nb_dof());
    getfem::compute_gradient(mf_U, mf_gradU, U, gradU);

    //mf_U is at this level 'FEM_PK(3,1)'


    // Compute DIV(U)
    getfem::mesh_fem mf_Ui(mesh);
    mf_Ui.set_classical_discontinuous_finite_element(0);
    size_type nb_dof_Ui = mf_Ui.nb_dof(); //= mf_gradU.nb_dof()/(N*N)! NOT mf_U.nb_dof()/N
    vector_type divU(nb_dof_Ui);    
    
    gmm::add(gmm::sub_vector(gradU, gmm::sub_interval(0*nb_dof_Ui, nb_dof_Ui)), divU);
    gmm::add(gmm::sub_vector(gradU, gmm::sub_interval(4*nb_dof_Ui, nb_dof_Ui)), divU);
    gmm::add(gmm::sub_vector(gradU, gmm::sub_interval(8*nb_dof_Ui, nb_dof_Ui)), divU);

This seems to be incorrect. the component of the gradient are consecutives in the vector gradU. So you should use a gmm::sub_slice instead.

OK, I will try to fix that. So, is the correct order
  GRAD(U) = [DxUx,1  DyUx,1 DzUx,1  DxUx,2  DyUx,2  DzUx,2 ... ]
or 
  GRAD(U) = [DxUx,1  DyUx,1 DzUx,1  DxUy,1  DyUy,1  DzUy,1 ... ] ?
I mean, also the components are stored in this "sliced" way?
And, in general, which approach do you think is better between the "algebraic" and the "assembly" ones?


    // 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
/// $ ||DIV(u) - f|| = sqrt( \int (DIV(u) - f)^2 dx ) $
template<typename VEC>
scalar_type
asm_div_error_L2_norm( const VEC &U, const mesh_fem &mf_u,
                       const VEC &F, const mesh_fem &mf_f,
                       const mesh_im &mim,
                       const mesh_region &rg = mesh_region::all_convexes() )
{
    GMM_ASSERT1(mf_u.get_qdim() > 1, "invalid data mesh fem (Qdim>1 required)");
    GMM_ASSERT1(mf_f.get_qdim() == 1, "invalid data mesh fem (Qdim=1 required)");
    GMM_ASSERT1(U.size() == mf_u.nb_dof(), "invalid vector data (size=mf.nb_dof() required)");
    GMM_ASSERT1(F.size() == mf_f.nb_dof(), "invalid vector data (size=mf.nb_dof() required)");
    GMM_ASSERT1(U.size() == (F.size()*mf_u.get_qdim()), "invalid vector data (U.size=Qdim*F.size required)" );

    generic_assembly
    assem("u=data$1(#1);"
          "f=data$2(#2);"
          "V()+=u(i).u(j).comp(vGrad(#1).vGrad(#1))(i,k,k,j,h,h)"
              "-u(i).f(j).comp(vGrad(#1).Base(#2))(i,k,k,j)"
              "-f(i).u(j).comp(Base(#2).vGrad(#1))(i,j,k,k)"
              "+f(i).f(j).comp(Base(#2).Base(#2))(i,j);");
    assem.push_mi(mim);
    assem.push_mf(mf_u);
    assem.push_mf(mf_f);
    assem.push_data(U);
    assem.push_data(F);
    std::vector<scalar_type> v(1);
    assem.push_vec(v);
    assem.assembly(rg);

    return sqrt(v[0]);
}


(?) 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.



You can also use here the high level generic assembly instead of the low-level one, this should be simpler. You can use the assembly string "sqr(Trace(Grad_u) - f)"

OK, I will try that approach too!
Apart from that, do you think this implementation should work?
I am not sure about the index reduction:
"V()+=u(i).u(j).comp(vGrad(#1).vGrad(#1))(i,k,k,j,h,h)".

Yves.

Many thanks.
Best,
Domenico
user id: domenico_notaro
---

Domenico Notaro

mobile:  (+1) 412 983 0940
address: 3765 Childs Street, Pittsburgh PA 15213
mailto:  address@hidden
Skype:   domenico.not



reply via email to

[Prev in Thread] Current Thread [Next in Thread]