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

# include "ppc_array.h"
# include "ppc_fetch_line.h"
# include "ppc_ll.h"
# include "ppc_nelder_mead.h"
# include "ppc_truss.h"
# include "ppc_truss_io.h"
# include "ppc_xmalloc.h"

# define BUFLEN 128

void   evaluate_reactions ( struct truss *truss );
void   evaluate_stresses ( struct truss *truss );
double link_stretch ( struct link *link );
int    solve_truss ( struct truss *truss, double h, double tol, int maxevals );
double stored_energy_function ( double lambda );
double stress_function ( double lambda );
double total_energy ( double *x, int n, void *params );

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

void evaluate_reactions ( struct truss *truss )

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

    evaluate_reactions() computes the reactions at support nodes.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double length;
  int np;
  conscell *p;
  conscell *q;
  double xi;
  double xj;
  double yi;
  double yj;

  for ( q = truss->nodes_list; q != NULL; q = q->next )
  {
    struct node *node = q->data;
/*
  Skip if node is not constrained.
*/
    if ( ( ! node->xfixed ) && ( ! node->yfixed ) )
    {
      continue;
    }
/*
  Find every link including this node as an endpoint.
*/
    else
    {
      np = node->n;
      xi = node->x;
      yi = node->y;
      for ( p = truss->links_list; p != NULL; p = p->next )
      {
        struct link *link = p->data;
        if ( link->np1->n == np )
        {
          xj = link->np2->x;
          yj = link->np2->y;
        }
        else if ( link->np2->n == np )
        {
          xj = link->np1->x;
          yj = link->np1->y;
        }
        else
        {
          continue;
        }
        length = sqrt ( pow ( xj - xi, 2 ) + pow ( yj - yi, 2 ) );
        node->rx = node->rx + link->A * link->stress * ( xj - xi ) / length;
        node->ry = node->ry + link->A * link->stress * ( yj - yi ) / length;
      }
    }
/*
  Add external force and reverse sign.
*/
    node->rx = - ( node->rx + node->fx );
    node->ry = - ( node->ry + node->fy );
  }

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

void evaluate_stresses ( struct truss *truss )

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

    evaluate_stresses() evaluates the stress at each link.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double lambda;
  conscell *p;
  double sigma;
  double value;
/*
  Traverse the link list.
*/
  for ( p = truss->links_list; p != NULL; p = p->next )
  {
    struct link *link = p->data;

    lambda = link_stretch ( link );
    sigma = stress_function ( lambda );
    value = link->E * sigma;
    link->stress = value;
  }
  return;
}
/******************************************************************************/

double link_stretch ( struct link *link )

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

    link_stretch() computes the relative stretching of a link.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double deformed_length;
  double dx;
  double dy;
  double reference_length;
  double value;

  reference_length = link->L;

  dx = link->np1->x - link->np2->x;
  dy = link->np1->y - link->np2->y;
  deformed_length = sqrt ( dx * dx + dy * dy );

  value = deformed_length / reference_length;

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

int solve_truss ( struct truss *truss, double h, double tol, int maxevals )

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

    solve_truss() computes the local minimum of the energy of the truss.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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 truss *truss: a structure defining the truss.

    double h: problem scale.

    double tol: error tolerance.

    int maxevals: maximum number of function evaluations.

  Output:

    int value: the local minimum of the truss energy.
*/
{
  int evalcount;
  int i;
  int j;
  int N = 2 * truss->nnodes;
  struct nelder_mead NM;
  conscell *p;
  double **s;
  double *x;
/*
  Copy the truss node locations into a vector x.
*/
  make_vector ( x, N );

  for ( j = 0, p = truss->nodes_list; p != NULL; p = p->next )
  {
    struct node *node = p->data;
    x[j++] = node->x;
    x[j++] = node->y;
  }
/*
  Create the simplex.
*/
  make_matrix ( s, N + 1, N );
  for ( i = 0; i < N + 1; i++ )
  {
    for ( j = 0; j < N; j++ )
    {
      s[i][j] = x[j];
    }
  }
  for ( j = 0; j < N; j++ )
  {
    s[j+1][j] = s[j+1][j] + h;
  }
/*
  Project the simplex over the constraint space.
*/
  for ( j = 0, p = truss->nodes_list; p != NULL; p = p->next )
  {
    struct node *node = p->data;
    if ( node->xfixed )
    {
      for ( i = 0; i < N + 1; i++ )
      {
        s[i][j] = node->X;
      }
    }
    j++;
    if ( node->yfixed )
    {
      for ( i = 0; i < N + 1; i++ )
      {
        s[i][j] = node->Y;
      }
    }
    j++;
  }
/*
  Define the Nelder-Mead structure.
*/
  NM.f = total_energy;
  NM.n = N;
  NM.s = s;
  NM.x = x;
  NM.h = h;
  NM.tol = tol;
  NM.maxevals = maxevals;
  NM.params = truss;
/*
  Invoke the Nelder-Mead function minimizer.
*/
  evalcount = nelder_mead ( &NM );

  free_vector ( x );
  free_matrix ( s );

  if ( maxevals < evalcount )
  {
    printf ( "\n" );
    printf ( "solve_truss(): Warning!\n" );
    printf ( "  No convergence after %d evaluations\n", evalcount );
    return 0;
  }
  else
  {
    printf ( "\n" );
    printf ( "solve_truss():\n" );
    printf ( "  Convergence after %d evaluations\n", evalcount );
    evaluate_stresses ( truss );
    evaluate_reactions ( truss );
    return 1;
  }
}
/******************************************************************************/

double stored_energy_function ( double lambda )

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

    stored_energy_function() evaluates the stored energy.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

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

  value = pow ( lambda, 4 ) / 24.0 + 1.0 / pow ( lambda, 2 ) / 12.0 - 1.0 / 8.0;

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

double stress_function ( double lambda )

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

    stress_function() evaluates the stress function.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    double lambda: the relative stretching of a link.

  Output:

    double value: the stress function for a given relative stretch.
*/
{
  double value;

  value = ( pow ( lambda, 3 ) - 1.0 / pow ( lambda, 3) );

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

double total_energy ( double *x, int n, void *params )

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

    total_energy() evaluates the total energy of a truss.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    04 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    double *x:

    int n:

    void *params:

  Output:

    double e: the total energy.
*/
{
  double e;
  int j;
  conscell *p;
  struct truss *truss = params;

  e = 0.0;
  j = 0;
  for ( p = truss->nodes_list; p != NULL; p = p->next )
  {
    struct node *node = p->data;
    node->x = x[j++];
    node->y = x[j++];
  }

  for ( p = truss->links_list; p != NULL; p = p->next )
  {
    struct link *link = p->data;
    e = e + link->E * link->A * link->L *
      stored_energy_function ( link_stretch ( link ) );
  }

  for ( p = truss->nodes_list; p != NULL; p = p->next )
  {
    struct node *node = p->data;
    e = e - ( node->x - node->X ) * node->fx
          - ( node->y - node->Y ) * node->fy;
  }

  return e;
}