Dear Domenico,
Le 05/03/2015 00:17, Domenico Notaro a écrit : > > I am sorry, I missed something. > > > --- > *Domenico Notaro* > *mobile*: (+1) 412 983 0940 > *address*: 3765 Childs Street, Pittsburgh PA 15213 > *mail**to: *[email protected] > <mailto:[email protected]> > *Skype: *domenico.not > > ------------------------------------------------------------------------ > *Da:* Domenico Notaro > *Inviato:* giovedì 5 marzo 2015 00.13 > *A:* [email protected] > *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 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". > > 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. > > // 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)" 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
