La fonction buidmatlaplace()

La construction de la matrice se fait, comme d'habitude, par assemblage par triangle. Pour les problèmes vectoriel le type T sera MatN et le type R sera VectN. Le cas n'est pas inclus ici pour plus de clarté. Le programme qui suit calcule donc l'intégrale sur le domaine du produit des gradients des fonctions chapeaux associées aux sommets de la triangulation. On ne tient pas compte des conditions aux limites dans ce calcul. Celle-ci seront prisent en compte plus tard par pénalisation.

template <class T, class R> 
      void buildmatlaplace (Grid& g, Bandmatrix<T,R>& aa, T alpha )
{
      
      const int next[4] = {1,2,0,1};
      T alph;
      int i,j,k,ip,jp,ipp,jpp,iloc,jloc;
      T dwidxa,dwjdxa,dwidya,dwjdya,aaloc;
      for ( k=0; k<aa.size*aa.bdth; k++) aa.cc[k] = 0.F;

      for ( k=0; k<g.nt; k++)
      {
            Triangle tk = g.t[k];
            alph = alpha * (tk.area / 12) ;
            for ( iloc=0; iloc<3; iloc++)
            {
                  i = g.no(tk.v[iloc]);
                  ip = g.no(tk.v[next[iloc]]);
                  ipp = g.no(tk.v[next[iloc+1]]);
                  for ( jloc=0; jloc<=iloc; jloc++)
                  {
                         j = g.no(tk.v[jloc]);                   
                         jp = g.no(tk.v[ next[jloc]]);                   
                         jpp = g.no(tk.v[next[jloc+1]]);
                         dwidxa = (g.v[ip].y - g.v[ipp].y)/(tk.area * 4); 
                         dwjdxa = g.v[jp].y - g.v[jpp].y;
                         dwidya = -(g.v[ip].x - g.v[ipp].x)/(tk.area * 4);
                         dwjdya = -(g.v[jp].x - g.v[jpp].x);
                  
                         aaloc =  dwidxa * dwjdxa + dwidya * dwjdya + alph;
                        aa(i,j) += aaloc;
                        if(i != j) aa(j,i) += aaloc; else aa(i,i) += alph;
                  }
            }
      }      
}