The first case is much easier to implement because the unknows at the center
of the triangles (bubbles) can be eliminated and the resulting linear system is of the form
In tnfem it is programmed exactly as the Lamé system except that N=3 instead of 2.
template <class T, class R>
void buildmatstokes (Grid& g, float nu, float alpha, float beta,
Bandmatrix<T,R>& aa )
{
const float alph = alpha/12.0, bet = beta*0.5/nu;
int i,j,k,kp,ip,jp,ipp,jpp,iloc,jloc;
float dwidxa,dwjdxa,dwidya,dwjdya;
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]]);
dwidxa = (g.v[ip].y - g.v[ipp].y)/tka;
dwidya = -(g.v[ip].x - g.v[ipp].x)/tka;
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]]);
dwjdxa = (g.v[jp].y - g.v[jpp].y)/tka;
dwjdya = -(g.v[jp].x - g.v[jpp].x)/tka;
aaloc(0,0) = nu*(dwidxa * dwjdxa + dwidya * dwjdya) + alph ;
aaloc(1,1) = nu*( dwidxa * dwjdxa + dwidya * dwjdya) + alph ;
aaloc(2,2) = bet*tk.area*( dwidxa * dwjdxa + dwidya * dwjdya) ;
aaloc(1,0) = aaloc(0,1) = 0;
aaloc(2,0) = aaloc(0,2) = dwidxa/3;
aaloc(2,1) = aaloc(1,2) = dwidya/3;
if(i==j) for(kp=0;kp<=1;kp++) aaloc(kp,kp) += alph;
aa(j,i) += aaloc * tk.area;
}
}
}
}