// fem_check22e_la.java  Galerkin FEM
//                       second order partial differential equation
//                       rectangular grid, boundary value problem
//                       solve for points on the grid, inside boundary
//                       sixth order Gauss Legendre integration
//                       uses  laphi.java, simeq.java, gaulegf.java
//            uxx(x,y)+uyy(x,y)+r(x,y)*u(x,y)=f(x,y)
//            r(x,y)=sin(x)+cos(y)
//            f(x,y)=exp((x+2*y)/5)/5 +(sin(x)+cos(y))*exp((x+2*y)/5)
//            solution u(x,y)=exp((x+2*y)/5)

class fem_check22e_la
{
  // i, j, nx, ii, jj, ny, xg, yg needed by galk or galf, available at call 
  int nx;
  int ny;
  int nxy;
  double xg[];
  double yg[];
  int i, j, ii, jj;
  laphi L;
  boolean debug;

  fem_check22e_la(int Anx, int Any, int Anp, boolean Adebug )
  {
    nx = Anx;
    ny = Any;
    debug = Adebug;
    nxy = nx*ny;
    int npx = Anp; // Gauss Legendre integration order
    int npy = Anp;
    xg = new double[nx];
    yg = new double[ny];
    L = new laphi();

    double x, y, hx, hy;
    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 = 1.0; // problem parameters 
    double xmax = 3.0;
    double ymin = 1.0;
    double ymax = 2.0;

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


    System.out.println("fem_check22e_la.java running  ");
    System.out.println("Given:");
    System.out.println("uxx(x,y)+uyy(x,y)+r(x,y)*u(x,y)=f(x,y)");
    System.out.println("r(x,y)=sin(x)+cos(y)");
    System.out.println("f(x,y)=exp((x+2*y)/5)/5 +(sin(x)+cos(y))*exp((x+2*y)/5)");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax Boundaries  ");
    System.out.println("Analytic solution u(x,y)=exp((x+2*y)/5)");

    System.out.println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+", ymax="+ymax);
    System.out.println("nx="+nx+", ny="+ny);
    if(debug) 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 
      if(debug) System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    }  
    if(debug) System.out.println(" ");
    if(debug) 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 
      if(debug) System.out.println("ii="+ii+", Ua("+yg[ii]+")="+Ua[ii]);
    }  
    if(debug) 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);
        if(debug) System.out.println("solution at i="+i+",x="+x+", ii="+
                                     ii+",y="+y+" is"+Ua[i*ny+ii]);
      }
    }
    if(debug) 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);
      if(debug) 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);
      if(debug) 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);
      if(debug) 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);
      if(debug) 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("calling gauleg ymin="+ymin+", ymax="+ymax+
                       ", npy="+npy);
    new gaulegf(ymin, ymax, yy, wy, npy); // yy and wy set up for integrations 

    // 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]);
	      }
            }
            if(debug) 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]);
          }
        }
        if(debug) System.out.println("Legendre integration="+val+", f at i="+i+
                           ", ii="+ii);
        fg[i*ny+ii] = val;
      }
    }
    System.out.println("k computed stiffness matrix");
    System.out.println("f computed forcing vector");
    System.out.println("solving  k u = f  for  u");
    new simeq(kg, fg, ug);

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("u computed Galerkin, U 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;
        if(debug || (i!=0 && i!=nx-1 && ii!=0 && ii!=ny-1))
        System.out.println("u["+i+","+ii+"]="+ug[i*ny+ii]+", U="+Ua[i*ny+ii]+
                           ", err="+(ug[i*ny+ii]-Ua[i*ny+ii]));
      }
    }
    System.out.println("nx="+nx+", ny="+ny+", np="+npx);
    System.out.println("                   maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx*ny));
    System.out.println(" ");
  } // end fem_check22e_la constructor

    // PDE problem definition functions
    double F(double x, double y)
    {
      return Math.exp((x+2.0*y)/5.0)/5.0+(Math.sin(x)+Math.cos(y))*
             Math.exp((x+2.0*y)/5.0);
    }

    double r(double x, double y)
    {
      return Math.sin(x)+Math.cos(y);
    }

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

    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) +
	     r(x,y)*L.phi(x,j,nx-1,xg)*    L.phi(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 

  public static void main (String[] args)
  {
    // study grid refinement ny,ny with integration order, np
    new fem_check22e_la(5,4,3, true);  // nx,ny,np, debug
    new fem_check22e_la(5,4,4, false); // nx,ny,np, debug
    new fem_check22e_la(6,6,3, false); // nx,ny,np, debug
    new fem_check22e_la(6,6,4, false); // nx,ny,np, debug
    new fem_check22e_la(6,6,5, false); // nx,ny,np, debug
    new fem_check22e_la(8,8,5, false); // nx,ny,np, debug
    new fem_check22e_la(8,8,6, false); // nx,ny,np, debug
    new fem_check22e_la(8,8,7, false); // nx,ny,np, debug
    new fem_check22e_la(8,8,8, false); // nx,ny,np, debug
    new fem_check22e_la(10,10,8, false); // nx,ny,np, debug
    new fem_check22e_la(10,10,9, false); // nx,ny,np, debug
  }
} // end class fem_check22e_la