La fonction rhs()

Pour calculer l'intégrale de f(x,y) sur le domaine on boucle sur les triangles et on utilise une formule de Gauss à 3 points, les milieux des arètes.

template <class T> Vector<T> rhs(Grid& g, const Vector<T>& f)
{
      
      const int next[4] = {1,2,0,1};
      const T tzero(0);
      int k;
      Vector<T> r(g.nv);
      for (k=0; k<g.nv; k++) r[k] = tzero;
      for (k=0; k<g.nt; k++)
      {
            Triangle& tk = g.t[k];
            int i = g.no(tk.v[0]);
            int ip = g.no(tk.v[1]);
            int ipp =  g.no(tk.v[2]);
            T f0 = f[i] + f[ip] + f[ipp];
            for (int iloc=0; iloc<3; iloc++)
                  r[i] += ( f[g.no(tk.v[iloc])] + f0 ) *  (tk.area/12) ;
      }
      return r;
}