Dear Domenico,
Le 06/03/2015 21:03, Domenico Notaro a écrit :
>
> 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 ... ] ?
the gradient on a finite element nodes is consecutively in the fortran
order, so it corresponds to your second expression.
> 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?
I think, the faster should be the high level generic assembly and the
corresponding interpolation facilities.
>>
>> // 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)".
This is correct.
Yves.
--
Yves Renard ([email protected]) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
