I am sorry, I missed something.
--- Domenico Notaro mobile: (+1) 412 983 0940 address: 3765 Childs Street, Pittsburgh PA 15213 mailto: [email protected]<mailto:[email protected]> Skype: domenico.not <mailto:[email protected]> ________________________________ 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). 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); // 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. Thank you, Domenico user: domenico_notaro
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
