Aの作成

ポアソン方程式を例に説明します。

弱形式です。
P1要素で離散化します。要素K内では一次関数です。

とでもしておきます。

を満たすので、

ですね。三角形要素が面積0でない限りこの逆行列は存在します。この逆行列をDとおきます。すると

であることがわかります。これをもとの式に戻すと

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

とおいてあげることにします。つまり

となります。ここでλは面積座標になっています。
fも要素内で一次関数に近似しておきます。離散化された弱形式は

このように変形できます。具体的に見ていきましょう。

まずAはAkから次のように組み立てることができます。

int i,j,k,l;
for(i=0;i<noe;i++){
 
 /*Akをもとめる部分*/
 
 for(j=0;j<3;j++)
  for(k=0;k<3;k++)
   for(l=row_ptr[ebn[i][j]];l<row_ptr[ebn[i][j]+1];l++)
    if(col_ind[l]==enb[i][k]){
     val[l]+=Ak[j][k];//valは0で初期化しておいたので、+=でok!
     break;
    }
}
 

/*Akをもとめる部分*/ですが、式を見る限り各局所節点のx座標とy座標の値が分かればよさそうです。

double D[3][3];
for(j=0;j<3;j++)
 for(k=0;k<3;k++)
  D[j][k]=k==2?1:nc[en[i][j]][k];
 

このDの逆行列を求めたいと思います。手計算でもできますが、転記ミスのチェックが大変なので関数を作っておくとよいと思います。3*3専用なのでクラメルの公式で作ってみました。

void m_mat(double a[][3],double m[][2],int i,int j){//minnor matrix
	int k,l,g[3][2]={{1,2},{2,0},{0,1}};
	for(k=0;k<2;k++)
		for(l=0;l<2;l++)
			m[k][l]=a[g[i][k]][g[j][l]];
}
double det2(double a[][2]){
	return a[0][0]*a[1][1]-a[0][1]*a[1][0];
}
double det3(double a[][3]){
	int i;
	double tmp=0,m[2][2];
	for(i=0;i<3;i++){
		m_mat(a,m,i,0);
		tmp+=a[i][0]*det2(m);
	}
	return tmp;
}
double inv3(double A[][3],double invA[][3]){
	double M[2][2],detA=det3(A);
	int i,j;
	for(i=0;i<3;i++)
		for(j=0;j<3;j++){
			m_mat(A,M,i,j);
			invA[j][i]=det2(M)/detA;
		}
	return detA;
}
 

さて三角形要素の面積はdetD/2であることも加味すると、

int i,j,k,l;
for(i=0;i<noe;i++){
 
 double D[3][3];
 for(j=0;j<3;j++)
  for(k=0;k<3;k++)
   D[j][k]=k==2?1:nc[en[i][j]][k];
 
 double invD[3][3];
 double s;//面積
 s=inv3(D,invD)/2.;
 
 doubke Ak[3][3];//成分は零リセットされているはずです
 for(j=0;j<3;j++)
  for(k=0;k<3;k++){
   for(l=0;l<2;l++)
    Ak[j][k]+=invD[l][j]*invD[l][k];
   Ak[j][k]*=s;
  }
 
 for(j=0;j<3;j++)
  for(k=0;k<3;k++)
   for(l=row_ptr[ebn[i][j]];l<row_ptr[ebn[i][j]+1];l++)
    if(col_ind[l]==enb[i][k]){
     val[l]+=Ak[j][k];//valは0で初期化しておいたので、+=でok!
     break;
    }
}
 

Aが完成しました。