// pde_nuderiv2dg_tri.java
//                    second order pde with second degree solution
//                    just proof of concept
//                    generalize to fourth order, fourth degree
//                    sufficient points must be available for each order,degree
//
//  given a PDE, a triangularization, boundary values,
//                                    e.g. 22_ts.tri .bound .coord
//  Uxx(x,y)+2Uxy(x,y)+3Uyy(x,y)+4Ux(x,y)+5Uy(x,y)+6U(x,y)=f(x,y)
//  f(x,y)=6x^2+12xy+18y^2+42x+68y+101
//  uana(x,y)=x^2+2xy+3y^2+4x+5y+6   soultion, generally not known
//                                   used to compute boundary values
//                                   used for error checking
// U(xf,yf) unknown at  nfree  vertices, DOF
// U(xb,yb) known at  nbound  boundary vertices, nvert  total vertices
// nuderiv2d determines cxx, cxy, cyy, cx, cy at a subset of vertices
//  Uxx(xf,yf)= cxx(xf,yf)  *U(xf,yf)+        unknown U(xf,yf)   m<=36 points 
//              cxx(xfi,yfi)*U(xfi,yfi)+...+  unknown U(xfi,yfi) for each xf,yf
//              cxx(xbj,ybj)*U(xbj,ybj)+...   known value     
//  Uxy(xf,yf)= cxy(xf,yf)  *U(xf,yf)+        unknown U(xf,yf)
//              cxy(xfi,yfi)*U(xfi,yfi)+...+  unknown U(xfi,yfi)
//              cxy(xbj,ybj)*U(xbj,ybj)+...   known value 
//  Uyy(xf,yf)= cyy(xf,yf)  *U(xf,yf)+        unknown U(xf,yf)
//              cyy(xfi,yfi)*U(xfi,yfi)+...+  unknown U(xfi,yfi)
//              cyy(xbj,ybj)*U(xbi,ybi)+...   known value 
//  Ux (xf,yf)= cx (xf,yf)  *U(xf,yf)+        unknown U(xf,yf)
//              cx (xfi,yfi)*U(xfi,yfi)+...+  unknown U(xfi,yfi)
//              cx (xbj,ybj)*U(xbi,ybi)+...   known value 
//  Uy (xf,yf)= cy (xf,yf)  *U(xf,yf)+        unknown U(xf,yf)
//              cy (xfi,yfi)*U(xfi,yfi)+...+  unknown U(xfi,yfi)
//              cy (xbj,ybj)*U(xbi,ybi)+...   known value 
//  f(xf,yf)    right hand side, known, computable
//
// substituting in the PDE, a set of linear simultaneous equations for U(xf,yf)
//  | a11  a12  ...  a1n |   |U(xf1,yf1)| = rhs(1)  n = nfree
//  | a21  a22  ...  a2n | * |U(xf2,yf2)| = rhs(2)
//              ...
//  | an1  an2  ...  ann |   |U(xfn,yfn)| = rhs(n)
//
//       for cxx, cxy, ... at point xf1,yf1  a free vertex, a DOF
// a11 = cxx(xf1,yf1) + 2*cxy(xf1,yf1) + 3*cyy(xf1,yf1) +
//     4*cx (xf1,yf1) + 5*cy (xf1,yf1) + 6
//
// a1j = cxx(xfj,yfj) + 2*cxy(xfj,yfj) + 3*cyy(xfj,yfj) +
//     4*cx (xfj,yfj) + 5*cy (xfj,yfj) + 6            ip[1]=j  any free xfj,yfj
//
// rhs(1) = f(xf1,yf1) - cxx(xbj,ybj)*U(xbj,ybj) - 2cxy(xbj,ybj) ...
//                                                   ip[1]=j  any bound xbj,ybj
//
// same code for a21, a22, ... new N.nu2d on next point
// then rest of DOF points
// solve for U(xfi,yfi)  i=1..nfree

import java.io.*;  // for reading triangle data files  22_t*.tri .coord .bound

class pde_nuderiv2dg_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 
  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 

  pde_nuderiv2dg_tri()
  {
    int debug = 0;
    int m = 0;
    int j;
    int stat;
    int order = 2;
    double cx[] = new double[36]; // present limit of nuderiv2dg
    double cy[] = new double[36];
    double cxx[] = new double[36];
    double cxy[] = new double[36];
    double cyy[] = new double[36];
    int ip[] = new int[36];
    double ap;
    double error, maxerror;
    int kind = 1;

    System.out.println("pde_nuderiv2dg_tri.java running");
    System.out.println(
      "Uxx(x,y)+2Uxy(x,y)+3Uyy(x,y)+4Ux(x,y)+5Uy(x,y)+6U(x,y)=f(x,y)");
    System.out.println("f(x,y)=6x^2+12xy+18y^2+42x+68y+101");
    System.out.println("uana(x,y)=x^2+2xy+3y^2+4x+5y+6   solution");
    stat = do_input();
    if(stat != 0) System.exit(1); // can not continue 

    nuderiv2dg N = new nuderiv2dg();

    double x[] = new double[nvert];
    double y[] = new double[nvert];
    double U[] = new double[nvert];
    double A[][] = new double[nfree][nfree]; // matrix 
    double rhs[] = new double[nfree];    // right hand side, forcing function 
    double ug[] = new double[nfree];     // solution computed 
    double Ua[] = new double[nfree];      // known solution for checking 

    // initialize vertices
    System.out.println("vertex, x, y");
    for(int i=0; i<nvert; i++)
    {
      x[i] = xc[i]; 
      y[i] = yc[i];
      U[i] = 0.0;  // known at boundaries
      System.out.println(i+"  x="+x[i]+", y="+y[i]);
    }

    // clear A and rhs
    for(int i=0; i<nfree; i++)
    {
      rhs[i] = 0.0;
      for(int k=0; k<nfree; k++) A[i][k] = 0.0;
    }

    // set boundary values (as used in fem_check22p_tri.c) 
    for(int i=0; i<nbound; i++)
    {
      j=uniqueb[i];
      U[j]=uana(x[j],y[j]);
      System.out.println("bound "+i+"  vertex j="+j+" U[j]="+U[j]);
    }

    maxerror = 0.0;

    for(int point=0; point<nfree; point++) // for DOF 
    {
      m = N.nu2d(kind, order, nvert, point, x, y, cxx, cxy, cyy, cx, cy, ip);
      System.out.println("order="+order+", m="+m+", point="+point);
      if(debug>0)
      {
        System.out.println("i  ip  cxx  cxy  cyy   cx   cy");
        for(int i=0; i<m; i++)
        {
	  System.out.println(+i+" "+ip[i]+" "+cxx[i]+" "+cxy[i]+" "+
                             cyy[i]+" "+cx[i]+" "+cy[i]);
        }
      }
      // build row for point
      // compact code
      rhs[point] = f(x[point],y[point]);
      A[point][point]+=aii(cxx[0],cxy[0],cyy[0],cx[0],cy[0]);
      for(int i=1; i<m; i++)
      {
        ap=aij(cxx[i],cxy[i],cyy[i],cx[i],cy[i]);
	j=ip[i];
        if(not_in(j, uniqueb, nuniqueb)) // free vertex
	{
	  A[point][j] += ap;
          if(debug>0) 
            System.out.println("point="+point+", i="+i+", j="+j+", ap="+ap);
        }
	else // boundary
	{
	  rhs[point] -= ap*U[j];
          if(debug>0)
            System.out.println("point="+point+", i="+i+", j="+j+", ap*U="+
	 		       (-ap*U[j]));
	}
      } // end i across cxx, cxy, cyy, ...
    } // end point
    if(debug>0)
    {
      System.out.println("A matrix and rhs");
      for(int i=0; i<nfree; i++)
      {
        for(int k=0; k<nfree; k++)
	  System.out.println("A["+i+"]["+k+"]="+A[i][k]);
        System.out.println("rhs["+i+"]="+rhs[i]);
      }
      System.out.println(" ");
    }
    // solve linear equations
    new simeq(A, rhs, ug);
    for(int i=0; i<nfree; i++)
    {
      error = Math.abs(ug[i]-uana(x[i],y[i]));
      maxerror = Math.max(maxerror, error);
      System.out.println("vertex "+i+" x="+x[i]+",y="+y[i]);
      System.out.println("    computed="+ug[i]+", actual="+uana(x[i],y[i])+
                         ", error = "+error);
    }
    System.out.println("maxerror = "+maxerror);
    System.out.println(" ");
    System.out.println("end pde_nuderiv2dg_tri.java");
  }

  // PDE definition functions
  double f(double x, double y)
  {
    return 6.0*x*x + 12.0*x*y + 18.0*y*y + 
          42.0*x   + 68.0*y   + 101.0;
  }

  // boundary and test function
  double uana(double x, double y)
  {
    return x*x + 2.0*x*y + 3.0*y*y + 
         4.0*x + 5.0*y   + 6.0;
  }

  // PDE coefficients
  double aii(double axx, double axy, double ayy, double ax, double ay)
  {
    //     Uxx       Uxy       Uyy       Ux       Uy      U
    return axx + 2.0*axy + 3.0*ayy + 4.0*ax + 5.0*ay + 6.0;
  }
  // PDE derivatives
  double aij(double axx, double axy, double ayy, double ax, double ay)
  {
    //     Uxx       Uxy       Uyy       Ux       Uy
    return axx + 2.0*axy + 3.0*ayy + 4.0*ax + 5.0*ay;
  }

  // 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";
    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_ts";
    if(args.length>0) root=args[0];
    new pde_nuderiv2dg_tri();
  }
} // end  class pde_nuderiv2dg_tri