// fem_check22p_tri.c  Galerkin FEM on triangular grid CM phi does not work
//                     triangle sides parallel to axes
// solve  uxx(x,y) + uyy(x,y) = f(x,y)
// boundary conditions computed using u(x,y)
// analytic solution may be given by  u(x,y) = x*x+2*x*y+3*y*y+4*x+5*y+6
// f(x,y) = 8
// input triangles and coordinates
// input geometry, root= 9t  .coord,  .tri,  .bound
// and others
// 
// solve for Uij numerical approximation U(x_i,y_j) of u(x_i,y_j)
// K * U = F, solve simultaneous equations for U

import java.io.*;
 
class fem_check22p_tri
{
  final int sz = 50;
  final int sz3 = 150;
  final int szsz = 2500;
  final int np = 3;
  final int szint = np*np;

  int    debug=3;                   // debug print level 
  static String   root;             // root name for .tri, .coord, .bound 
  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 
  double cm1[] = new double[36];     // coefficients for phi functions 
  double cm2[] = new double[36];     // specific to second order two dimensions 
  double cm3[] = new double[36];     // for three vertices, in order 
  int uniquev[] = new int[sz3];     // all vertices 
  int nvert;                        // number of unique vertices 
  int uniqueb[] = new int[sz3];     // bound vertices 
  int nbound;                       // number of boundary vertices 
  int uniquef[] = new int[sz3];     // free vertices, 0 to nfree-1 
  int nfree;                        // degrees of freedom  ntri-nbound 
  int nuniquev, nuniqueb, nuniquef; // vertices counts 
  double xs[] = new double[szint];  // integration coordinates, this triangle 
  double ys[] = new double[szint];  // integration coordinates, this triangle 
  double ws[] = new double[szint];  // integration coordinates, this triangle 
  double ph1[] = new double[szint]; // phi22(xs,ys) for this triangle cm1 
  double ph2[] = new double[szint]; // phi22(xs,ys) for this triangle cm2 
  double ph3[] = new double[szint]; // phi22(xs,ys) for this triangle cm3 
  double gs1[] = new double[szint]; // galk(cm1,xs,ys) 
  double gs2[] = new double[szint]; // galk(cm2,xs,ys) 
  double gs3[] = new double[szint]; // galk(cm3,xs,ys) 
  double fs[] = new double[szint];  // F(xs,ys)        
  
  fem_check22p_tri() // main constructor
  {
    double x, y;
    double x1, y1, x2, y2, x3, y3;
    int    i1,     i2,     i3;
    double a, b, sum, intg;
    double c[] = new double[np*np]; // max coefficients 
    double val, err, avgerr, maxerr;
    int stat, i, j, ii, k, kk, kdebug;
    int nv; // number of vertices for numerical integration 
    int nuvert; // number of unknown vertices after boundary applied 
    double T1[] = new double[6]; // x1,y1 x2,y2 x3,y3  this triangle

    System.out.println("fem_check22p_tri.c running ");
    System.out.println("uxx(x,y)+uyy(x,y)=f(x,y)");
    System.out.println("u(x,y)=x*x+2*x*y+3*y*y+4*x+5*y+6");
    System.out.println("f(x,y)=8");

    new triquad(); // get verification for tint

    stat = do_input();
    if(stat != 0) System.exit(1); // can not continue 


    // initialize kg, fg and ug  first cut, not using sparse matrix
    for(i=0; i<nvert; i++)
    {
      for(j=0; j<nvert; j++)
      {
        kg[i*nvert+j] = 0.0;
      }
    }

    for(i=0; i<nvert; i++)
    {
      fg[i] = 0.0; // zero except at boundary 
      ug[i] = 0.0;
    }
    // set boundary values (as used in fem_check22p_tri.c) 
    for(i=0; i<nbound; i++)
    {
      j=uniqueb[i];
      fg[j]=uana(xc[j],yc[j]);
      kg[j*nvert+j]=1.0;
    }

    // compute entries in stiffness matrix 
    System.out.println("compute global stiffness matrix ");
    for(k=0; k<ntri; k++) // generate contribution for each triangle 
    {
      System.out.println("computing local stiffness matrix for triangle "+k);
      i1=t1[k];
      i2=t2[k];
      i3=t3[k];
      x1=T1[0]=xc[i1];
      y1=T1[1]=yc[i1];
      x2=T1[2]=xc[i2];
      y2=T1[3]=yc[i2];
      x3=T1[4]=xc[i3];
      y3=T1[5]=yc[i3];
      kdebug=debug;  // (k<1?1:0); 
      if(kdebug>0)
      {
        System.out.println("nodes i1="+i1+", i2="+i2+", i3="+i3);
        System.out.println("coord ("+x1+","+y1+") ("+x2+","+y2+
                           ") ("+x3+","+y3+") ");
      }

      // determine coefficients for the three vertices 
      phi_tri_cm22(x1, y1, x2, y2, x3, y3, kdebug, cm1);
      phi_tri_cm22(x2, y2, x3, y3, x1, y1, kdebug, cm2);
      phi_tri_cm22(x3, y3, x1, y1, x2, y2, kdebug, cm3);
      // build integration points and values
      nv = tri_int1a(np, T1, kdebug);

      if(not_in(i1, uniqueb, nuniqueb)) // x1,y1 contributions if not boundary 
      {
        intg=tri_int2a(nv, gs1, ph1);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i1+", j="+i1);
        kg[i1*nvert+i1] += intg;
        intg=tri_int2a(nv, fs, ph1);
        if(kdebug>0) System.out.println("intg fg ="+intg+", at i="+i1);
        fg[i1] += intg;
        intg=tri_int2a(nv, gs2, ph1);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i2+", j="+i1);
        kg[i1*nvert+i2] += intg;
        // if(not_in(i2, uniqueb, nuniqueb)) kg[i2*nvert+i1] += intg; 
        intg=tri_int2a(nv, gs3, ph1);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i3+", j="+i1);
        kg[i1*nvert+i3] += intg;
        // if(not_in(i3, uniqueb, nuniqueb)) kg[i3*nvert+i1] += intg; 
      }

      if(not_in(i2, uniqueb, nuniqueb)) // x2,y2 contributions if not boundary 
      {
        intg=tri_int2a(nv, gs2, ph2);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i2+", j="+i2);
        kg[i2*nvert+i2] += intg;
        intg=tri_int2a(nv, fs, ph2);
        if(kdebug>0) System.out.println("intg fg ="+intg+", at i="+i2);
        fg[i2] += intg;
        intg=tri_int2a(nv, gs1, ph2);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i1+", j="+i2);
        kg[i2*nvert+i1] += intg;
        // if(not_in(i1, uniqueb, nuniqueb)) kg[i1*nvert+i2] += intg; 
        intg=tri_int2a(nv, gs3, ph2);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i3+", j="+i2);
        kg[i2*nvert+i3] += intg;
        // if(not_in(i3, uniqueb, nuniqueb)) kg[i3*nvert+i2] += intg; 
      }

      if(not_in(i3, uniqueb, nuniqueb)) // x3,y3 contributions if not boundary 
      {
        intg=tri_int2a(nv, gs3, ph3);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i3+", j="+i3);
        kg[i3*nvert+i3] += intg;
        intg=tri_int2a(nv, fs, ph3);
        if(kdebug>0) System.out.println("intg fg ="+intg+", at i="+i3);
        fg[i3] += intg;
        intg=tri_int2a(nv, gs1, ph3);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i1+", j="+i3);
        kg[i3*nvert+i1] += intg;
        // if(not_in(i1, uniqueb, nuniqueb)) kg[i1*nvert+i3] += intg; 
        intg=tri_int2a(nv, gs2, ph3);
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i2+", j="+i3);
        kg[i3*nvert+i2] += intg;
        // if(not_in(i2, uniqueb, nuniqueb)) kg[i2*nvert+i3] += intg; 
      }
      System.out.println("finished k="+k);
    } // end k loop making local stifness matrices, inserted into global  

    if(debug>1) for(i=0; i<ntri; i++)
    {
      System.out.println("tri "+i+" has vertices " +t1[i]+", "+t2[i]+", "+t3[i]);
    }

    if(debug>1) 
    {
      System.out.println("vertices, coordinates and analytic values ");
      for(i=0; i<nvert; i++)
      {
        Ua[i]=uana(xc[i],yc[i]);
        System.out.println("coordinate "+i+" at "+xc[i]+", "+yc[i]+
                           ", uana="+Ua[i]);
      }
    }

    if(debug>2)
    {
      System.out.println("k computed stiffness matrix, if debug ");
      for(i=0; i<nvert; i++)
        for(j=0; j<nvert; j++)
          System.out.println("K("+i+","+j+")="+kg[i*nvert+j]);
    }

    if(debug>2)
    {
      System.out.println("f computed forcing function and boundary, if debug ");
      for(i=0; i<nvert; i++)
        System.out.println("F("+i+")="+fg[i]+", ug["+i+"]="+ug[i]);
    }

    nuvert=nvert;

    new simeq(nvert, kg, fg, ug); 

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error ");
    for(i=0; i<nuvert; i++)
    {
      err = Math.abs(ug[i]-Ua[i]);
      if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    System.out.println("ug["+i+"]="+ug[i]+", Ua="+Ua[i]+", err="+(ug[i]-Ua[i]));
  }
  System.out.println("                 maxerr="+maxerr+", avgerr="+
                     (avgerr/(double)(nvert)));
  System.out.println("");

  } // end fem_check22p_tri constructor 


  // user provided, a specific differential equation 

  double F(double x, double y) // forcing function 
  {
    return 8.0;
  }

  double uana(double x, double y) // analytic solution and boundaries 
  {
    return x*x+2.0*x*y+3.0*y*y+4.0*x+5.0*y+6.0;
  }

  // this is the encoding of the differential equation 
  // L(u(x,y))=F(x,y), L(u(x,y)) becomes Galerkin L(phi(x,y)) 
  // u=phi22  ux=phi22x  uy=phi22y uxx=phi22xx uxy=phi22xy uyy=phi22yy 
  double galk(double c[], double x, double y)
  {                         // Galerkin k stiffness function for this problem 
    return phi22xx(c,x,y)+phi22yy(c,x,y);
  } // end galk used for integration, c for phi_i, phi_j separate 


  int tri_int1a(int np, double T[], int debug)
  {
    int nv, i;

    nv = triquad.tint(np, T, xs, ys, ws);
    if(debug>0) System.out.println("tri_int1a running with np="+np+", nv="+nv);
    for(i=0; i<nv; i++)
    {
      ph1[i]=phi22(cm1,xs[i],ys[i]);
      ph2[i]=phi22(cm2,xs[i],ys[i]);
      ph3[i]=phi22(cm3,xs[i],ys[i]);
      gs1[i]= galk(cm1,xs[i],ys[i]);
      gs2[i]= galk(cm2,xs[i],ys[i]);
      gs3[i]= galk(cm3,xs[i],ys[i]);
      fs[i]=F(xs[i],ys[i]);
      if(debug>2)
      {
	System.out.println("i="+i+", xs="+xs[i]+", ys="+ys[i]+", fs="+fs[i]);
	System.out.println("ph1="+ph1[i]+", ph2="+ph2[i]+", ph3="+ph3[i]);
	System.out.println("gs1="+gs1[i]+", gs2="+gs2[i]+", gs3="+gs3[i]);
      }
    }
    return nv;
  } // end tri_int1 

  double tri_int2a(int nv, double A[], double B[])
  {
    double sum;

    sum = 0.0;
    for(int i=0; i<nv; i++)
    {
      sum=sum+ws[i]*A[i]*B[i]; // area build into ws
      System.out.println("i="+i+", A="+A[i]+", B="+B[i]+", ws="+ws[i]);
    }
    return sum; // this is integral of dot product
  }

  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 

  //  phi22(c,x,y) is second order FEM phi function in 2 dimensions
  //  phi22x(c,x,y) is first derivative of phi21(c,x) wrt x
  //  phi22y(c,x,y) is first derivative of phi21(c,x) wrt y
  //  phi22xx(c,x,y) is second derivative of phi21(c,x) wrt x
  //  phi22yy(c,x,y) is second derivative of phi21(c,x) wrt y

  double phi22(double c[], double x, double y)
  {
    return c[0]+c[1]*x+c[2]*y+c[3]*x*x+c[4]*x*y+c[5]*y*y;
  }
  double phi22x(double c[], double x, double y)
  {
    return c[1]+2.0*c[3]*x+c[4]*y;
  }
  double phi22y(double c[], double x, double y)
  {
    return c[2]+c[4]*x+2.0*c[5]*y;
  }
  double phi22xx(double c[], double x, double y)
  {
    return 2.0*c[3];
  }
  double phi22xy(double c[], double x, double y)
  {
    return c[4];
  }
  double phi22yy(double c[], double x, double y)
  {
    return 2.0*c[5];
  }


  void phi_tri_cm22(double x1, double y1, double x2, double y2,
                    double x3, double y3, int kdebug, double cm[])
  {
    // specific to second order in two dimensions for triangles 
    // phi12(x,y)=cm[0]+cm[1]*x+cm[2]*y+cm[3]*x*x+cm[4]*x*y+cm[5]*y*y  
    double am[] = new double[36];
    double ym[] = new double[6];
    double sum;
    int i, j;

    am[0]=1.0;
    am[1]=x1;
    am[2]=y1;
    am[3]=x1*x1;
    am[4]=x1*y1;
    am[5]=y1*y1;

    am[6]=1.0;
    am[7]=x2;
    am[8]=y2;
    am[9]=x2*x2;
    am[10]=x2*y2;
    am[11]=y2*y2;

    am[12]=1.0;
    am[13]=x3;
    am[14]=y3;
    am[15]=x3*x3;
    am[16]=x3*y3;
    am[17]=y3*y3;

    double x12=(x1+x2)/2.0;
    double y12=(y1+y2)/2.0;
    am[18]=1.0;
    am[19]=x12;
    am[20]=y12;
    am[21]=x12*x12;
    am[22]=x12*y12;
    am[23]=y12*y12;

    double x13=(x1+x3)/2.0;
    double y13=(y1+y3)/2.0;
    am[24]=1.0;
    am[25]=x13;
    am[26]=y13;
    am[27]=x13*x13;
    am[28]=x13*y13;
    am[29]=y13*y13;

    double x23=(x2+x3)/2.0;
    double y23=(y2+y3)/2.0;
    am[30]=1.0;
    am[31]=x23;
    am[32]=y23;
    am[33]=x23*x23;
    am[34]=x23*y23;
    am[35]=y23*y23;

    ym[0]=1.0; // 1.0 for first, all others zero 
    ym[1]=0.0;
    ym[2]=0.0;
    ym[3]=0.0;
    ym[4]=0.0;
    ym[5]=0.0;
    if(kdebug>5)
    {
      System.out.println("am matrix  of  am cm = ym   solve for cm ");
      for(i=0; i<6; i++)
        for(j=0; j<6; j++)
          System.out.println("am["+i+","+j+"]= "+am[6*i+j]);
    }

    new simeq(6,am,ym,cm);
    if(kdebug>0)
    {
      System.out.println("cm matrix  of  am cm = ym   solve for cm ");
      for(i=0; i<6; i++) System.out.print(" "+cm[i]);
      System.out.println("");
    }
    // check solution 
    for(i=0; i<6; i++)
    {
      sum=0.0;
      for(j=0; j<6; j++)
        sum=sum+am[i*6+j]*cm[j];
      if(i==0 && Math.abs(sum-1.0)>1.0e-6)
          System.out.println("BAD cm i="+i+", sum="+sum);
      if(i!=0 && Math.abs(sum)>1.0e-6)
          System.out.println("BAD cm i="+i+", sum="+sum);
    }
    if(Math.abs(phi22(cm,x1,y1)-1.0)>1.0e-6)
      System.out.println("BAD cm x1="+x1+", y1="+y1);
    if(Math.abs(phi22(cm,x2,y2))>1.0e-6)
      System.out.println("BAD cm x2="+x2+", y2="+y2);
    if(Math.abs(phi22(cm,x3,y3))>1.0e-6)
      System.out.println("BAD cm x3="+x3+", y3="+y3);
    if(Math.abs(phi22(cm,x12,y12))>1.0e-6)
      System.out.println("BAD cm x12="+x12+", y12="+y12);
    if(Math.abs(phi22(cm,x13,y13))>1.0e-6)
      System.out.println("BAD cm x13="+x13+", y13="+y13);
    if(Math.abs(phi22(cm,x23,y23))>1.0e-6)
      System.out.println("BAD cm x23="+x23+", y23="+y23);
  } // end phi_tri_cm22 

  int do_input() // three files root .tri  root.bound  root.coord
  {
    String fname;
    int index, last, len;
    // 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)
      {
	 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);
         System.out.println("index="+index+", last="+last);
         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);
         System.out.println("index="+index+", last="+last);
         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;
         System.out.println("index="+index+", last="+last);
         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);

    System.out.println("subtracting minvert="+minvert+
                       " from all vertices, using base zero");
    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)
      {
	 System.out.println(input_line);
         // 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);
         System.out.println("index="+index+", last="+last);
         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;
         System.out.println("index="+index+", last="+last);
         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);

    System.out.println("subtracting minvert="+minvert+
                       " from all vertices, using base zero");
    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)
      {
	 System.out.println(input_line);
         // 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);
         System.out.println("index="+index+", last="+last);
         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;
         System.out.println("index="+index+", last="+last);
         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;

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