[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Getfem-users] interpolated fem: gradient of the shape functions
From: |
Yves Renard |
Subject: |
Re: [Getfem-users] interpolated fem: gradient of the shape functions |
Date: |
Tue, 5 Jun 2007 17:16:35 +0200 |
User-agent: |
KMail/1.9.5 |
Le lundi 28 mai 2007 19:17, Andriy Andreykiv a écrit :
> Dear Yves and Julien,
>
> Would you be so kind to consult me on the following:
>
> I'm using the interpolated fem (new_interpolated_fem) in an FE
> model, based on fictitious domain method. For instance I was already
> able to successfully build a mass matrix between the mesh_fem's defined
> on two overlapping meshes that was used to impose Dirichlet boundary
> conditions in a week sense:
> ___________________________________________________________________________
>_____________________ getfem::mesh_fem mq_interpolate(line_mesh);
> //LINE_MESH='GT_PK(1,1)' getfem::pfem
> fem=getfem::new_interpolated_fem(mf_QuadMesh,mim);
> //MF_QUADMESH='FEM_QK(2,1)';
> mq_interpolate.set_finite_element(line_mesh.convex_index(),ifem);
> size_type n=mq_interpolate.nb_dof(),m=mf_line.nb_dof(); /MF_LINE =
> 'FEM_PK(1,1)';
> sparse_matrix M(n,m);
> getfem::generic_assembly assem;
> assem.set("M(#1,#2)+=comp(Base(#1).Base(#2))");
> assem.push_mi(mim);
> assem.push_mf(mq_interpolate);
> assem.push_mf(mf_line);
> assem.push_mat(M);
> assem.assembly();
> ___________________________________________________________________________
>______________________
>
> The above works perfect. However, now I want to interpolate the gradient
> of the shape functions from the quad mesh on the line mesh and calculate
> something like this:
>
> 1) assem.set("M(#1,#2)+=comp(Grad(#1).Base(#2))(:,i,:)");
> 2) assem.set("M(#1,#1)+=comp(Grad(#1).Grad(#1))(:,i,:,i)"); // with
> M(n,n)
> 3) assem.set("M(#1,#1)+=comp(vGrad(#1).vGrad(#1))(:,i,j,:,i,j)"); //
> with M(n,n)
>
> and so on
>
>
>
> In all cases I'm getting an error message:
>
>
> ===========================================================
>
> | An error has been detected !!!
>
> ===========================================================
> Error in getfem_mat_elem.cc, line 352:
> Internal error
>
> *** Exited with status: 1 ***
>
>
> Is this a bug? Or I'm doing something wrong? If this is a bug, could
> you, please, suggest a fix or a different way to obtain those matrices?
>
> Thank you in advance,
> Andriy
Yes, you are right, there is a bug in the computation of the gradient in
interpolated_fem. This simply has not been tested when the dimension of the
two fems are not the same.
Julien made the correction. you can get the corrected version of the file
getfem_interpolated_fem.cc on the gna.org web site (but, if you use the
version 2.2, you will get the version 3.0 which is a little bit different).
Otherwise, you can just replace the function
void interpolated_fem::real_grad_base_value
(const fem_interpolation_context& c, base_tensor &t, bool) const { ... }
in your version of getfem_interpolated_fem.cc
by the following one:
void interpolated_fem::real_grad_base_value
(const fem_interpolation_context& c, base_tensor &t, bool) const {
size_type N0 = mf.linked_mesh().dim();
elt_interpolation_data &e = elements.at(c.convex_num());
size_type nbdof = e.nb_dof, cv;
mi3[2] = N0; mi3[1] = target_dim(); mi3[0] = nbdof;
t.adjust_sizes(mi3);
std::fill(t.begin(), t.end(), scalar_type(0));
if (nbdof == 0) return;
if (c.have_pgp() &&
(&c.pgp()->get_point_tab() ==
&e.pim->approx_method()->integration_points()))
{
gausspt_interpolation_data &gpid = e.gausspt.at(c.ii());
if (gpid.iflags & 1) {
cv = gpid.elt;
pfem pf = mf.fem_of_element(cv);
if (gpid.iflags & 4) { t = gpid.grad_val; return; }
actualize_fictx(pf, cv, gpid.ptref);
pf->real_grad_base_value(fictx, taux);
if (pif) {
pif->grad(c.xreal(), trans);
for (size_type i = 0; i < pf->nb_dof(cv); ++i)
if (gpid.local_dof[i] != size_type(-1))
for (size_type j = 0; j < target_dim(); ++j)
for (size_type k = 0; k < N0; ++k) {
scalar_type ee(0);
for (size_type l = 0; l < N0; ++l)
ee += trans(l, k) * taux(i, j, l);
t(gpid.local_dof[i], j, k) = ee;
}
}
else {
for (size_type i = 0; i < pf->nb_dof(cv); ++i)
if (gpid.local_dof[i] != size_type(-1))
for (size_type j = 0; j < target_dim(); ++j)
for (size_type k = 0; k < N0; ++k)
t(gpid.local_dof[i], j, k) = taux(i, j, k);
if (store_values) { gpid.grad_val = t; gpid.iflags |= 4; }
}
}
}
else {
cerr << "NON PGP OU MAUVAIS PTS sz=" << elements.size() << ", cv=" <<
c.convex_num() << " ";
cerr << "ii=" << c.ii() << ", sz=" << e.gausspt.size() << "\n";
if (find_a_point(c.xreal(), ptref, cv)) {
pfem pf = mf.fem_of_element(cv);
actualize_fictx(pf, cv, ptref);
pf->real_grad_base_value(fictx, taux);
for (size_type i = 0; i < nbdof; ++i)
ind_dof.at(elements.at(cv).inddof.at(i)) = i;
if (pif) {
pif->grad(c.xreal(), trans);
for (size_type i = 0; i < pf->nb_dof(cv); ++i)
for (size_type j = 0; j < target_dim(); ++j)
for (size_type k = 0; k < N0; ++k)
if (ind_dof[mf.ind_dof_of_element(cv)[i]] != size_type(-1)) {
scalar_type ee(0);
for (size_type l = 0; l < N0; ++l)
ee += trans(l, k) * taux(i, j, l);
t(ind_dof[mf.ind_dof_of_element(cv)[i]],j,k) = ee;
}
}
else {
for (size_type i = 0; i < pf->nb_dof(cv); ++i)
for (size_type j = 0; j < target_dim(); ++j)
for (size_type k = 0; k < N0; ++k)
if (ind_dof[mf.ind_dof_of_element(cv)[i]] != size_type(-1))
t(ind_dof[mf.ind_dof_of_element(cv)[i]],j,k) = taux(i,j,k);
}
for (size_type i = 0; i < nbdof; ++i)
ind_dof[e.inddof[i]] = size_type(-1);
}
}
}
Could you please test this version and say us if it works now.
With kind regards,
Yves.
-------------------------------------------------------------------------
Yves Renard (address@hidden) tel : (33) 04.72.43.80.11
Pole de Mathematiques, INSA de Lyon fax : (33) 04.72.43.85.29
Institut Camille Jordan - CNRS UMR 5208
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
-------------------------------------------------------------------------
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- Re: [Getfem-users] interpolated fem: gradient of the shape functions,
Yves Renard <=