La classe Problem

Un problème de résolution d'une EDP est défini par

• Un operateur, ici le Laplacien,
• une solution "sol" de type P1.
• un second membre: rhs,
• des conditions aux bords ici de Dirichlet, stockées dans "sol" en entrée,
• une méthode de résolution: ici le gradient conjugué
• les méthodes nécessaire au gradient conjugué: l'évaluation de la dérivée de lénergie par rapport aux degrés de liberté de la solution

class Laplace{ public:
  Grid g;
  P1 sol; // contains bdy cond on entry
  P1 rhs; // right hand side

  Laplace(const char* meshfile, const char* rhsfile, 
   const char* dirichletfile):
   g(meshfile), rhs(&g), sol(&g)
   { rhs.load(rhsfile), sol.load(dirichletfile);}
  void solvegradconj(const int niter, const float precise); 
  Vector derivener() const;
};

//-----------------------------------------------------
  void Laplace::solvegradconj(const int niter, const float eps)
{
	int i,m;
	float g2hm=1, g2h, eps1, ro, gamma, E1;
	P1 hconj(&g);
	Vector hconjm(g.nv), gradE(g.nv);
	
	for(m=0; m<niter; m++)
	{
		gradE = derivener();
		for(i=0;i<g.nv;i++) if(g.v[i].where) gradE[i] = 0;
        g2h = gradE.scal(gradE);
		if(m==0) { eps1 =  eps*g2h; gamma = 0;}
		  else  gamma = g2h / g2hm;	
  for(i=0;i<g.nv;i++) 
    hconj[i] = -gradE[i] + gamma * hconjm[i];
  ro = - gradE.scal(hconj) / hconj.intgrad2();
  cout << m <<"\t"<<ro<<"\t"<<g2h<<endl;
  if(g2h<eps1) break;
  for(i=0;i<g.nv;i++)
  {
    sol[i] +=  ro* hconj[i];
    hconjm[i] = hconj[i];
	 }
  g2hm = g2h;
 }
}

Vector Laplace::derivener() const
{
	int i,j,k,jloc;
	float x[3], y[3], phi[3];
	float gradFx,gradFy,gradGx, gradGy;
	Vector r(g.nv);
	for(i=0;i<g.nv;i++) r[i]=0;
	for (k=0; k<g.nt; k++)
	{
	  Triangle& tk = g.t[k];
	  float a1 = 0.5/ tk.area;
	  float f0 = rhs[g.no(tk.v[0])] + rhs[g.no(tk.v[1])] 
	             + rhs[g.no(tk.v[2])];
	  for (jloc=0; jloc<3; jloc++)
	  {
     j = g.no(tk.v[jloc]);
	  	 r[j] -= ( rhs[j] + f0 ) * (tk.area/12);
    	x[jloc] = g.v[j].x;  y[jloc]= g.v[j].y;                                           
    	phi[jloc] = sol[j];
   }
   gradFx = a1 * deter(phi[0],phi[1],phi[2],y[0],y[1],y[2]);
   gradFy = a1 * deter(x[0],x[1],x[2],phi[0],phi[1],phi[2]);
   for(i=0;i<3;i++)
   {  int i0= i==0, i1= i==1, i2= i==2,
     	gradGy = a1 * deter(x[0],x[1],x[2],i0,i1,i2);
     	gradGx = a1 * deter(i0,i1,i2,y[0],y[1],y[2]);
     	r[g.no(tk.v[i])] += tk.area * (gradFx*gradGx + gradFy*gradGy);  
	   }
	}
	return r;
}

void myexit();

void myexit()
{
  cout << "fin" << endl;
}

void main()
{
	atexit(myexit);
	Laplace p("circle.msh", "one.dta", "zero.dta");
	p.solvegradconj(1000, 1e-4);
	p.sol.save("sol.dta");
}