I See that the only

solution we have from what you have described is

recursively refine each triangle to 4 or 16 …

and draw the T3 triangles instead! Yes this is an

options but this is tricky

Just guessing here, but I think you are looking for something like this:

```
void DrawPolynomialSurface(
float P[3],
float DirX[3],
float DirY[3],
float DirZ[3],
float x,
float y,
float dx,
float dy,
float (*f)(float x,float y),
int SubdivisionDepth)
{
int i,j,k;
glBegin(GL_TRIANGLES);
for(i=0;i<SubdivisionDepth;i++)
for(j=0;j<SubdivisionDepth;j++)
{
float a0=i/(float)SubdivisionDepth;
float b0=j/(float)SubdivisionDepth;
float a1=(i+1)/(float)SubdivisionDepth;
float b1=(j+1)/(float)SubdivisionDepth;
float x0=x+a0*dx;
float y0=y+b0*dy;
float x1=x+a1*dx;
float y1=y+b1*dy;
float Disp00[3];//Disp00=P+a0*DirX+b0*DirY+f(x00,y00)*DirZ;
float Disp01[3];//Disp01=P+a0*DirX+b1*DirY+f(x01,y01)*DirZ;
float Disp10[3];//Disp10=P+a1*DirX+b0*DirY+f(x10,y10)*DirZ;
float Disp11[3];//Disp11=P+a1*DirX+b1*DirY+f(x11,y11)*DirZ;
for(k=0;k<3;k++)
{
Disp00[k]=P[k]+a0*DirX[k]+b0*DirY[k]+f(x0,y0)*DirZ[k];
Disp01[k]=P[k]+a0*DirX[k]+b1*DirY[k]+f(x0,y1)*DirZ[k];
Disp10[k]=P[k]+a1*DirX[k]+b0*DirY[k]+f(x1,y0)*DirZ[k];
Disp11[k]=P[k]+a1*DirX[k]+b1*DirY[k]+f(x1,y1)*DirZ[k];
}
glVertex3fv(Disp00);
glVertex3fv(Disp01);
glVertex3fv(Disp10);
glVertex3fv(Disp01);
glVertex3fv(Disp11);
glVertex3fv(Disp10);
}
glEnd();
}
// Calling it like this should give you a nice surface.
float f(float x,float y)
{
return x*x*x+y*y;
}
//...
glMatrixMode(GL_PROJECTION);
glOrtho(-5.0f,5.0f,-5.0f,5.0f,-20.0f,20.0f);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
//...
float P[3]={0.0f,0.0f,0.0f};
float DirX[3]={+1.0f,1.0f,1.0f};
float DirY[3]={-1.0f,1.0f,1.0f};
float DirZ[3]={0.0f,1.0f,-1.0f};
DrawPolynomialSurface(P,DirX,DirY,DirZ,-1.0f,-1.0f,2.0f,2.0f,f,10);
```

Note: this is illustrative (read: really slow) code.