# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <ctime>
# include <cmath>

using namespace std;

int main ( );
double exact ( double x, double y );
double r8_abs ( double x );
double *r8ge_fs ( int n, double a[], double b[] );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main routine of the finite element program FEM2D_PROJECT_FUNCTION.
//
//  Discussion:
//
//    This program seeks the continuous piecewise linear function U(X,Y) which
//    minimizes the least squares approximation error to a given function W(X,Y),
//    while satisfying given boundary conditions.
//
//    For basis functions V(x,y), we seek U(x,y) such that
//
//      ( U(x,y) - W(x,y), V(x,y) ) = 0
//
//    in a rectangular region in the plane.
//
//    At nodes on the boundary, exact conditions are imposed:
//
//      U(x,y) = W(x,y)
//
//    The code uses continuous piecewise linear basis functions on
//    triangles determined by a uniform grid of NX by NY points.
//
//    In this version of the program, the function W is defined as:
//
//      W(x,y)  = sin ( pi * x ) * sin ( pi * y ) + x
//
//  THINGS YOU CAN EASILY CHANGE:
//
//    1) Change NX or NY, the number of nodes in the X and Y directions.
//    2) Change XL, XR, YB, YT, the left, right, bottom and top limits of the rectangle.
//    3) Change the exact solution in the EXACT routine.
//
//  HARDER TO CHANGE:
//
//    4) Change from "linear" to "quadratic" triangles;
//    5) Change the region from a rectangle to a general triangulated region;
//    6) Store the matrix as a sparse matrix so you can solve bigger systems.
//    7) Handle Neumann boundary conditions.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    02 June 2009
//
//  Author:
//
//    John Burkardt
//
{
  int nx = 17;
  int ny = 17;

  double *a;
  double area;
  double *b;
  double *c;
  int e;
  int *element_node;
  int element_num;
  int i;
  int i1;
  int i2;
  int i3;
  int info;
  int j;
  int j2;
  int k;
  int node_num;
  int nq1;
  int nq2;
  int nti1;
  int nti2;
  int nti3;
  int ntj1;
  int ntj2;
  int ntj3;
  double pi = 3.141592653589793;
  int q;
  int q1;
  double q16[6] = { 0.0, 0.5, 0.5, 4.0/6.0, 1.0/6.0, 1.0/6.0 };
  int q2;
  double q26[6] = { 0.5, 0.0, 0.5, 1.0/6.0, 4.0/6.0, 1.0/6.0 };
  double q36[6] = { 0.5, 0.5, 0.0, 1.0/6.0, 1.0/6.0, 4.0/6.0 };
  double qi;
  double qj;
  int ti1;
  int ti2;
  int ti3;
  int tj1;
  int tj2;
  int tj3;
  double rhs;
  double u;
  double u_norm;
  double uw_norm;
  double w;
  double w_norm;
  double wq;
  double wq6[6] = { 1.0/30.0, 1.0/30.0, 1.0/30.0, 9.0/30.0, 9.0/30.0, 9.0/30.0 };
  double *x;
  double xl = 0.0;
  double xq;
  double xr = 1.0;
  double *y;
  double yb = 0.0;
  double yq;
  double yt = 1.0;

  cout << "\n";
  cout << "FEM2D_PROJECT_FUNCTION\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "  Seek U(x,y), the solution of the least squares equation:\n";
  cout << "\n";
  cout << "  Minimize L2 norm of U(x,y) - W(x,y), for W(x,y) given,\n";
  cout << "  with U(x,y) a piecewise linear function in the interior,\n";
  cout << "  and matching W(x,y) on the boundary.\n";
  cout << "\n";
  cout << "  Reformulate this in terms of a finite element problem:\n";
  cout << "\n";
  cout << "  ( U(x,y) - W(x,y), V(x,y) ) = 0 inside the region,\n";
  cout << "    U(x,y)                    = W(x,y) on the boundary\n";
  cout << "\n";
  cout << "  The region is a rectangle, defined by:\n";
  cout << "\n";
  cout << "  " << xl << " = XL<= X <= XR = " << xr << "\n";
  cout << "  " << yb << " = YB<= Y <= YT = " << yt << "\n";
  cout << "\n";
  cout << "  The finite element method is used, with piecewise\n";
  cout << "  linear basis functions on 3 node triangular\n";
  cout << "  elements.\n";
  cout << "\n";
  cout << "  The corner nodes of the triangles are generated by an\n";
  cout << "  underlying grid whose dimensions are\n";
  cout << "\n";
  cout << "  NX =                       " << nx << "\n";
  cout << "  NY =                       " << ny << "\n";
//
//  NODE COORDINATES
//
//  Numbering of nodes is suggested by the following 5x10 example:
//
//    J=5 | K=41  K=42 ... K=50
//    ... |
//    J=2 | K=11  K=12 ... K=20
//    J=1 | K= 1  K= 2     K=10
//        +--------------------
//          I= 1  I= 2 ... I=10
//
  node_num = nx * ny;

  cout << "  Number of nodes =          " << node_num << "\n";

  x = new double[node_num];
  y = new double[node_num];

  k = 0;
  for ( j = 1; j <= ny; j++ )
  {
    for ( i = 1; i <= nx; i++ )
    {
      x[k] = ( ( double ) ( nx - i     ) * xl   
             + ( double ) (      i - 1 ) * xr ) 
             / ( double ) ( nx     - 1 );

      y[k] = ( ( double ) ( ny - j     ) * yb   
             + ( double ) (      j - 1 ) * yt ) 
             / ( double ) ( ny     - 1 );
      k = k + 1;
    }
  }
//
//  ELEMENT array
//
//  Organize the nodes into a grid of 3-node triangles.
//  Here is part of the diagram for a 5x10 example:
//
//    |  \ |  \ |  \ |
//    |   \|   \|   \|
//   21---22---23---24--
//    |\ 8 |\10 |\12 |
//    | \  | \  | \  |
//    |  \ |  \ |  \ |  \ |
//    |  7\|  9\| 11\|   \|
//   11---12---13---14---15---16---17---18---19---20
//    |\ 2 |\ 4 |\ 6 |\  8|                   |\ 18|
//    | \  | \  | \  | \  |                   | \  |
//    |  \ |  \ |  \ |  \ |      ...          |  \ |
//    |  1\|  3\|  5\| 7 \|                   |17 \|
//    1----2----3----4----5----6----7----8----9---10
//
  element_num = 2 * ( nx - 1 ) * ( ny - 1 );

  cout << "  Number of elements =       " << element_num << "\n";

  element_node = new int[3*element_num];

  k = 0;

  for ( j = 1; j <= ny - 1; j++ )
  {
    for ( i = 1; i <= nx - 1; i++ )
    {
      element_node[0+k*3] = i     + ( j - 1 ) * nx - 1;
      element_node[1+k*3] = i + 1 + ( j - 1 ) * nx - 1;
      element_node[2+k*3] = i     +   j       * nx - 1;
      k = k + 1;

      element_node[0+k*3] = i + 1 +   j       * nx - 1;
      element_node[1+k*3] = i     +   j       * nx - 1;
      element_node[2+k*3] = i + 1 + ( j - 1 ) * nx - 1;
      k = k + 1;
    }
  }
//
//  Assemble the coefficient matrix A and the right-hand side B of the
//  finite element equations, ignoring boundary conditions.
//
  a = new double[node_num * node_num ];
  b = new double[node_num];

  for ( i = 0; i < node_num; i++ )
  {
    b[i] = 0.0;
  }
  for ( j = 0; j < node_num; j++ )
  {
    for ( i = 0; i < node_num; i++ )
    {
      a[i+j*node_num] = 0.0;
    }
  }

  for ( e = 0; e < element_num; e++ )
  {
    i1 = element_node[0+e*3];
    i2 = element_node[1+e*3];
    i3 = element_node[2+e*3];
    area = 0.5 * 
      ( x[i1] * ( y[i2] - y[i3] ) 
      + x[i2] * ( y[i3] - y[i1] ) 
      + x[i3] * ( y[i1] - y[i2] ) );
//
//  Consider each quadrature point.
//  Here, we use the midside nodes as quadrature points.
//
    for ( q1 = 0; q1 < 3; q1++ )
    {
      q2 = ( q1 + 1 ) % 3;

      nq1 = element_node[q1+e*3];
      nq2 = element_node[q2+e*3];

      xq = 0.5 * ( x[nq1] + x[nq2] );
      yq = 0.5 * ( y[nq1] + y[nq2] );
      wq = 1.0 / 3.0;
//
//  Consider each test function in the element.
//
      for ( ti1 = 0; ti1 < 3; ti1++ )
      {
        ti2 = ( ti1 + 1 ) % 3;
        ti3 = ( ti1 + 2 ) % 3;

        nti1 = element_node[ti1+e*3];
        nti2 = element_node[ti2+e*3];
        nti3 = element_node[ti3+e*3];

        qi = 0.5 * ( 
            ( x[nti3] - x[nti2] ) * ( yq - y[nti2] ) 
          - ( y[nti3] - y[nti2] ) * ( xq - x[nti2] ) ) / area;

        w = exact ( xq, yq );

        b[nti1] = b[nti1] + area * wq * ( w * qi );
//
//  Consider each basis function in the element.
//
        for ( tj1 = 0; tj1 < 3; tj1++ )
        {
          tj2 = ( tj1 + 1 ) % 3;
          tj3 = ( tj1 + 2 ) % 3;

          ntj1 = element_node[tj1+e*3];
          ntj2 = element_node[tj2+e*3];
          ntj3 = element_node[tj3+e*3];

          qj = 0.5 * ( 
              ( x[ntj3] - x[ntj2] ) * ( yq - y[ntj2] ) 
            - ( y[ntj3] - y[ntj2] ) * ( xq - x[ntj2] ) ) / area;

          a[nti1+ntj1*node_num] = a[nti1+ntj1*node_num] + area * wq * ( qi * qj );
        }
      }
    }
  }
//
//  BOUNDARY CONDITIONS
//
//  If the K-th variable is at a boundary node, replace the K-th finite
//  element equation by a boundary condition that sets the variable to U(K).
//
  k = 0;
  for ( j = 1; j <= ny; j++ )
  {
    for ( i = 1; i <= nx; i++ )
    {
      if ( i == 1 || i == nx || j == 1 || j == ny )
      {
        w = exact ( x[k], y[k] );
        for ( j2 = 0; j2 < node_num; j2++ )
        {
          a[k+j2*node_num] = 0.0;
        }
        a[k+k*node_num] = 1.0;
        b[k] = w;
      }
      k = k + 1;
    }
  }
//  SOLVE the linear system A * C = B.
//
//  The solution is returned in C.
//
  c = r8ge_fs ( node_num, a, b );
//
//  COMPARE U and W at the grid points only.
//  Unless W is itself a finite element function, we can't expect these values to
//  be equal.  But we aren't trying to match the pointwise data.
//
  cout << "\n";
  cout << "     K     I     J          X           Y        U(x,y)          W(x,y)\n";
  cout << "\n";

  k = 0;

  for ( j = 1; j <= ny; j++ )
  {
    for ( i = 1; i <= nx; i++ )
    {
      w = exact ( x[k], y[k]  );

      cout << "  " << setw(4) << k
           << "  " << setw(4) << i
           << "  " << setw(4) << j
           << "  " << setw(10) << x[k]
           << "  " << setw(10) << y[k]
           << "  " << setw(14) << c[k]
           << "  " << setw(14) << w
           << "  " << setw(14) << r8_abs ( w - c[k] ) << "\n";

      k = k + 1;
    }
    cout << "\n";
  }
//
//  COMPUTE the L2 norm of U - W over the region.  
//  This is the quanity we want to be small.
//
//  Here, use a 6 point quadrature rule.
//
  u_norm = 0.0;
  w_norm = 0.0;
  uw_norm = 0.0;

  for ( e = 0; e < element_num; e++ )
  {
    i1 = element_node[0+e*3];
    i2 = element_node[1+e*3];
    i3 = element_node[2+e*3];

    area = 0.5 * 
      ( x[i1] * ( y[i2] - y[i3] ) 
      + x[i2] * ( y[i3] - y[i1] ) 
      + x[i3] * ( y[i1] - y[i2] ) );
//
//  Consider each quadrature point.
//
    for ( q = 0; q < 6; q++ )
    {
      xq = q16[q] * x[i1] + q26[q] * x[i2] + q36[q] * x[i3];
      yq = q16[q] * y[i1] + q26[q] * y[i2] + q36[q] * y[i3];
//
//  Inside element E, W is the sum of nodal values B times the basis functions.
//
      u = 0.0;

      for ( ti1 = 0; ti1 < 3; ti1++ )
      {
        ti2 = ( ti1 + 1 ) % 3;
        ti3 = ( ti1 + 2 ) % 3;

        nti1 = element_node[ti1+e*3];
        nti2 = element_node[ti2+e*3];
        nti3 = element_node[ti3+e*3];

        qi = 0.5 * ( 
            ( x[nti3] - x[nti2] ) * ( yq - y[nti2] ) 
          - ( y[nti3] - y[nti2] ) * ( xq - x[nti2] ) ) / area;

        u = u + c[i1] * qi;
      }
      w = exact ( xq, yq );
//
//  Add the value of ( U - W )^2 to the quadrature sum.
//
      u_norm  = u_norm  + area * wq6[q] * u * u;
      w_norm  = w_norm  + area * wq6[q] * w * w;
      uw_norm = uw_norm + area * wq6[q] * ( u - w ) * ( u - w );
    }
  }

  u_norm  = sqrt ( u_norm );
  w_norm  = sqrt ( w_norm );
  uw_norm = sqrt ( uw_norm );

  cout << "\n";
  cout << "  ||U||   = " << u_norm << "\n";
  cout << "  ||W||   = " << w_norm << "\n";
  cout << "  ||U-W|| = " << uw_norm << "\n";
//
//  TERMINATE the program.
//
//  Deallocate memory.
//
  delete [] a;
  delete [] b;
  delete [] c;
  delete [] element_node;
  delete [] x;
  delete [] y;

  cout << "\n";
  cout << "FEM2D_PROJECT_FUNCTION:\n";
  cout << "  Normal end of execution.\n";

  return 0;
}
//****************************************************************************80

double exact ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    EXACT calculates the exact solution and its first derivatives.
//
//  Discussion:
//
//    The function specified here depends on the problem being
//    solved.  The user must be sure to change both EXACT and RHS
//    or the program will have inconsistent data.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    02 December 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the coordinates of a point
//    in the region, at which the exact solution is to be evaluated.
//
//    Output, double EXACT, the value of the exact solution.
//
{
  double pi = 3.141592653589793;
  double u;

  u = sin ( pi * x ) * sin ( pi * y ) + x;
       
  return u;
}
//****************************************************************************80

double r8_abs ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    R8_ABS returns the absolute value of an R8.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    14 November 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the quantity whose absolute value is desired.
//
//    Output, double R8_ABS, the absolute value of X.
//
{
  double value;

  if ( 0.0 <= x )
  {
    value = x;
  } 
  else
  {
    value = - x;
  }
  return value;
}
//****************************************************************************80

double *r8ge_fs ( int n, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8GE_FS factors and solves a R8GE system.
//
//  Discussion:
//
//    The R8GE storage format is used for a "general" M by N matrix.  
//    A physical storage space is made for each logical entry.  The two 
//    dimensional logical array is mapped to a vector, in which storage is 
//    by columns.
//
//    R8GE_FS does not save the LU factors of the matrix, and hence cannot
//    be used to efficiently solve multiple linear systems, or even to
//    factor A at one time, and solve a single linear system at a later time.
//
//    R8GE_FS uses partial pivoting, but no pivot vector is required.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 December 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//    N must be positive.
//
//    Input/output, double A[N*N].
//    On input, A is the coefficient matrix of the linear system.
//    On output, A is in unit upper triangular form, and
//    represents the U factor of an LU factorization of the
//    original coefficient matrix.
//
//    Input, double B[N], the right hand side of the linear system.
//
//    Output, double R8GE_FS[N], the solution of the linear system.
//
{
  int i;
  int ipiv;
  int j;
  int jcol;
  double piv;
  double t;
  double *x;

  x = new double[n];

  for ( i = 0; i < n; i++ )
  {
    x[i] = b[i];
  }

  for ( jcol = 1; jcol <= n; jcol++ )
  {
//
//  Find the maximum element in column I.
//
    piv = r8_abs ( a[jcol-1+(jcol-1)*n] );
    ipiv = jcol;
    for ( i = jcol+1; i <= n; i++ )
    {
      if ( piv < r8_abs ( a[i-1+(jcol-1)*n] ) )
      {
        piv = r8_abs ( a[i-1+(jcol-1)*n] );
        ipiv = i;
      }
    }

    if ( piv == 0.0 )
    {
      cout << "\n";
      cout << "R8GE_FS - Fatal error!\n";
      cout << "  Zero pivot on step " << jcol << "\n";
      return NULL;
    }
//
//  Switch rows JCOL and IPIV, and X.
//
    if ( jcol != ipiv )
    {
      for ( j = 1; j <= n; j++ )
      {
        t                 = a[jcol-1+(j-1)*n];
        a[jcol-1+(j-1)*n] = a[ipiv-1+(j-1)*n];
        a[ipiv-1+(j-1)*n] = t;
      }
      t         = x[jcol-1];
      x[jcol-1] = x[ipiv-1];
      x[ipiv-1] = t;
    }
//
//  Scale the pivot row.
//
    t = a[jcol-1+(jcol-1)*n];
    a[jcol-1+(jcol-1)*n] = 1.0;
    for ( j = jcol+1; j <= n; j++ )
    {
      a[jcol-1+(j-1)*n] = a[jcol-1+(j-1)*n] / t;
    }
    x[jcol-1] = x[jcol-1] / t;
//
//  Use the pivot row to eliminate lower entries in that column.
//
    for ( i = jcol+1; i <= n; i++ )
    {
      if ( a[i-1+(jcol-1)*n] != 0.0 )
      {
        t = - a[i-1+(jcol-1)*n];
        a[i-1+(jcol-1)*n] = 0.0;
        for ( j = jcol+1; j <= n; j++ )
        {
          a[i-1+(j-1)*n] = a[i-1+(j-1)*n] + t * a[jcol-1+(j-1)*n];
        }
        x[i-1] = x[i-1] + t * x[jcol-1];
      }
    }
  }
//
//  Back solve.
//
  for ( jcol = n; 2 <= jcol; jcol-- )
  {
    for ( i = 1; i < jcol; i++ )
    {
      x[i-1] = x[i-1] - a[i-1+(jcol-1)*n] * x[jcol-1];
    }
  }

  return x;
}