# include <cmath>
# include <cstdlib>

using namespace std;

# include "poisson_2d_exact.hpp"

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

void poisson_2d_exact ( int n, double x[], double y[], double u[], double ux[], 
  double uy[], double uxx[], double uxy[], double uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    poisson_2d_exact() evaluates an exact solution of the Poisson equation.
//
//  Discussion:
//
//    The equation is:
//      - Uxx - Uyy = F(x,y)
//
//    The domain is the unit square:
//      Omega = 0 <= x, y <= 1
//
//    The solution is:
//      U(X,Y) = 2 * ( 1 + Y ) / ( ( 3 + X )^2 + ( 1 + Y )^2 )
//
//    The boundary conditions are determined by evaluating the exact
//    solution U(X,Y).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of points.
//
//    real X, Y: the coordinates of the points.
//
//  Output:
//
//    real U, Ux, Uy, Uxx, Uxy, Uyy: the solution and first and second
//    partial derivatives evaluated at the points (X,Y).
//
{
  int i;
 
 for ( i = 0; i < n; i++ )
  {
    u[i] = 2.0 * ( 1.0 + y[i] ) 
      / ( pow ( 3.0 + x[i], 2 ) + pow ( 1.0 + y[i], 2 ) );
    ux[i] = ( - 2 * x[i] - 6 ) * ( 2 * y[i] + 2 ) 
      / pow ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ), 2 );
    uy[i] = ( - 2 * y[i] - 2 ) * ( 2 * y[i] + 2 ) 
      / pow ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ), 2 ) 
      + 2 / ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ) );
    uxx[i] = 4 * ( y[i] + 1 ) * ( 4 * pow ( x[i] + 3, 2 ) 
      / ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ) ) - 1 ) 
      / pow ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ), 2 );
    uxy[i] = 4 * ( x[i] + 3 )* ( 4 * pow ( y[i] + 1, 2 ) 
      / ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ) ) - 1 ) 
      / pow ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ), 2 );
    uyy[i] = 4 * ( y[i] + 1 ) * ( 4 * pow ( y[i] + 1, 2 ) 
      / ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ) ) - 3 ) 
      / pow ( pow ( x[i] + 3, 2 ) + pow ( y[i] + 1, 2 ), 2 );
  }

  return;
}