/* test_nuderiv2dgc.c  corner test
 *
 *      B           2    five(5) known boundary points B
 *        P           5  five(5) unknown points P, DOF = 5
 *      B P         1 4
 *        P P P       3 7 8
 *      B   B   B   0   6   9
 */
/* solve Uxx(x,y) + 2*Uyy(x,y) = 14 */
/* boundary and solution U(x,y) = x^2+2xy+3y^2+4x+5y+6 */

#include "nuderiv2dg.h"
#include <stdio.h>
#include <stdlib.h>

double f(double x, double y) /* RHS */
{
  return 14.0;
}

double ub(double x, double y) /* solution for boundary and check */
{
  return x*x + 2.0*x*y + 3.0*y*y + 4.0*x + 5.0*y + 6.0;
}

int main(int argc, char *argv[])
{
  double xp[10] = {0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.5, 2.0};
  double yp[10] = {0.0, 1.0, 2.0, 0.5, 1.0, 1.5, 0.0, 0.5, 0.5, 0.0};
  int bound[10] = {  0,   0,   0,   1,   1,   1,   0,   1,   1,   0};
  double U[10];
  double cxx[10];
  double cyy[10];
  double uxx, uyy, sum, err;
  int ipx[10], ipy[10]; /* index of point used */
  int npoint = 10;
  int nupx, nupy; /* number of points used in cxx, cyy, etc. */
  int point;
  int ndof;
  int i, j;

  printf("test_nuderiv2dgc.c on corner\n");
  printf("solve Uxx(x,y) + 2*Uyy(x,y) = 14 \n");
  printf("boundary and solution U(x,y) = x^2+2xy+3y^2+4x+5y+6 \n");
  printf("B           2       B = boundary \n");
  printf("  P           5     P = DOF \n");
  printf("B P         1 4       \n");
  printf("  P P P       3 7 8   \n");
  printf("B   B   B   0   6   9 \n");

  for(i=0; i<npoint; i++)
  {
    U[i] = ub(xp[i],yp[i]);
    printf("i=%2d, b=%2d, xp=%f, yp=%f, U=%f\n",
	   i, bound[i], xp[i], yp[i], U[i]);
  }

  ndof = 0;
  for(i=0; i<npoint; i++)
  {
    if(bound[i]) /* not a boundary */
    {
      printf("\nbuild derivative for Point=%d \n", i);
      ndof++;
      nupx = nu2dxx(2, npoint, i, xp, yp, cxx, ipx);
      printf("nu2dxx at %d, nupx=%d \n", i, nupx);
      uxx = 0.0;
      for(j=0; j<nupx; j++)
      {
        uxx += cxx[j]*U[ipx[j]];
        printf("     j=%2d, ipx[j]=%2d, cxx[j]=%f \n", j, ipx[j], cxx[j]);
      }
      nupy = nu2dyy(2, npoint, i, xp, yp, cyy, ipy);
      printf("nu2dyy at %d, nupy=%d \n", i, nupy);
      uyy = 0.0;
      for(j=0; j<nupy; j++)
      {
        uyy += cyy[j]*U[ipy[j]];
        printf("     j=%2d, ipy[j]=%2d, cyy[j]=%f \n", j, ipy[j], cyy[j]);
      }
      sum = uxx+2.0*uyy;
      printf("at %2d sum=%f, f[i]=%f, err=%e\n",
             i, sum, f(xp[i],yp[i]), f(xp[i],yp[i])-sum); 
    } /* end not bound */
  } /* end i=point */
  printf("end test_nuderiv2dgc.c ndof=%d\n", ndof);
  return 0;
} /* end test_nuderiv2dgc.c */