Re: [Getfem-users] advection term in getfem

[Yves Renard]

Dear Mattia Bozzola,

My answer was not complete since it depends on the form of your function F(u).
What I described in my previous email was the manner when it is suitable to 
describe F(u) on a finite element. If you do not want to interpole it on a 
finite element, you can also use what is called a nonlinear term in getfem. 
The use is a little bit more complicated because it needs to use more 
internal tools of Getfem. Here is an example of use.


template<typename VECT1> class your_nonlinear_term
        : public getfem::nonlinear_elem_term {
  
  const getfem::mesh_fem &mf;
  const VECT1 &U;
  size_type N;
  getfem::base_vector coeff, locU;
  getfem::base_matrix gradU;
  bgeot::multi_index sizes_;
  
public:
  your_nonlinear_term(const getfem::mesh_fem &mf_, const VECT1 &U_) 
    : mf(mf_), U(U_), N(mf_.linked_mesh().dim()), locU(mf_.get_qdim())
      gradU(mf_.get_qdim(), N) {
         if (mf_.get_qdim() == 1) { sizes_.resize(1); sizes_[0] = N; }
         else { sizes_.resize(2); sizes_[0] = mf_.get_qdim(); sizes_[1] = N; }
     }
  
  const bgeot::multi_index &sizes() const { return sizes_; }
  
  virtual void compute(getfem::fem_interpolation_context& ctx,
                                    bgeot::base_tensor &t) {
    size_type cv = ctx.convex_num();
    coeff.resize(mf.nb_dof_of_element(cv));
    gmm::copy(gmm::sub_vector(U, gmm::sub_index(mf.ind_dof_of_element(cv))),
              coeff);
    ctx.pf()->interpolation(ctx, coeff, locU, mf.get_qdim());
    ctx.pf()->interpolation_grad(ctx, coeff, gradU, mf.get_qdim());
    
   // ... your computation here. locU and gradU are the value of U and the
   // gradient of U on the current gauss point, respectively.
   // The result is to be put in t[i], 0 <= i <= N-1.    

  }
};

template<typename MAT, typename VECT>
  void asm_advection
  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
   const VECT &U, const mesh_region &rg = mesh_region::all_convexes()) {

    your_nonlinear_term<VECT> nl(mf, U);

     getfem::generic_assembly
       assem("M$1(#1,#1)+=comp(Base(#1).NonLin(#1).Grad(#1))(:,j,:,j)");
     assem.push_mi(mim);
     assem.push_mf(mf);
     assem.push_nonlinear_term(&nl);
     assem.push_mat(const_cast<MAT&>(M));
     assem.assembly(rg);
}


> > Dear professor,
> > I am an italian undergraduate student of  mathematics, and I am using
> > getfem in my thesis.
> > I would like to know if it is possible to assemble advection terms like
> > integral of (F(u)grad(u) wdx).
> > where u is the unknown, grad is the gradient, w is the test function and
> > F(u) is a function of u (that eventually I can treat in explicit...so it
> > can be considered simply as a coefficient).
> > I searched on the help file, but I was able to find only assembler for
> > mas matrix, stiffness matrix or assembler for spacific problems.
> > Looking foreward for your reply,
> > I thank you very much for your help.
> > Mattia Bozzola
>
> Here is an example of assembly procedure for an advection term. Vector A
> should contain the vectorial parameter you called F(u) traited explicitely.
>
>   template<typename MAT, typename VECT>
>   void asm_advection
>   (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
>    const mesh_fem &mf_data,
>    const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
>
>     generic_assembly
>       assem("a=data$1(#2); M$1(#1,#1)+="
>           "comp(Base(#1).vBase(#2).Grad(#1))(:,i,j,:,j).a(i)");
>     assem.push_mi(mim);
>     assem.push_mf(mf);
>     assem.push_mf(mf_data);
>     assem.push_data(A);
>     assem.push_mat(const_cast<MAT&>(M));
>     assem.assembly(rg);
>   }
>
>
> If you want to treat explicitely the term F(u) then you have either to make
                                   /\/\/\ implicitely of course
> iterations (fixed point strategy) or to use a Newton method. With the
> Newton method you have to compute the derivative of F(u) with respect to u.
>


Yves.


-------------------------------------------------------------------------
  Yves Renard (address@hidden)       tel : (33) 04.72.43.80.11
  Pole de Mathematiques,                       fax : (33) 04.72.43.85.29
  Institut Camille Jordan - CNRS UMR 5208
  INSA de Lyon, Universite de Lyon
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard
-------------------------------------------------------------------------