// fem_laplace_la.java  Galerkin FEM
// solve uxx+uyy=0
//       Analytic solution u(x,y)=exp(x)*sin(y)
//       inside 6 by 6 grid, Gauss Legendre integration

import java.io.*;

class fem_laplace_la
{
  // i, j, nx, ii, jj, ny, xg, yg needed by galk or galf, available at call 
  int i, j, ii, jj, nx, ny;
  double xg[]; // X grig
  double yg[]; // Y grid
  double hx, hy;
  laphi L = new laphi();
  double ug[]; // solution

  fem_laplace_la(double xmax, double ymax, int anx, int any, int npx, int npy)
  {
    nx = anx;
    ny = any;
    int nxy = nx*ny;
    xg = new double[nx];
    yg = new double[ny];
    double x, y;
    // npx, npy  Gauss Legendre integration order
    double xx[] = new double[npx+1];
    double wx[] = new double[npx+1]; // for Gauss-Legendre 
    double yy[] = new double[npy+1];
    double wy[] = new double[npy+1]; // for Gauss-Legendre 
    double val, err, avgerr, maxerr;
    int px, py;

    double xmin = 0.0; // problem parameters 
    double ymin = 0.0;

    double kg[][] = new double[nxy][nxy];  
    double fg[] = new double[nxy];
    ug = new double[nxy];
    double Ua[] = new double[nxy];


    System.out.println("fem_laplace_la.java running  ");
    System.out.println("Given:");
    System.out.println("uxx+uyy=0");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax Boundaries  ");
    System.out.println("Analytic solution u(x,y)=exp(x)*sim(y)");

    System.out.println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+
                       ", ymax="+ymax);
    System.out.println("nx="+nx+", ny="+ny);
    System.out.println("x grid and analytic solution at ymin  ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      Ua[i] = uana(xg[i],ymin); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    }  
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin ");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      Ua[ii] = uana(xmin, yg[ii]); // temp 
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+Ua[ii]);
    }  
    System.out.println(" ");
    // put solution in Ua vector 
    for(i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        Ua[i*ny+ii] = uana(x,y);
        // System.out.println("solution at i="+i+",x="+x+", ii="+ii+",y="+y+
        //                    " is"+Ua[i*ny+ii]);
      }
    }
    System.out.println(" ");

    // initialize k and f 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=0; ii<ny; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            kg[(i*ny+ii)][(j*ny+jj)] = 0.0;
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        fg[i*ny+ii] = 0.0;
      }
    }
    // set boundary conditions, four sides 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii=0;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = uana(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
      ii=ny-1;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = uana(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
    }
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i=0;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = uana(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
      i=nx-1;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = uana(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
    }

    System.out.println("calling gauleg xmin="+xmin+", xmax="+xmax+
                       ", npx="+npx);
    new gaulegf(xmin, xmax, xx, wx, npx); // xx and wx set up for integrations 
    System.out.println("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                       ", wx[1]="+wx[1]+", wx[2]="+wx[2]);
    System.out.println("calling gauleg ymin="+ymin+", ymax="+ymax+
                       ", npy="+npy);
    new gaulegf(ymin, ymax, yy, wy, npy); // yy and wy set up for integrations 
    System.out.println("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                       ", wy[1]="+wy[1]+", wy[2]="+wy[2]);
    i=1; ii=1;  j=1; jj=1;
    val = galk(xx[2],yy[2]);
    System.out.println("galk(xx[2],yy[2])="+val);
    val = galf(xx[2],yy[2]);
    System.out.println("galf(xx[2],yy[2])="+val);


    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix  ");
    for(i=1; i<nx-1; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=1; ii<ny-1; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            // compute integral for k(i,j)  
            val = 0.0;
            for(px=1; px<=npx; px++)
            {
              for(py=1; py<=npy; py++)
              {
                val = val + wx[px]*wy[py]*galk(xx[px],yy[py]);
	      }
            }
            // System.out.println("Legendre integration="+val+", at i="+i+
            //                    ", j="+j+", ii="+ii+", jj="+jj);
            kg[(i*ny+ii)][(j*ny+jj)] = val;
          } // end jj loop 
        } // end ii loop 
      } // end j loop 
    } // end i loop 

    // compute integral for f(i)  
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        val = 0.0;
        for(px=1; px<=npx; px++)
        {
          for(py=1; py<=npy; py++)
          {
            val = val + wx[px]*wy[py]*galf(xx[px],yy[py]);
          }
        }
        // System.out.println("Legendre integration="+val+", f at i="+i+
        //                    ", ii="+ii);
        fg[i*ny+ii] = val;
      }
    }
    System.out.println("k computed stiffness matrix, see above  ");
    System.out.println("f computed forcing function, see above  ");
    new simeq(kg, fg, ug);

    check_soln(); // against PDE definition, not knowing solution
    write_soln(); // to file x y soln

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        err = Math.abs(ug[i*ny+ii]-Ua[i*ny+ii]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        System.out.println("ug["+i+","+ii+"]="+ug[i*ny+ii]+", Ua="+Ua[i*ny+ii]+
                           ", err="+(ug[i*ny+ii]-Ua[i*ny+ii]));
      }
    }
    System.out.println("xmax="+xmax+", ymax="+ymax+", nx="+nx+
                       ", ny="+ny+", npx="+npx+", npy="+npy);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx*ny));
    System.out.println("\n");
  } // end fem_laplace_la constructor

  // PDE problem definition functions
  double F(double x, double y)
  {
    return 0.0;
  }

  double uana(double x, double y) // u 
  {
    return Math.exp(x)*Math.sin(y);
  }

  double galk(double x, double y)
  {                        // Galerkin k stiffness function for this problem 
    return ( L.phipp(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg) +
             L.phi(x,j,nx-1,xg)*    L.phipp(y,jj,ny-1,yg) )*
             L.phi(x,i,nx-1,xg)*    L.phi(y,ii,ny-1,yg);
  } // end galk used by integration 

  double galf(double x, double y) // Galerkin f function for this problem 
  {
    return F(x,y)*L.phi(x,i,nx-1,xg)*L.phi(y,ii,ny-1,yg);
  } // end galf used by integration 

  void write_soln()
  {
    try
    {
      FileWriter fout = new FileWriter("fem_laplace_la.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing fem_laplace_la.dat");

      for(int i=0; i<nx; i++)
      {
        for(int ii=0; ii<ny; ii++)
        {
          fileout.write(" "+xg[i]+" "+yg[ii]+" "+ug[i*ny+ii]+" \n");
        } // end y
        fileout.write(" \n");
      } // end x
      fileout.close();
      System.out.println("finished writing fem_laplace_la.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("fem_laplace_la.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("fem_laplace_la.java output-file error");
    } // end try
  } // end write_soln

  void check_soln_nu() // for all interior points
  {
    double err, sumerr, sumsqerr, rmserr, avgerr, maxerr;
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double eval_deriv;
    int cnt;

    System.out.println("check solution nu against PDE Uxx + Uyy = 0");
                       // -F(x,y) in general
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;

    for(int i=1; i<nx-1; i++)
    {
      // compute derivative coefficients in x direction
      new nuderiv(2,nx,i,xg,cxx);

      for(int ii=1; ii<ny-1; ii++)
      {
	// computer derivative coefficients in x direction
	new nuderiv(2,ny,ii,yg,cyy);
        eval_deriv = 0.0;
        
        for(int j=0; j<nx; j++)
        {
	  eval_deriv += cxx[j]*ug[j*ny+ii];
	} 
        for(int jj=0; jj<ny; jj++)
        {
	  eval_deriv += cyy[jj]*ug[i*ny+jj];
	}
        err = eval_deriv; // -F(xg[i],yg[ii]) not needed here
        sumerr += Math.abs(err);
        sumsqerr += err*err;
        maxerr = Math.max(Math.abs(err),maxerr);
        cnt++;  
      } // end y
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("nu maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln_nu

  void check_soln() // for all interior points
  {
    double err, sumerr, sumsqerr, rmserr, avgerr, maxerr;
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double eval_deriv;
    int cnt;

    System.out.println("check solution against PDE Uxx + Uyy = 0");
                       // -F(x,y) in general
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;
    nderiv D = new nderiv();

    for(int i=1; i<nx-1; i++)
    {
      // compute derivative coefficients in x direction
      D.rderiv(2,nx,i,hx,cxx);

      for(int ii=1; ii<ny-1; ii++)
      {
	// computer derivative coefficients in x direction
	D.rderiv(2,ny,ii,hy,cyy);
        eval_deriv = 0.0;
        
        for(int j=0; j<nx; j++)
        {
	  eval_deriv += cxx[j]*ug[j*ny+ii];
	} 
        for(int jj=0; jj<ny; jj++)
        {
	  eval_deriv += cyy[jj]*ug[i*ny+jj];
	}
        err = eval_deriv; // -F(xg[i],yg[ii]) not needed here
        sumerr += Math.abs(err);
        sumsqerr += err*err;
        maxerr = Math.max(Math.abs(err),maxerr);
        cnt++;  
      } // end y
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln 

  public static void main (String[] args)
  {
    new fem_laplace_la(1.0, 1.0, 4, 4, 3, 3); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 4, 4, 4, 4); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 4, 4, 6, 6); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 6, 6, 4, 4); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 6, 6, 6, 6); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 6, 6, 8, 8); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 8, 8, 6, 6); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.0, 1.0, 8, 8, 8, 8); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(1.5, 1.5, 8, 8, 8, 8); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(3.0, 3.0, 8, 8, 8, 8); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(3.0, 3.0, 10, 10, 8, 8); // xmax, ymax, nx, ny, npx, npy
    new fem_laplace_la(3.0, 6.3, 10, 10, 8, 8); // xmax, ymax, nx, ny, npx, npy
    try
    {
      Runtime runtime = Runtime.getRuntime();
      Process proc = runtime.exec("gnuplot -persist fem_la.plot");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("fem_laplace_la.java runtime error");
    } // end try
  }
} // end class fem_laplace_la