# include <math.h>
# include <stdlib.h>

# include "poisson_2d_exact.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  Purpose:

    poisson_2d_exact() evaluates an exact solution of the Poisson equation.

  Discussion:

    The equation is:
      - Uxx - Uyy = 0

    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;
}