# include "fem.h"


femProblem *advdiffNew(double epsilon, double beta, double zeta, int size, int sizePlot)
{
    femProblem *myProblem = malloc(sizeof(femProblem));
    myProblem->epsilon = epsilon;
    myProblem->beta    = beta;
    myProblem->zeta    = zeta;

    myProblem->size = size;
    myProblem->X = malloc(sizeof(double) * size);
    myProblem->U = malloc(sizeof(double) * size);

    myProblem->sizePlot = sizePlot;
    myProblem->x = malloc(sizeof(double) * sizePlot);
    myProblem->u = malloc(sizeof(double) * sizePlot);
    myProblem->uregime = malloc(sizeof(double) * sizePlot);

    advdiffReset(myProblem);


    return myProblem;
}

void advdiffReset(femProblem *myProblem)
{
    double h;
    int i;

    h = 1.0/(myProblem->size-1);
    for (i=0; i < myProblem->size; i++) {
        myProblem->X[i] = i*h;
        myProblem->U[i] = 0.0; }
    myProblem->U[0] = 1.0;


    h = 1.0/(myProblem->sizePlot-1);
    for (i=0; i < myProblem->sizePlot; i++) {
        myProblem->x[i] = i*h;
        myProblem->u[i] = 0.0;
        myProblem->uregime[i] = 0.0; }
}

void advdiffFree(femProblem *myProblem)
{
    free(myProblem->X);
    free(myProblem->U);
    free(myProblem->x);
    free(myProblem->u);
    free(myProblem->uregime);
    free(myProblem);
}

void advdiffSteadySolution(femProblem *myProblem)
{
    int  size = myProblem->sizePlot, i;
    double *u = myProblem->uregime;
    double *x = myProblem->x;
    double beta    = myProblem->beta;
    double epsilon = myProblem->epsilon;
    double pe  = beta/epsilon;


    if (beta != 0.0)
       for (i=0; i < size; i++)
          u[i] = (exp(pe) - exp(pe*x[i]) )/(exp(pe) - 1.0);
    else
       for (i=0; i < size; i++)
          u[i] = 1 - x[i];

}

# ifndef NOANALYTIC

void advdiffSolution(femProblem *myProblem, double time)
{
    int  size = myProblem->sizePlot, i,j;
    double *u = myProblem->u;
    double *x = myProblem->x;
    double *uregime = myProblem->uregime;
    double beta = myProblem->beta;
    double epsilon = myProblem->epsilon;
    advdiffSteadySolution(myProblem);
    double prec=1.0;


  //
  //  Solution diffusion pure :-)
  //  A modifier !!!!
  //
  if (beta==0){
    for (i=0; i < size; i++)
        u[i] = 1 - x[i];
    for (j=1; j < 500; j++) {
        double k = M_PI * j;
        double alpha = k*k*epsilon;
        double coeff = 2/k;
        for (i=0; i < size; i++) {
                u[i] = u[i] - coeff * sin(k*x[i]) * exp(- alpha*time); }}
  }
  else
    {
        for (i=0; i < size; i++)
        u[i] = uregime[i];
        for (j=1; j < 500 && prec>1e-9; j++) {
            double k = M_PI * j;
            double coeff = beta/2.0;
            double coeff2 = coeff/epsilon;
            double coeff3 = (coeff*coeff)/epsilon;
            prec=2.0 * k * exp(coeff2 - coeff3*time - k*k*epsilon*time)/(k*k + coeff2*coeff2);
            for (i=0; i < size; i++) {
                double num = 2.0 * k * sin(k*x[i]) * exp(coeff2*x[i] - coeff3*time - k*k*epsilon*time);
                double den = k*k + coeff2*coeff2;
                u[i] = u[i] - (num/den);
            }
        }
    }


}

# endif
# ifndef NOTIMESTEP

double advdiffComputeTimeStep(femProblem *myProblem)
{
    int size = myProblem->size;
    double beta    = myProblem->beta;
    double epsilon = myProblem->epsilon;
    double zeta = myProblem->zeta ;

  //
  //  Solution pour Euler explicite :-)
  //  A modifier pour Runge-Kutta !!!!
  //


    double h = 1.0/(size-1);
    double dt = fmin( (zeta*h*beta + 2*epsilon)/(beta*beta), (h*h)/(zeta*h*beta + 2*epsilon));

    return 1.35*dt;
}

# endif
# ifndef NORIGHTHANDSIDE


void advdiffComputeRightHandSide(femProblem *myProblem, double *U, double *F)
{
    int  size = myProblem->size, i;
    double epsilon = myProblem->epsilon;
    double beta = myProblem->beta;
    double zeta = myProblem->zeta;
    double h = 1.0/(size-1);

  // A modifier !!
    F[0] = 0.0;
    F[size-1] = 0.0;
    for (i=1; i < size-1; i++)
        {
            F[i] = (-(1-zeta)*beta*((U[i+1]-U[i-1])/(2.0*h))) - (zeta*beta*((U[i]-U[i-1])/h)) + (epsilon*((U[i+1]-(2.0*U[i])+U[i-1])/(h*h)));
        }
}


# endif
# ifndef NOEULER

void advdiffUpdateEuler(femProblem *myProblem, double dt)
{
    double *U = myProblem->U;
    int size = myProblem->size;
    double *F = malloc(sizeof(double)*size);
    advdiffComputeRightHandSide(myProblem, U, F);

    int i;
    for (i=0;i<size;i++){
        F[i] = F[i]*dt;
        U[i] = U[i] + F[i];
        }
    free(F);

 }

# endif
# ifndef NORK

void advdiffUpdateRungeKutta(femProblem *myProblem, double dt)
{
    double *U = myProblem->U;
    int size = myProblem->size;
    double *K1 = malloc(sizeof(double)*size);
    double *Uu = malloc(sizeof(double)*size);
    advdiffComputeRightHandSide(myProblem, U, K1);
    int i;
    for (i=0;i<size;i++){
        Uu[i] = U[i] + dt*K1[i]/2.0;
        }
    double *K2 = malloc(sizeof(double)*size);
    advdiffComputeRightHandSide(myProblem, Uu, K2);
    for (i=0;i<size;i++){
        Uu[i] = U[i] + dt*K2[i]/2.0;
        }
    double *K3 = malloc(sizeof(double)*size);
    advdiffComputeRightHandSide(myProblem, Uu, K3);
    for (i=0;i<size;i++){
        Uu[i] = U[i] + dt*K3[i];
        }
    double *K4 = malloc(sizeof(double)*size);
    advdiffComputeRightHandSide(myProblem, Uu, K4);
    for (i=0;i<size;i++){
        U[i] = U[i] + (dt/6.0)*(K1[i]+2.0*K2[i]+2.0*K3[i]+K4[i]);
        }

    free(Uu);
    free(K1);
    free(K2);
    free(K3);
    free(K4);

}

# endif