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

# include "ppc_array.h"
# include "ppc_bvp.h"
# include "ppc_gauss_quad.h"
# include "ppc_twb_quad.h"

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

void apply_dirichlet_bc ( struct elem *ep, double (*g)(double x, double y ), 
  double k[3][4] )

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

    apply_dirichlet_bc() applies a Dirichlet boundary condition.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    12 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 elem *ep:

    double (*g)(double x, double y):

    double k[3][4]: The stiffness matrix.

  Output:

    double k[3][4]: The stiffness matrix after the Dirichlet boundary
    condition has been applied.
*/
{
  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 apply_neumann_bc ( struct elem *ep, double (*h)(double x, double y ), 
  double k[3][4], int gauss_d, struct Gauss_qdat *gauss_data )

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

    apply_neumann_bc() applies a Neumann boundary condition.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    12 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 elem *ep:

    double (*h)(double x, double y):

    double k[3][4]: The stiffness matrix.

    int gauss_d: the quadrature strength.

    struct Gauss_qdat *gauss_data: a structure defining a Gauss quadrature rule
    of a particular order.

  Output:

    double k[3][4]: The stiffness matrix after the Neumann boundary
    condition has been applied.
*/
{
  int i;
  int j;
  double L;
  double sum1;
  double sum2;
  double x;
  double x1;
  double x2;
  double xsi;
  double y;
  double y1;
  double y2;

  for ( i = 0; i < 3; i++ )
  {
    if ( ep->e[i]->bc == FEM_BC_NEUMANN )
    {
      x1 = ep->e[i]->n[0]->x;
      y1 = ep->e[i]->n[0]->y;
      x2 = ep->e[i]->n[1]->x;
      y2 = ep->e[i]->n[1]->y;
      L = sqrt ( pow ( x2 - x1, 2 ) + pow ( y2 - y1, 2 ) );

      sum1 = 0.0;
      sum2 = 0.0;
      for ( j = 0; j < gauss_d; j++ )
      {
        xsi = gauss_data[j].p;
        x = 0.5 * ( ( 1.0 - xsi ) * x1 + ( 1.0 + xsi ) * x2 );
        y = 0.5 * ( ( 1.0 - xsi ) * y1 + ( 1.0 + xsi ) * y2 );
        sum1 = sum1 + gauss_data[j].w
          * 0.5 * ( 1.0 - xsi ) * h ( x, y );
        sum2 = sum2 + gauss_data[j].w
          * 0.5 * ( 1.0 + xsi ) * h ( x, y );
      }

      k[(i+1)%3][3] = k[(i+1)%3][3] + 0.5 * L * sum1;
      k[(i+2)%3][3] = k[(i+2)%3][3] + 0.5 * L * sum2;
    }
  }
  return;
}
/******************************************************************************/

void bvp_solve ( struct problem_spec *spec, struct mesh *mesh, int twb_d,
  int gauss_d )

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

    bvp_solve() solves a boundary value 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 twb_d: the desired TWB quadrature strength.

    int gauss_d: the desired Gauss quadrature strength.
*/
{
  int *Ai;
  int *Ap;
  double *Ax;
  double *F;
  struct Gauss_qdat *gauss_data = gauss_qdat ( &gauss_d );
  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 *twb_data = twb_qdat ( &twb_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, twb_d, twb_data, spec->f, spec->eta, k );
    apply_neumann_bc ( ep, spec->h, k, gauss_d, gauss_data );
    apply_dirichlet_bc ( ep, spec->g, 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;
}
/******************************************************************************/

void compute_element_stiffness ( struct elem *ep, int twb_d, 
  struct TWB_qdat *twb_data, double (*f)( double x, double y ),
  double (*eta)( 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:

    12 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 elem *ep: a structure definining the triangle as a finite element.

    int twb_d: the quadrature strength.

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

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

    double (*eta)( double x, double y ): evaluates the diffusivity or
    conductivity.

  Output:

    double k[3][4]: the element stiffness matrix, and the force vector 
    appended as the final column.
*/
{
  double eta_average;
  double lambda[3];
  int i;
  int j;
  double x[3];
  double X;
  double y[3];
  double Y;
/*
  Compute average eta over element.
*/
  eta_average = 0.0;

  for ( j = 0; j < twb_d; j++ )
  {
    lambda[0] = twb_data[j].lambda1;
    lambda[1] = twb_data[j].lambda2;
    lambda[2] = twb_data[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];
    }
    eta_average = eta_average + twb_data[j].weight * eta ( X, Y );
  }
  eta_average = eta_average / ep->area;
/*
  Set the stiffness matrix entries.
*/
  for ( i = 0; i < 3; i++ )
  {
    for ( j = 0; j < 3; j++ )
    {
      k[i][j] = 0.25 * eta_average 
        * ( 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 < twb_d; j++ )
  {
    lambda[0] = twb_data[j].lambda1;
    lambda[1] = twb_data[j].lambda2;
    lambda[2] = twb_data[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] + twb_data[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 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_bvp(): 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, 
  int twb_d, struct TWB_qdat *twb_data, struct elem *ep )

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

    element_errors() computes the error in an element.

  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;

    int twb_d: the strength of the quadrature rule.

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

    struct elem *ep:

  Output:

    struct errors element_errors(): the error norms comparing the computed 
    and exact solutions within a given element.
*/
{
  struct errors errors;
  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;

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

    w = twb_data[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 );

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

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

struct errors eval_errors ( struct problem_spec *spec, struct mesh *mesh,
  int twb_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 twb_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 *twb_data = twb_qdat ( &twb_d, NULL );
 
  errors.Linfty = 0.0;
  errors.L2norm = 0.0;

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

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

  errors.L2norm = sqrt ( errors.L2norm );

  return errors;
}