bの作成

さてAの作成で出てきた

ですが、面積座標の積分公式

を使うと

であることがわかります。またFは

であるのでコードは次のようになります。

double *b=calloc(non,sizeof(double));//零リセットしておきます
double F[3];
for(j=0;j<3;i++)
 F[i]=f(nbc[ebn[i][j]][0],nbc[ebn[i][j]][1]);
double sum;//一時的に使う変数です
double bk[3];
for(j=0;j<3;j++){
 sum=0;
 for(k=0;k<3;k++)
  sum+=s*(j==k?1/12.:1/6.)*F[k];
 bk[j]=sum;
}
for(j=0;j<3;j++)
 b[ebn[i][j]]+=bk[j];