// test_nuderiv2dg2_tri.java   test  nuderiv2dg2  with 28 points 
//                    order 1  minimum points 3       
//                    order 2  minimum points 6       
//                    order 3  minimum points 10      
//                    order 4  minimum points 15      
//                    order 5  minimum points 28      
//                    going beyond minimum points does not help much 
import java.io.*;

class test_nuderiv2dg2_tri
{
  // variables for reading triangles
  static String   root;             // root name for .tri, .coord, .bound 
  final int sz = 512;
  final int szsz = 1024;
  int debug = 0;
  int    ntri;                      // number of triangles  
  int    minvert;                   // minimum vertex number, reduced to zero 
  int    t1[] = new int[sz];        // vertex numbers on triangles 
  int    t2[] = new int[sz];        // vertex numbers on triangles 
  int    t3[] = new int[sz];        // vertex numbers on triangles 
  int    ncoord;                   // number of coordinates, also equals nvert 
  double xc[] = new double[sz];     // coordinates in vertex order 
  double yc[] = new double[sz];     // coordinates in vertex order 
  int    b1[] = new int[sz];        // dirichlet boundary edge vertices 
  int    b2[] = new int[sz];        // dirichlet boundary edge vertices 
  double xmin, xmax, ymin, ymax;    // determined by input 
  double kg[] = new double[szsz];   // global stiffness matrix 
  double fg[] = new double[sz];     // right hand side, forcing function 
  double ug[] = new double[sz];     // solution computed by fem 
  double Ua[] = new double[sz];     // known solution for checking 
  int uniquev[] = new int[sz];      // all vertices 
  int nvert;                        // number of unique vertices 
  int uniqueb[] = new int[sz];      // bound vertices 
  int nbound;                       // number of boundary vertices 
  int uniquef[] = new int[sz];      // free vertices, 0 to nfree-1 
  int nfree;                        // degrees of freedom  nvert-nbound 
  int nuniquev, nuniqueb, nuniquef; // vertices counts 

  test_nuderiv2dg2_tri()
  {
    int debug = 0;
    int i, j, n;
    int m = 0;
    int order = 3;
    int stat;
    final int siz = 1024;
    double x[] = new double[siz];
    double y[] = new double[siz];
    double U[] = new double[siz];
    double est, actual, dx, dy;
    int ip[] = new int[36];
    double cx[] = new double[36];
    double cy[] = new double[36];
    double cxx[] = new double[36];
    double cxy[] = new double[36];
    double cyy[] = new double[36];
    double cxxx[] = new double[36];
    double cxxy[] = new double[36];
    double cxyy[] = new double[36];
    double cyyy[] = new double[36];
    double error, maxerror;

    nuderiv2dg N = new nuderiv2dg();

    System.out.println("test_nuderiv2dg2_tri.java running");
    stat = do_input();
    if(stat != 0) System.exit(1); // can not continue 

    n = ncoord;
    // initialize vertices
    for(i=0; i<n; i++)
    {
      x[i] = xc[i]; 
      y[i] = yc[i];
    }

    // set boundary values (as used in fem_check22p_tri.c) 
    for(i=0; i<nbound; i++)
    {
      j=uniqueb[i];
      // U[j]=f(xc[j],yc[j]);
    }

    for(i=0; i<n; i++)
    {
      U[i] = f(x[i],y[i]);
      dx = fx(x[i],y[i]);
      dy = fy(x[i],y[i]);
      System.out.println("point "+i+"  x="+x[i]+", y="+y[i]+
                         ", U="+U[i]+", dx="+dx+", dy="+dy);
    } 

    maxerror = 0.0;

    for(int point=0; point<n; point++) // Ux 
    {
      m = N.nu2dx(order, n, point, x, y, cx, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
      actual = fx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Ux 

    for(int point=0; point<n; point++) // Uy 
    {
      m = N.nu2dy(order, n, point, x, y, cy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
      actual = fy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uy 

    for(int point=0; point<n; point++) // Uxx 
    {
      m = N.nu2dxx(order, n, point, x, y, cxx, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
      actual = fxx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uxx 

    for(int point=0; point<n; point++) // Uxy 
    {
      m = N.nu2dxy(order, n, point, x, y, cxy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
      actual = fxy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uxy 

    for(int point=0; point<n; point++) // Uyy 
    {
      m = N.nu2dyy(order, n, point, x, y, cyy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
      actual = fyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uyy 


    for(int point=0; point<n; point++) // Uxxx 
    {
      m = N.nu2dxxx(order, n, point, x, y, cxxx, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
      actual = fxxx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxxx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uxxx 

    for(int point=0; point<n; point++) // Uxxy 
    {
      m = N.nu2dxxy(order, n, point, x, y, cxxy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
      actual = fxxy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxxy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uxxy 

    for(int point=0; point<n; point++) // Uxyy 
    {
      m = N.nu2dxyy(order, n, point, x, y, cxyy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
      actual = fxyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uxyy 

    for(int point=0; point<n; point++) // Uyyy 
    {
      m = N.nu2dyyy(order, n, point, x, y, cyyy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      est = 0.0;
      for(i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
      actual = fyyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cyyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
    } // end Uyyy 

    System.out.println("maxerror for order="+order+", m="+m+
                         ", is "+maxerror);

    // now, get all at once, for first four points
    order = 3;
    int kind = 1;
    for(int point=0; point<4; point++) // all 
    {
	m = N.nu2d3(kind, order, n, point, x, y, cxxx, cxxy, cxyy, cyyy,
                    cxx, cxy, cyy, cx, cy, ip);
      System.out.println("point="+point+", kind="+kind+
                         " cxxx, cxyy, cxyy, cyyy,cxx,cxy,cyy,cx,cy");
      est = 0.0;
      for(i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
      actual = fxxx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxxx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
      actual = fxxy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxxy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
      actual = fxyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
      actual = fyyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cyyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
      actual = fxx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
      actual = fxy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cxy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
      actual = fyy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cyy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
      actual = fx(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cx "+point+", est="+est+", actual="+actual+
                         ", error="+ error);
      est = 0.0;
      for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
      actual = fy(x[point],y[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      System.out.println("cy "+point+", est="+est+", actual="+actual+
                         ", error="+ error);

    }
      System.out.println("maxerror for order="+order+", m="+m+
                         ", is "+maxerror);
  } // end test_nuderiv2gg_tri

  // test functions
  double f(double x, double y)
  {
    return x*x*x + 2.0*x*x*y + 3.0*x*y*y + 4.0*y*y*y +
             x*x +   2.0*x*y + 3.0*y*y + 
           4.0*x +     5.0*y + 6.0;
  }

  double fx(double x, double y)
  {
    return 3.0*x*x + 4.0*x*y + 3.0*y*y +
           2.0*x   + 2.0*y   + 4.0;
  }

  double fy(double x, double y)
  {
    return 2.0*x*x + 6.0*x*y + 12.0*y*y +
             6.0*y +   2.0*x + 5.0;
  }

  double fxx(double x, double y)
  {
    return 6.0*x + 4.0*y + 
           2.0;
  }

  double fxy(double x, double y)
  {
    return 4.0*x + 6.0*y +
	   2.0;
  }

  double fyy(double x, double y)
  {
    return 6.0*x + 24.0*y +
           6.0;
  }

  double fxxx(double x, double y)
  {
    return 6.0;
  }

  double fxxy(double x, double y)
  {
    return 4.0;
  }

  double fxyy(double x, double y)
  {
    return 6.0;
  }

  double fyyy(double x, double y)
  {
    return 24.0;
  }

  // routines for reading triangle data
  int uniqueint(int nold, int a[])
  {
    int i, j;
    if(debug>2) System.out.println("finding unique of nold="+nold);
    sortint(nold,a);
    j=0;
    for(i=1; i<nold; i++)
    {
      if(a[i]>a[i-1])
      {
        j++;
        a[j]=a[i];
      }
    }
    j++;
    if(debug>2) System.out.println("found unique of n="+j);
    return j;
  } // end uniqueint 

  void sortint(int n, int a[])
  {
    int i, j, temp;
    for(i=0; i<n-1; i++)
      for(j=i+1; j<n; j++)
        if(a[i]>a[j]) {temp=a[i]; a[i]=a[j]; a[j]=temp;}
  } // end sortint 

  int freevert(int c[], int na, int a[], int nb, int b[])
  {
    // a and b sorted integer arrays, c will be from a, not in b 
    int ia=0;
    int ib=0;
    int ic=0;
    if(debug>2) System.out.println("freevert na="+na+", nb="+nb);
    while(true)
    {
      if(ib>=nb || a[ia]!=b[ib]) {c[ic]=a[ia]; ic++; ia++;}
      else             {ia++; ib++;}
      if(ia>=na) break;
    }
    if(debug>2) System.out.println("freevert nc free = "+ic);
    return ic;
  }

  boolean not_in(int i, int a[], int n)
  {
    boolean found = false;
    int j;

    for(j=0; j<n; j++)
    {
      if(i==a[j]) { found=true; break; }
    }
    return !found;
  } // end not_in 

  int do_input() // three files root .tri  root.bound  root.coord
  {
    String fname;
    int index, last, len;
    double tmp;

    // read in triangles, set ntri 
    fname = root+".tri";
    if(debug>2) System.out.println("about to read triangles from "+fname);
    ntri=0;
    minvert=999999;
    nvert=0;

    try
    {
      BufferedReader in = new BufferedReader(new FileReader(fname));
      String input_line;
      input_line = in.readLine();
      while(input_line != null)
      {
	if(debug>0) System.out.println(input_line);
        // have a line, extract 3 integers
        len = input_line.length();
        if(len<5)  // not enough data, ignore blank lines
	{
	  input_line = in.readLine();
          continue;
        }
        index = 0;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        t1[ntri] = Integer.parseInt(input_line.substring(index,last)); 
        index = last+1;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        t2[ntri] = Integer.parseInt(input_line.substring(index,last)); 
        index = last+1;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        if(last<1) last = len;
        t3[ntri] = Integer.parseInt(input_line.substring(index,last)); 

        if(debug>0) System.out.println("tri "+ntri+" has vertices "+t1[ntri]+
                                        "  "+t2[ntri]+"  "+t3[ntri]);
        if(t1[ntri]>nvert) nvert=t1[ntri];
        if(t2[ntri]>nvert) nvert=t2[ntri];
        if(t3[ntri]>nvert) nvert=t3[ntri];
        if(t1[ntri]<minvert) minvert=t1[ntri];
        if(t2[ntri]<minvert) minvert=t2[ntri];
        if(t3[ntri]<minvert) minvert=t3[ntri];
        ntri++;
	input_line = in.readLine();
      }
    }
    catch(IOException exception)
    {
      System.out.println("can not open for reading "+fname);
      System.out.println(exception);
      return 1;
    }

    if(debug>0) System.out.println(ntri+" triangles read from "+fname);

    nvert=nvert+1-minvert; // still must be dense sequence of numbers 
    nuniquev=0;
    for(int i=0; i<ntri; i++)
    {
      t1[i]=t1[i]-minvert;
      t2[i]=t2[i]-minvert;
      t3[i]=t3[i]-minvert;
      uniquev[nuniquev]=t1[i];
      uniquev[nuniquev+1]=t2[i];
      uniquev[nuniquev+2]=t3[i];
      nuniquev=nuniquev+3;
    }
    nuniquev=uniqueint(nuniquev, uniquev);

    // read in dirichlet boundary, set nbound and nuniqueb 
    fname = root+".bound";
    if(debug>2) System.out.println("about to read boundary from "+fname);
    nbound=0;

    try
    {
      BufferedReader in = new BufferedReader(new FileReader(fname));
      String input_line;
      input_line = in.readLine();
      while(input_line != null)
      {
        // have a line, extract 2 integers
        len = input_line.length();
        if(len<4)  // not enough data, ignore blank lines
	{
	  input_line = in.readLine();
          continue;
        }
        index = 0;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        b1[nbound] = Integer.parseInt(input_line.substring(index,last)); 
        index = last+1;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        if(last<1) last = len;
        b2[nbound] = Integer.parseInt(input_line.substring(index,last)); 
        if(debug>0) System.out.println("boundary segment nbound="+nbound+
                              " has vertices "+b1[nbound]+", "+b2[nbound]);
        if(b1[nbound]>nvert) System.out.println("bad boundary vertex "+
                                                b1[nbound]);
        if(b2[nbound]>nvert) System.out.println("bad boundary vertex "+
                                                b2[nbound]);
        if(b1[nbound]<minvert) System.out.println("bad boundary vertex "+
                                                  b1[nbound]);
        if(b2[nbound]<minvert) System.out.println("bad boundary vertex "+
                                                  b2[nbound]);
        nbound++;
	input_line = in.readLine();
      }
    }
    catch(IOException exception)
    {
      System.out.println("can not open for reading "+fname);
      System.out.println(exception);
      return 1;
    }
    if(debug>0) System.out.println("read Dirichlet boundaries from "+fname);

    nuniqueb=0;
    for(int i=0; i<nbound; i++)
    {
      b1[i]=b1[i]-minvert;
      b2[i]=b2[i]-minvert;
      uniqueb[nuniqueb]=b1[i];
      uniqueb[nuniqueb+1]=b2[i];
      nuniqueb=nuniqueb+2;
    }
    nuniqueb=uniqueint(nuniqueb, uniqueb);
    nuniquef=freevert(uniquef, nuniquev, uniquev, nuniqueb, uniqueb);
    if(debug>0)
      for(int i=0; i<nuniqueb; i++)
        System.out.println("unique boundary "+uniqueb[i]);
    if(debug>0) System.out.println("number of free vertices is "+nuniquef);
    if(debug>0) 
      for(int i=0; i<nuniquef; i++) 
        System.out.println("free vertex "+uniquef[i]);

    // read in coordinates, check consistency 
    fname = root+".coord";
    if(debug>2) System.out.println("about to read coordinates from "+ fname);

    ncoord=0;
    try
    {
      BufferedReader in = new BufferedReader(new FileReader(fname));
      String input_line;
      input_line = in.readLine();
      while(input_line != null)
      {
        // have a line, extract 2 double
        len = input_line.length();
        if(len<5)  // not enough data, ignore blank lines
	{
	  input_line = in.readLine();
          continue;
        }
        index = 0;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        xc[ncoord] = Double.parseDouble(input_line.substring(index,last)); 
        index = last+1;
        while(input_line.substring(index,index+1).equals(" ")){index++;}
        last = input_line.indexOf(' ',index+1);
        if(last<1) last = len;
        yc[ncoord] = Double.parseDouble(input_line.substring(index,last)); 
        if(debug>0) System.out.println("coordinate "+ncoord+" at  "+
                                       xc[ncoord]+", "+yc[ncoord]);
        if(ncoord==0)
        {
          xmax=xc[ncoord];
          xmin=xc[ncoord];
          ymax=yc[ncoord];
          ymin=yc[ncoord];
        }
        if(xc[ncoord]>xmax) xmax=xc[ncoord];
        if(xc[ncoord]<xmin) xmin=xc[ncoord];
        if(yc[ncoord]>ymax) ymax=yc[ncoord];
        if(yc[ncoord]<ymin) ymin=yc[ncoord];
        ncoord++;
	input_line = in.readLine();
      }
    }
    catch(IOException exception)
    {
      System.out.println("can not open for reading "+fname);
      System.out.println(exception);
      return 1;
    }
    if(debug>0) System.out.println("coordinates read from "+fname);


    if(ncoord!=nvert)
       System.out.println("**** wrong number of coordinates="+ncoord+
                          ", should be"+nvert);
    nfree=nvert-nuniqueb;
    // make free vertices first
    if(debug>0)
    {
      System.out.println("vertex, x, y  before free first");
      for(int i=0; i<nvert; i++)
	  System.out.println(i+"  x="+xc[i]+", y="+yc[i]);
    }
    for(int i=0; i<nfree; i++)
    {
      if(not_in(i, uniqueb, nuniqueb)) continue; // OK so far
      // find a free vertex to swap
      for(int j=i+1; j<nvert; j++)
      {
        if(not_in(j, uniqueb, nuniqueb))
	{
	  // swap coord
          tmp = xc[i]; xc[i]=xc[j]; xc[j]=tmp;
          tmp = yc[i]; yc[i]=yc[j]; yc[j]=tmp;
          // swap bound
          for(int k=0; k<nuniqueb; k++)
	  {
	     if(uniqueb[k]==i)
	     {
               uniqueb[k]=j;
	     } 
	     else if(uniqueb[k]==j)
	     {
	       uniqueb[k]=i;
	     } 
	  }
          // swap tri
          for(int k=0; k<ntri; k++)
          {
	     if(t1[k]==i)
	     {
               t1[k]=j;
	     } 
	     else if(t1[k]==j)
	     {
	       t1[k]=i;
	     } 
	     if(t2[k]==i)
	     {
               t2[k]=j;
	     } 
	     else if(t2[k]==j)
	     {
	       t2[k]=i;
	     } 
	     if(t3[k]==i)
	     {
               t3[k]=j;
	     } 
	     else if(t3[k]==j)
	     {
	       t3[k]=i;
	     } 
          }
	} // end "if" no swap
      } // end j of free first
    } // end i of free first

    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    System.out.println("nvert="+nvert+", nbound="+nbound+", nuniqueb="+
                        nuniqueb+", nfree="+nfree+", ntri="+ntri);
    return 0;
  } // end do_input 


  public static void main (String[] args)
  {
    root = "22_tss";
    if(args.length>0) root=args[0];
    new test_nuderiv2dg2_tri();
  }
} // end  class test_nuderiv2dg2_tri