Principle: P1 interpolation

If the domain is triangulated (divided into triangles which cover the domain , and do not intersect in other thing than a vertex or a whole side) the easiest is to approximate the function f by its "P1 interpolate". It is the unique function which

• is continuous,
• is linear within each triangle
• is equal to f at the vertices.

Proof

If a point x is within a triangle T with vertices then it can represented by its 3 "barycentric coordinates":

This is 3 equations for 3 unknowns lambda, the determinant of this linear system is twice the area of the triangle so it is non-zero and the lambdas are unique. Then the interpolate is

If x is on an edge, the interpolate one of the barycentric coordinate is zero so it is

The "level lines" of f are

When f is replaced by its interpolate, the level line are polygonal and on each triangle it is a line segment.

An intersection point q between a level line and an edge of the triangulation will be found if for some i,j,a we have:

So the principle for the drawing of a portion of level line is as follows.

for(int k = 0; k<nt; k++)
{ m=0;
  for(int i1=0;i1<3;i1++)
  {
    int j1 = (i1+1) % 3;
    int i = t[k].v[i1]->no, j = t[k].v[j1]->no;
    float a = (mu - f[j]) / (f[i] - f[j]);
    if((a>=0)&&(a<=1)) 
    {  m++; if(m==1) move2(a * q[i].x + (1-a) * q[j].x, 
                   a * q[i].y + (1-a) * q[j].y);
        else     line2(a * q[i].x + (1-a) * q[j].x, 
                   a * q[i].y + (1-a) * q[j].y);
    }
  }
}

It is necessary to check that the denominators are not zero and the drawing of the border of the domain must be added.

Write a program which reads a triangulation in the freefem format. Write a program which reads a file of values of f at the vertices and then generates a gnuplot file for the display of level lines. The presence of 2 parameters in the unix command must be tested. If there is only one parameter; if only one parameter is there, then the program will decide that it is a mesh plotting request and do it.