La matrice du système d'élasticité

A titre d'exemple nous donnons la fonction pour la construction de la matrice issue des équations de Lamé. Le principe est le même mais le calcul est plus complexe. Peut-être auront on avantage a définir des classes de vecteurs de fonctions chapeaux avec des opérateurs vectoriels pour simplifier ce genre de calcul

template <class T, class R> 
            void buildmatlame (Grid& g, float lam, float mu, float alpha,
                                                            Bandmatrix<T,R>& aa )
{      
      const int next[4] = {1,2,0,1};
      int i,j,k,kp,hp,ip,jp,ipp,jpp,iloc,jloc;
      float dwidxa,dwjdxa,dwidya,dwjdya;
      const float lammu = (lam+mu)*2, mu2 = mu*2, lam2 = lam*2, 
                  alph = alpha/12.0;;
      T aaloc;
      
      aa.zero();
      for ( k=0; k<g.nt; k++)
      {
            Triangle& tk = g.t[k];
            float tka = tk.area * 2;
            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<3; 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)/tka;
                         dwjdxa = (g.v[jp].y - g.v[jpp].y)/tka;
                        dwidya = -(g.v[ip].x - g.v[ipp].x)/tka;
                         dwjdya = -(g.v[jp].x - g.v[jpp].x)/tka;
            aaloc(0,0) = lammu * dwidxa * dwjdxa + mu2 * dwidya * dwjdya + alph;
            aaloc(1,1) = mu2 * dwidxa * dwjdxa + lammu * dwidya * dwjdya + alph;
            aaloc(1,0) = lam2 * dwidya * dwjdxa + mu2 * dwidxa * dwjdya;
            aaloc(0,1) = lam2 * dwidxa * dwjdya + mu2 * dwidya * dwjdxa;
                         if(i==j) for(kp=0;kp<=1;kp++) aaloc(kp,kp) += alph; 
                        aa(i,j) += aaloc * tk.area;
                  }
            }
      }      
}