Re: [Getfem-users] interpolated fem: gradient of the shape functions

[Yves Renard]

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
-------------------------------------------------------------------------