// fem_check21_tri.c  Galerkin FEM on triangular grid
// solve 3 ux(x,y)  + 2 uy(x,y)  + u(x,y) = f(x,y)
// boundary conditions computed using u(x,y)
// analytic solution may be given by u(x,y) = 1+ 2*x + 3*y
// f(x,y)=2.0*x+3.0*y+13.0
// input triangles and coordinates
// input geometry, root= 22_t  .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_check21_tri
{
  final int sz = 50;
  final int sz3 = 150;
  final int szsz = 2500;
  final int szint = 6000;

  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[6];     // coefficients for phi functions 
  double cm2[] = new double[6];     // specific to second order two dimensions 
  double cm3[] = new double[6];     // 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 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_check21_tri() // main constructor
  {
    double x, y, darea;
    double x1, y1, x2, y2, x3, y3;
    int    i1,     i2,     i3;
    double a, b, sum, intg;
    int np = 30;   // number of divisions for points to integrate 
    double c[] = new double[6]; // 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 

    System.out.println("fem_check21_tri.c running ");
    System.out.println("Given 3 ux + 2 uy + u = 2 x + 3 y + 13 ");
    System.out.println("Analytic solution u(x,y) =  1 + 2 x + 3 y ");

    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_check21_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=xc[i1];
      y1=yc[i1];
      x2=xc[i2];
      y2=yc[i2];
      x3=xc[i3];
      y3=yc[i3];
      kdebug=debug;  // (k<1?1:0); 
      darea=area(x1, y1, x2, y2, x3, y3);
      if(kdebug>0)
      {
        System.out.println("nodes i1="+i1+", i2="+i2+", i3="+i3+", area="+
                           darea);
        System.out.println("coord ("+x1+","+y1+") ("+x2+","+y2+
                           ") ("+x3+","+y3+") ");
      }
      // determine coefficients for the three vertices 
      phi_tri_cm12(x1, y1, x2, y2, x3, y3, kdebug, cm1);
      phi_tri_cm12(x2, y2, x3, y3, x1, y1, kdebug, cm2);
      phi_tri_cm12(x3, y3, x1, y1, x2, y2, kdebug, cm3);
      // build integration points and values
      nv = tri_int1(np, x1, y1, x2, y2, x3, y3, kdebug); 

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

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

      if(not_in(i3, uniqueb, nbound)) // x3,y3 contributions if not boundary 
      {
        intg=tri_int2(nv, gs3, ph3)*darea;
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i3+", j="+i3);
        kg[i3*nvert+i3] += intg;
        intg=tri_int2(nv, fs, ph3)*darea;
        if(kdebug>0) System.out.println("intg fg ="+intg+", at i="+i3);
        fg[i3] += intg;
        intg=tri_int2(nv, gs1, ph3)*darea;
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i1+", j="+i3);
        kg[i3*nvert+i1] += intg;
        // if(not_in(i1, uniqueb, nbound)) kg[i1*nvert+i3] += intg; 
        intg=tri_int2(nv, gs2, ph3)*darea;
        if(kdebug>0) System.out.println("integral="+intg+", at i="+i2+", j="+i3);
        kg[i3*nvert+i2] += intg;
        // if(not_in(i2, uniqueb, nbound)) 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); // bad news, indexing, solve large matrix 

    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_check21_tri constructor 


  // user provided, a specific differential equation 

  double F(double x, double y) // forcing function 
  {
    return 2.0*x+3.0*y+13.0;
  }

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

  // 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  3.0 * phi12x (c,x,y) +
            2.0 * phi12y (c,x,y) +
                  phi12  (c,x,y);
  } // end galk used for integration, c for phi_i, phi_j separate 



  int tri_int1(int np, double x1, double y1, double x2, double y2,
               double x3, double y3, int debug)
  {
    int nv, i;

    nv=phi_tri_pts(np, x1, y1, x2, y2, x3, y3, xs, ys);
    if(debug>0) System.out.println("tri_int1 running with np="+np+", nv="+nv);
    for(i=0; i<nv; i++)
    {
      ph1[i]=phi12(cm1,xs[i],ys[i]);
      ph2[i]=phi12(cm2,xs[i],ys[i]);
      ph3[i]=phi12(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]);
    }
    return nv;
  } // end tri_int1 

  double tri_int2(int nv, double A[], double B[])
  {
    double sum, weight;
    double w1 = 1.0;
    double w2 = 2.0;
    double w3 = 3.0;
    int i;

    sum=w1*(A[0]*B[0]+A[1]*B[1]+A[2]*B[2]);
    weight=3.0*w1;
    sum=sum+w2*(A[3]*B[3]+A[4]*B[4]+A[5]*B[5]);
    weight=weight+3.0*w2;
    for(i=6; i<nv; i++)
    {
      sum=sum+w3*A[i]*B[i];
      weight=weight+w3;
    }
    return sum/weight;
  }

  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;
  }

  double area(double P1x, double P1y,
              double P2x, double P2y,
              double P3x, double P3y)
  {
    double a;

    a = P1x*P2y - P2x*P1y + P2x*P3y - P3x*P2y + P3x*P1y - P1x*P3y;
    return Math.abs(a)/2.0;
  } // end area

  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 

  //  phi12(c,x,y) is first order FEM phi function in 2 dimensions
  //  phi12x(c,x,y) is first derivative of phi21(c,x) wrt x
  //  phi12y(c,x,y) is first derivative of phi21(c,x) wrt y

    double phi12(double c[], double x, double y)
  {
    return c[0]+c[1]*x+c[2]*y;
  }
    double phi12x(double c[], double x, double y)
  {
    return c[1];
  }
    double phi12y(double c[], double x, double y)
  {
    return c[2];
  }


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

    am[0]=1.0;
    am[1]=x1;
    am[2]=y1;
    am[3]=1.0;
    am[4]=x2;
    am[5]=y2;
    am[6]=1.0;
    am[7]=x3;
    am[8]=y3;

    ym[0]=1.0; // 1.0 for first, all others zero 
    ym[1]=0.0;
    ym[2]=0.0;
    new simeq(3,am,ym,cm);
    if(kdebug>0)
    {
      System.out.println("cm matrix  of  am cm = ym   solve for cm ");
      for(i=0; i<3; i++) System.out.print(" "+cm[i]);
      System.out.println("");
    }
    // check solution 
    for(i=0; i<3; i++)
    {
      sum=0.0;
      for(j=0; j<3; j++)
        sum=sum+am[i*3+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(phi12(cm,x1,y1)-1.0)>1.0e-6)
        System.out.println("BAD cm x1="+x1+", y1="+y1);
     if(Math.abs(phi12(cm,x2,y2))>1.0e-6)
        System.out.println("BAD cm x2="+x2+", y2="+y2);
     if(Math.abs(phi12(cm,x3,y3))>1.0e-6)
        System.out.println("BAD cm x3="+x3+", y3="+y3);
  } // end phi_tri_cm12 


  int phi_tri_pts(int np, double x1, double y1, double x2, double y2,
	          double x3, double y3, double xs[], double ys[])
  {
    double x1a, y1a, x2a, y2a, x3a, y3a, x12, y12, x32, y32, nf;
    int p, j, nv;

    nv=0;
    nf=(double)np;
    xs[nv]=x1a=((2*np-2)*x1+x2+x3)/(2.0*nf); // move in three points 
    ys[nv]=y1a=((2*np-2)*y1+y2+y3)/(2.0*nf);
    nv++;
    xs[nv]=x2a=((2*np-2)*x2+x3+x1)/(2.0*nf);
    ys[nv]=y2a=((2*np-2)*y2+y3+y1)/(2.0*nf);
    nv++;
    xs[nv]=x3a=((2*np-2)*x3+x1+x2)/(2.0*nf);
    ys[nv]=y3a=((2*np-2)*y3+y1+y2)/(2.0*nf);
    nv++;

    for(p=1; p<=np-1; p++)     // 3 edges, n points each 
    { 
      xs[nv]=(p/nf)*x1a+((np-p)/nf)*x2a;
      ys[nv]=(p/nf)*y1a+((np-p)/nf)*y2a;
      nv++;
      xs[nv]=(p/nf)*x2a+((np-p)/nf)*x3a;
      ys[nv]=(p/nf)*y2a+((np-p)/nf)*y3a;
      nv++;
      xs[nv]=(p/nf)*x3a+((np-p)/nf)*x1a;
      ys[nv]=(p/nf)*y3a+((np-p)/nf)*y1a;
      nv++;
    }

    for(j=1; j<=np-2; j++)
    {
      x12=(j/nf)*x2a+((np-j)/nf)*x1a;
      y12=(j/nf)*y2a+((np-j)/nf)*y1a;
      x32=((np-j)/nf)*x3a+(j/nf)*x2a;
      y32=((np-j)/nf)*y3a+(j/nf)*y2a;
      for(p=1; p<=np-1-j; p++)
      {    
        xs[nv]=(p/(nf-j))*x12+((np-p-j)/(nf-j))*x32;
        ys[nv]=(p/(nf-j))*y12+((np-p-j)/(nf-j))*y32;
        nv++;
      }
    }
    return nv;
  } // end phi_tri_pts 



  int do_input()
  {
    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_check21_tri();
  }
} // end fem_check21_tri class