// pde22js_la  Galerkin FEM  
// solve ccxx*uxx+cxy*uxy+ccyy*uyy+ccx*ux+ccy*uy+cc*u=F(x,y)
//       given f(x,y) =
//       Analytic solution u(x,y) = boundary and check
//       inside nx by ny grid, Gauss Legendre integration

class pde22js_la
{
  laphi L = new laphi();
pde22js_la()
{
  double start, now;
  boolean debug = false;
  double ccxx = 1.0; // PDE definition
  double ccxy = 2.0;
  double ccyy = 3.0;
  double ccx  = 4.0;
  double ccy  = 5.0;
  double cc   = 6.0;
  final int nx = 5;
  final int ny = 6;
  int nxy = nx*ny;
  double x, y, hx, hy, tmp;
  int npx = 10; // Gauss Legendre integration order
  int npy = 11;
  double xmin = 0.0; // problem parameters 
  double xmax = 1.0;
  double ymin = 0.1;
  double ymax = 1.1;
  double val, err, avgerr, maxerr;
  double xg[] = new double[nx]; 
  double yg[] = new double[ny];
  double xx[] = new double[npx+1]; // for Gauss-Legendre 
  double wx[] = new double[npx+1]; 
  double yy[] = new double[npy+1];
  double wy[] = new double[npy+1];
  double kg[][] = new double[nxy][nxy];  
  double fg[] = new double[nxy];
  double ug[] = new double[nxy];
  double Ua[] = new double[nxy];

  System.out.print("pde22js_la.java running \n");
  System.out.print("solve ccxx*uxx+cxy*uxy+ccyy*uyy+ccx*ux+ccy*uy+cc*u=F(x,y)\n");
  System.out.print("given  F(x,y) \n");
  System.out.print("  xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  System.out.print("  boundary and possibly check uana(x,y) \n");
  System.out.print("ccxx="+ccxx+", ccxy="+ccxy+", ccyy="+ccyy+" \n");
  System.out.print("ccx="+ccx+", ccy="+ccy+", cc="+cc+" \n");
  System.out.print("xmin="+xmin+", xmax="+xmax+", nx="+nx+" \n");
  System.out.print("ymin="+ymin+", ymax="+ymax+", ny="+ny+" \n");
  System.out.print("npx="+npx+", npy="+npy+" \n");
  start = System.currentTimeMillis()/1000.0;

  hx = (xmax-xmin)/(nx-1);
  for(int i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
  }  
  hy = (ymax-ymin)/(ny-1);
  for(int ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
  }  
  System.out.print("xg and yg initialized \n");
  // put solution in Ua vector 
  for(int i=0; i<nx; i++)
  {
    x = xg[i];
    for(int ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      Ua[S(i,ii,ny)] = uana(x,y);
      if(debug)
        System.out.print("solution at i="+i+",x="+x+", ii="+ii+",y="+y+
	   	         " is "+Ua[S(i,ii,ny)]+"\n");
    }
  }
  if(debug) System.out.print("xg and yg initialized \n");

  // initialize k and f 
  for(int i=0; i<nx; i++)
  {
    for(int j=0; j<nx; j++)
    {
      for(int ii=0; ii<ny; ii++)
      {
        for(int jj=0; jj<ny; jj++)
        {
	  kg[S(i,ii,ny)][S(j,jj,ny)] = 0.0;
	}
      }
    }
  }
  for(int i=0; i<nx; i++)
  {
    for(int ii=0; ii<ny; ii++)
    {
      fg[S(i,ii,ny)] = 0.0;
    }
  }
  // set boundary conditions, four sides 
  for(int i=0; i<nx; i++)
  {
    x = xg[i];
    for(int ii=0; ii<ny; ii=ii+ny-1) // ii=0, ii=ny-1
    {
      y = yg[ii];
      kg[S(i,ii,ny)][S(i,ii,ny)] = 1.0;
      fg[S(i,ii,ny)] = uana(x,y);
      System.out.print("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
 		       " is "+fg[S(i,ii,ny)]+"\n");
    }
  }
  for(int ii=0; ii<ny; ii++)
  {
    y = yg[ii];
    for(int i=0; i<nx; i=i+nx-1) // i=0, i=nx-1
    {
      x = xg[i];
      kg[S(i,ii,ny)][S(i,ii,ny)] = 1.0;
      fg[S(i,ii,ny)] = uana(x,y);
      System.out.print("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
                       " is "+fg[S(i,ii,ny)]+"\n");
    }
  }
  // initialize xx, wx, yy, wy
  for(int k=0; k<=npx; k++)
  {
    xx[k] = 0.0;
    wx[k] = 0.0;
  }
  for(int k=0; k<=npy; k++)
  {
    yy[k] = 0.0;
    wy[k] = 0.0;
  }
  System.out.print("calling gauleg xmin="+xmin+", xmax="+xmax+
                     ", npx="+npx+"\n");
  new gaulegf(xmin, xmax, xx, wx, npx); // xx and wx set up for integrations 
  System.out.print("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                     ", wx[1]="+wx[1]+", wx[2]="+wx[2]+"\n");
  System.out.print("calling gauleg ymin="+ymin+", ymax="+ymax+
                 ", npy="+npy+"\n");
  new gaulegf(ymin, ymax, yy, wy, npy); // yy and wy set up for integrations 
  System.out.print("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                 ", wy[1]="+wy[1]+", wy[2]="+wy[2]+"\n");
  System.out.print("calling F( xmax="+xmax+", ymax="+ymax+"\n");
  val = F(xmax, ymax);
  System.out.print("F(xmax,ymax)="+val+"\n");
  System.out.print("calling phi(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phi(xmin, 0, nx-1, xg);
  System.out.print("phi(xmin,j=0,nx-i,xg)="+val+"\n");
  System.out.print("calling phip(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phip(xmin, 0, nx-1, xg);
  System.out.print("phip(xmin,j=0,nx-i,xg)="+val+"\n");
  System.out.print("calling phipp(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phipp(xmin, 0, nx-1, xg);
  System.out.print("phipp(xmin,j=0,nx-i,xg)="+val+"\n");

  val = galk(xx[2],yy[2],1,1,1,1,ccxx,ccxy,ccyy,ccx,ccy,cc,
             nx,ny,xg,yg);
  System.out.print("galk(xx[2],yy[2])="+val+"\n");
  val = galf(xx[2],yy[2],1,1,nx,ny,xg,yg);
  System.out.print("galf(xx[2],yy[2])="+val+"\n");
  // compute entries in stiffness matrix 
  System.out.print("compute stiffness matrix  \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int j=0; j<nx; j++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int jj=0; jj<ny; jj++)
        {
          // compute integral for k(i,j)  
          val = 0.0;
          for(int px=1; px<=npx; px++)
          {
            for(int py=1; py<=npy; py++)
            {
	      val = val + wx[px]*wy[py]*
		  galk(xx[px],yy[py],i,ii,j,jj,
                       ccxx,ccxy,ccyy,ccx,ccy,cc,
                       nx,ny,xg,yg);
            } // end py
          } // end px
          if(debug)
            System.out.print("Legendre integration="+val+", at i="+i+
                             ", j="+j+", ii="+ii+", jj="+jj+"\n");
          kg[S(i,ii,ny)][S(j,jj,ny)] = val;
        } // end jj 
      } // end ii 
    } // end j 
  } // end i 

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

  avgerr = 0.0;
  maxerr = 0.0;
  System.out.print("ug computed Galerkin, Ua analytic, error \n");
  for(int i=0; i<nx; i++)
  {
    for(int ii=0; ii<ny; ii++)
    {
      err = Math.abs(ug[S(i,ii,ny)]-Ua[S(i,ii,ny)]);
      if(err>maxerr) maxerr = err;
      avgerr = avgerr + err;
      System.out.print("ug["+i+","+ii+"]="+ug[S(i,ii,ny)]+", Ua="+
                       Ua[S(i,ii,ny)]+", err="+
                       (ug[S(i,ii,ny)]-Ua[S(i,ii,ny)])+"\n");
    }
  }
  System.out.print("                  maxerr="+maxerr+", avgerr="+
                 (avgerr/(nx*ny))+"\n");

  now = System.currentTimeMillis()/1000.0;
  System.out.println("total time "+(now-start)+" seconds");
  System.out.print("pde22js_la.java finished\n");
} // end pde22js_la

int S(int i, int j, int ny) // zero based index includes boundary
{                           // this is row of k[][]
    return i*ny+j; // index into us and ut, RHS
} // end S

// from uana22_fem.js  problem definition
double F(double x, double y) // RHS  
{
    return 6.0*x*x*x+6.0*y*y*y+12.0*x*x+15.0*y*y+6.0*x*y+11.0*x+22.0*y+8.0;
} // end f

//       Analytic solution
double uana(double x, double y)
{
  return x*x*x + y*y*y + x*y + 1.0;
} // end u

// needs laphi.java available
// for finite element method (not needed if coefficients available)
double galk(double x, double y, int i, int ii, int j, int jj,
            double ccxx, double ccxy, double ccyy, double ccx,
            double ccy,  double cc,
            int nx, int ny, double xg[], double yg[])
{           // Galerkin k stiffness function for this problem 
  return (ccxx * L.phipp(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg) +
          ccxy *  L.phip(x,j,nx-1,xg)* L.phip(y,jj,ny-1,yg) +
          ccyy *   L.phi(x,j,nx-1,xg)*L.phipp(y,jj,ny-1,yg) +
          ccx  *  L.phip(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg) +
          ccy  *   L.phi(x,j,nx-1,xg)* L.phip(y,jj,ny-1,yg) +
          cc   *   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, int i, int ii, int nx, int ny,
            double xg[], double yg[])
{           // 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) // "main" required
{
  new pde22js_la(); // construct and execute
} // end main

} // end class pde22js_la