# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <suitesparse/umfpack.h>

# include "ppc_array.h"
# include "ppc_poisson.h"
# include "ppc_twb_quad.h"

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

void compute_element_stiffness ( struct elem *ep, struct TWB_qdat *qdat,
  int d, double (*f)( double x, double y ), double k[3][4] )

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

    compute_element_stiffness() computes an element stiffness matrix K.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct elem *ep: a structure definining the triangle as a finite element.

    struct TWB_qdat *qdat: a structure defining a TWB quadrature rule of
    a particular order.

    int d: the quadrature strength.

    double (*f)( double x, double y ): evaluates the integrand.

  Output:

    double k[3][4]: the element stiffness matrix, and the force vector 
    appended as the final column.
*/
{
  double lambda[3];
  int i;
  int j;
  double x[3];
  double X;
  double y[3];
  double Y;
/*
  Set the stiffness matrix entries.
*/
  for ( i = 0; i < 3; i++ )
  {
    for ( j = 0; j < 3; j++ )
    {
      k[i][j] = 0.25 * ( ep->ex[i] * ep->ex[j] + ep->ey[i] * ep->ey[j] ) 
        / ep->area;
    }
  }
/*
  Set the force vector.
*/
  for ( i = 0; i < 3; i++ )
  {
    x[i] = ep->n[i]->x;
    y[i] = ep->n[i]->y;
  }

  for ( i = 0; i < 3; i++ )
  {
    k[i][3] = 0.0;
  }

  for ( j = 0; j < d; j++ )
  {
    lambda[0] = qdat[j].lambda1;
    lambda[1] = qdat[j].lambda2;
    lambda[2] = qdat[j].lambda3;
    X = 0.0;
    Y = 0.0;
    for ( i = 0; i < 3; i++ )
    {
      X = X + lambda[i] * x[i];
      Y = Y + lambda[i] * y[i];
    }

    for ( i = 0; i < 3; i++ )
    {
      k[i][3] = k[i][3] + qdat[j].weight * f ( X, Y ) * lambda[i];
    }
  }

  for ( i = 0; i < 3; i++ )
  {
    k[i][3] = k[i][3] * ( ep->area / TWB_STANDARD_AREA );
  }

  return;
}
/******************************************************************************/

void enforce_zero_dirichlet_bc ( struct elem *ep, double k[3][4] )

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

    enforce_zero_dirichlet_bc() modifies an element matrix for a boundary condition.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    25 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int col;
  int i;
  int row;

  for ( i = 0; i < 3; i++ )
  {
    if ( ep->n[i]->bc == FEM_BC_DIRICHLET )
    {
      for ( col = 0; col <= 3; col++ )
      {
        k[i][col] = 0.0;
      }
      for ( row = 0; row <= 2; row++ )
      {
        k[row][i] = 0.0;
      }
      k[i][i] = 1.0;
    }
  }
  return;
}
/******************************************************************************/

void error_and_exit ( int status, const char *file, int line )

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

    error_and_exit() reports an UMFPACK error and terminates execution.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    25 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    int status: the UMFPACK error status code.

    const char *file: the file in which the error occurred.

    int line: the line number in which the error occurred.
*/
{
  fprintf ( stderr, "\n" );
  fprintf ( stderr, "ppc_poisson(): Fatal error!\n" );
  fprintf ( stderr, "  umfpack() returned error status code %d\n", status );
  fprintf ( stderr, "  See file %s, line %d\n", file, line );

  exit ( EXIT_FAILURE );
}
/******************************************************************************/

struct errors element_errors ( struct problem_spec *spec, 
  struct TWB_qdat *qdat, struct elem *ep, int d )

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

    element_errors() computes the error in an element.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: a specification of the problem;

    struct TWB_qdat *qdat: the TWB quadrature rule for a triangle.

    struct elem *ep: an element structure.

    int d: the strength of the quadrature rule.

  Output:

    struct errors element_errors(): the error norms comparing the computed 
    and exact solutions within a given element.
*/
{
  double energy;
  struct errors errors;
  double fxy;
  int i;
  int j;
  double L2norm;
  double lambda[3];
  double Linfty;
  double uex;
  double uh;
  double w;
  double x;
  double y;

  L2norm = 0.0;
  Linfty = 0.0;
  energy = 0.0;

  for ( j = 0; j < d; j++ )
  {
    lambda[0] = qdat[j].lambda1;
    lambda[1] = qdat[j].lambda2;
    lambda[2] = qdat[j].lambda3;

    w = qdat[j].weight;

    x = 0.0;
    y = 0.0;
    uh = 0.0;
    for ( i = 0; i < 3; i++ )
    {
      x =  x  + lambda[i] * ep->n[i]->x;
      y =  y  + lambda[i] * ep->n[i]->y;
      uh = uh + lambda[i] * ep->n[i]->z;
    }

    uex = spec->u_exact ( x, y );
    fxy = spec->f ( x, y );

    L2norm = L2norm + w * pow ( uex - uh, 2 );
    Linfty = fmax ( Linfty, fabs ( uex - uh ) );
    energy = energy + w * ( uex - uh ) * fxy;
  }
  
  errors.L2norm = L2norm * ep->area;
  errors.Linfty = Linfty;
  errors.energy = energy * ep->area;

  return errors;
}
/******************************************************************************/

struct errors eval_errors ( struct problem_spec *spec, struct mesh *mesh,
  int d )

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

    eval_errors() reports the error between exact and computed solutions.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: a specification of the problem;

    struct mesh *mesh: a definition of the mesh;

    int d: the desired quadrature strength.

  Output:

    struct errors eval_errors(): the error norms comparing the computed 
    and exact solutions.
*/
{
  struct errors errors;

  struct TWB_qdat *qdat = twb_qdat ( &d, NULL );
 
  errors.Linfty = 0.0;
  errors.L2norm = 0.0;
  errors.energy = 0.0;

  for ( int i = 0; i < mesh->nelems; i++ )
  {
    struct elem *ep = &mesh->elems[i];
    struct errors elem_errors = element_errors ( spec, qdat, ep, d );

    errors.Linfty = fmax ( errors.Linfty, elem_errors.Linfty );
    errors.L2norm = errors.L2norm + elem_errors.L2norm;
    errors.energy = errors.energy + elem_errors.energy;
  }

  errors.L2norm = sqrt ( errors.L2norm );
  errors.energy = sqrt ( errors.energy );

  return errors;
}
/******************************************************************************/

void poisson_solve ( struct problem_spec *spec, struct mesh *mesh, int d )

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

    poisson_solve() solves a poisson problem using the finite element method.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: a specification of the problem;

    struct mesh *mesh: a definition of the mesh;

    int d: the desired quadrature strength.
*/
{
  int *Ai;
  int *Ap;
  double *Ax;
  double *F;
  int i;
  int j;
  double k[3][4];
  void *Numeric = NULL;
  int r;
  int s;
  int status;
  void *Symbolic = NULL;
  int *Ti;
  int *Tj;
  struct TWB_qdat *qdat = twb_qdat ( &d, NULL );
  double *Tx;
  double *U;

  make_vector ( Ti, 3 * 3 * mesh->nelems );
  make_vector ( Tj, 3 * 3 * mesh->nelems );
  make_vector ( Tx, 3 * 3 * mesh->nelems );

  make_vector ( Ap, 1     + mesh->nnodes );
  make_vector ( Ai, 3 * 3 * mesh->nelems );
  make_vector ( Ax, 3 * 3 * mesh->nelems );

  make_vector ( F, mesh->nnodes );
  make_vector ( U, mesh->nnodes );
/*
  Define right hand side F,
  and i, j, Aij for finite element matrix.
*/
  for ( i = 0; i < mesh->nnodes; i++ )
  {
    F[i] = 0.0;
  }

  s = 0;
  for ( r = 0; r < mesh->nelems; r++ )
  {
    struct elem *ep = &mesh->elems[r];
    compute_element_stiffness ( ep, qdat, d, spec->f, k );
    enforce_zero_dirichlet_bc ( ep, k );

    for ( i = 0; i < 3; i++ )
    {
      int I = ep->n[i]->nodeno;
      for ( j = 0; j < 3; j++ )
      {
        if ( k[i][j] != 0.0 )
        {
          int J = ep->n[j]->nodeno;
          Ti[s] = I;
          Tj[s] = J;
          Tx[s] = k[i][j];
          s = s + 1;
        }
      }
      F[I] = F[I] + k[i][3];
    }
  }
/*
  Convert i, j, Aij to sparse column storage.
*/
  status = umfpack_di_triplet_to_col ( mesh->nnodes, mesh->nnodes, s,
    Ti, Tj, Tx, Ap, Ai, Ax, NULL );

  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }
/*
  Print matrix statistics.
*/
  printf ( "  Stiffness matrix has %d rows and %d columns\n", 
    mesh->nnodes, mesh->nnodes );
  printf ( "  and has %d nonzero entries\n", Ap[mesh->nnodes] );
/*
  Symbolic analysis.
*/
  status = umfpack_di_symbolic ( mesh->nnodes, mesh->nnodes, 
    Ap, Ai, Ax, &Symbolic, NULL, NULL );

  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }
/*
  Numeric analysis.
*/
  status = umfpack_di_numeric ( Ap, Ai, Ax, Symbolic, &Numeric, NULL, NULL );

  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }
/*
  Solve system.
*/
  status = umfpack_di_solve ( UMFPACK_A, Ap, Ai, Ax, U, F, Numeric, NULL, NULL );

  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }

  for ( i = 0; i < mesh->nnodes; i++ )
  {
    mesh->nodes[i].z = U[i];
  }
/*
  Release memory.
*/
  free_vector ( Ti );
  free_vector ( Tj );
  free_vector ( Tx );
  free_vector ( Ap );
  free_vector ( Ai );
  free_vector ( Ax );
  free_vector ( F );
  free_vector ( U );
  umfpack_di_free_symbolic ( &Symbolic );
  umfpack_di_free_numeric ( &Numeric );

  return;
}