|
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 (address@hidden) 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
---------
|