Fonctions P1

Une fonction P1 est définie par ses valeurs aux sommets de la triangulation. Le coeur de la classe P1 sera donc une sous classe de la classe Vector

Les fonctions P1 sont susceptibles d'autres opérations:

• Lecture et écriture au format freefem,
• integration, dérivation,
• tracé graphique...

Mais compte tenu de la suite nous n'introduisons que les fonctions de lecture d'écriture sur fihier, et d'integration pour calculer:

float deter(const float a,const float b,const float c,
        const float d,const float e,const float f) 
      { float g = (a-c)*(e-f) - (b-c)*(d-f);
          return g;
      } 
/***************************************************/
class P1: public Vector
{ 
public:
  Grid* g;
  P1(Grid* gg): g(gg), Vector(gg->nv){}
  void load(const char* filename);
  void save(const char* filename) const;
  float intgrad2() const; //int gradf.gradf dx
};

//---------------------------------------------------------
void P1::save(const char* path) const
{  // writes a P1 continuous function in freefem format
 	ofstream file(path);
 	int nv = size;
 	assert(!file.fail());
 	file << nv << endl;
 	for(int i=0; i<nv; i++) file << cc[i] << endl;
 	file.close();
}

//---------------------------------------------------------
void P1::load(const char* path)
{  // read a P1 continuous function in freefem format
 	ifstream file(path);
 	int nv ;
 	assert(!file.fail());
 	file >> nv;
 	assert(nv==size);
 	for(int i=0; i<nv; i++) file >> cc[i];
 	file.close();
}

//-----------------------------------------------------
float P1::intgrad2() const
{
 float x[3], y[3], phi[3],gradFx, gradFy, ee;
 int k, iloc, j;

  ee = 0.;
  for(k=0;k<g->nt;k++)	
  { float a1 = 0.5/g->t[k].area;
    for(int jloc=0;jloc<3;jloc++)
    {   j = g->no(g->t[k].v[jloc]);
    	x[jloc] = g->v[j].x;  y[jloc]= g->v[j].y;                                           
    	phi[jloc] = cc[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]);
    ee +=  g->t[k].area * (gradFx*gradFx + gradFy*gradFy);  
   }
 return ee;
}