// chladni.edp --- solve eigenvalue problems of the biharmonic operator
//
// 参考: http://ljll.math.upmc.fr/pipermail/freefempp/2012-March/001737.html
// 使用法: FreeFem++ chladni.edp project_name [number of ev]
//   FreeFem++ chladni.edp 6gon80 とすると、
//     領域のメッシュデータ 6gon80.msh を読み、
//     200個の固有値、固有関数を計算して
//     6gon80/ev.dat, 6gon80/ef{0,1,...,199}.eps, 6gon80/ef{0,1,...,199}.dat
//   を出力する
//   FreeFem++ chladni.edp 6gon80 100 とすると、固有値は100個のみ計算
//   FreeFem++ chladni.edp 6gon80 200 -nw とすると、図は描かない(バッチむき)
//   当面、正多角形の固有値・固有関数の計算に用いるが、機能は一般的
//
//   試しに
//     gcc -o polygon polygon.c
//     ./polygon 6 20 6gon20.msh
//     mkdir 6gon20
//     FreeFem++ chladni.edp 6gon20

load "Morley"

verbosity=1;
int n, NN, i;
int nev=200;
string projectname, dirname;

// ポアソン比はアルミニウムの値
real mu=0.345;

if (ARGV.n >= 3) {
  projectname = ARGV[2];
  if (ARGV.n >= 4)
    nev = atoi(ARGV[3]);
  cout << "projectname: " + projectname + ", nev=" + nev << endl;
}
else {
  cout << "usage: " + ARGV[0] + " " + ARGV[1] + " meshfile_name [nev [-nw]]" << endl;
  exit(0);
}
// ディレクトリ作成
dirname=projectname;
exec("mkdir " + dirname);

// 三角形分割
mesh Th=readmesh(projectname+".msh");
plot(Th,wait=false,ps=dirname+"/mesh.eps");

fespace Vh(Th, P2Morley);
Vh [u,ux,uy], [v,vx,vy];

varf J([u,ux,uy], [v,vx,vy]) = int2d(Th)((dxx(u)+dyy(u))*(dxx(v)+dyy(v))
  -(1-mu)*(dxx(u)*dyy(v)+dyy(u)*dxx(v)-2.0*dxy(u)*dxy(v)));
varf K([u,ux,uy], [v,vx,vy]) = int2d(Th)(u*v); 

matrix A=J(Vh,Vh,solver=UMFPACK); 
matrix B=K(Vh,Vh,solver=UMFPACK); 

real[int] ev(nev); // Stockage des valeurs propres
Vh[int] [eVu,eVux,eVuy](nev);  // Stockage des vecteurs propres

int k=EigenValue(A,B,sym=true,value=ev,vector=eVu,tol=1e-10,maxit=0,ncv=0);

ofstream f(dirname+"/ev.dat");
// 15桁表示する
f.precision(15);
for (i=0; i< nev; i++)
  f << ev[i] << endl;

real[int] levels=[0.0];

for (i=0;i<nev;i++) {
  cout << ev[i] << endl;
  ofstream f(dirname+"/ef"+i+".dat");
  f << eVu[i][] << eVux[i][] << eVuy[i][] << endl;

  // ,viso=levels の代わりに ,nbiso=64 のようにするとレベルはお任せ
  // ,fill={false,true} で塗るか塗らないか
  plot(eVu[i],viso=levels,cmm=i+"th eigenfunction, eigenvalue="+ev[i],
       wait=false,ps=dirname+"/ef"+i+".eps", bw=true);
}