#include <iostream>
#include <fstream>
#include <string>
#include <math.h>
#include <vector>

using namespace std;

double norme(vector<double> u,vector<double> v, int n)
{
    double p(0); 
    for(int i=1;i<=n;i++) p=(u[i]-v[i])*(u[i]-v[i])+p;
    return sqrt(p); 
}

int main()
{
    const int nn(1000);
    vector< double > x(nn+1);
    vector< double > u1(nn+1);
    vector< double > u0(nn+1);
    vector< double > f(nn+1);
    vector< double > vit(nn+1);
    vector< double > uxx(nn+1);
    vector< double > xnew(nn+1);
    vector< double > hloc(nn+1);
    vector< double > u1new(nn+1);
    string const nomFichier("results.dat");
    ofstream monFlux(nomFichier.c_str());

    int n(1);
    double xnu=0.01, xlam=1, time=10;
    double xmax=1, xmin=0;
    int iadaption(1), iteradap(30);
    double hmin(0.01), hmax(0.2), err(0.0001);

        for(int iadap = 1; iadap <= iteradap; iadap++)
        {
            if(iadap == 1) for (int i = 1; i <= n; i++) u1[i] = 0;
            if(iadaption == 0 || iadap == 1)
            {
                n = n+2;
                double h;
                h = (xmax-xmin)/(n-1);
                for(int i = 1; i <= n; i++) x[i] = xmin + (i-1)*h;
            }
            else if(iadaption == 1 && iadap > 1)
            {
                int nnew(1);
                xnew[nnew] = xmin;
                u1new[nnew] = u1[1];
                do
                {
                    for(int i = 1; i <= n-1; i++)
                    {
                        if(xnew[nnew] >= x[i] && xnew[nnew] <= x[i+1])
                        {
                            double hll(0);
                            hll = (hloc[i]*(x[i+1]-xnew[nnew])+hloc[i+1]*(xnew[nnew]-x[i]))/(x[i+1]-x[i]);
                            nnew = nnew + 1;
                            xnew[nnew] = min(xmax, xnew[nnew-1]+hll);
                            for(int j=1; j <= n-1; j++)
                            {
                                if(xnew[nnew] >= x[j] && xnew[nnew] <= x[j+1])
                                {
                                    u1new[nnew] = (u1[j]*(x[j+1]-xnew[nnew])+u1[j+1]*(xnew[nnew]-x[j]))/(x[j+1]-x[j]);
                                    goto fg;
                                }
                            }
                       }
                   }
                   fg:
                    {
                        continue;
                    }

                }while(xnew[nnew] < xmax);
                n = nnew;
                int w;
                for(w = 1; w <= nnew; w++)
                {
                    x[w] = xnew[w];
                    u1[w] = u1new[w];
                }
                cout << "iteradap = " << iadap << "new mesh = " << nnew << endl;

            }
            float f1, f2;
            for(int i = 1; i <= n; i++)
            {
                f1 = exp(-2000*(x[i]-0.3)*(x[i]-0.3));
                f2 = exp(-2000*(x[i]-0.6)*(x[i]-0.6));
                f[i] = 10*f1+20*f2;
                vit[i] = 1; //+fabs(sin(3*x[i]));
                //cout << "vit[i] = " << vit[i] << endl;
            }
            u1[1] = 1;
            for(int i = 1; i <= n; i++) u0[i] = u1[i];
            double t(0), kt(0);
            double som;

            do
            {
                kt = kt+1;
                double dt(1.e10);
                for(int i=2; i<=n-1; i++)
                {
                    double xip, xim, hh, dtloc;
                    xip = 0.5*(x[i+1]+x[i]);
                    xim = 0.5*(x[i]+x[i-1]);
                    hh = xip-xim;
                    //cout << "u0[" << i << "] = " << u0[i] << endl;
                    dtloc = 0.5*hh*hh/(hh*fabs(vit[i])+xnu+hh*hh*(fabs(f[i])+fabs(xlam*u0[i])));
                    dt = min(dt, dtloc);

                }
            t = t + dt;
            double ux;
            for(int i = 2; i <= n-1; i++)
            {
                double xip, xim;//, hh, dtloc;
                double uxp, uxm, uip, uim, xnuv;
                xip = 0.5*(x[i+1]+x[i]);
                //cout << 'x(2)'<< x[2] << endl;
                xim = 0.5*(x[i]+x[i-1]);
                uxp = (u0[i+1]-u0[i])/(x[i+1]-x[i]);
                uxm = (u0[i]-u0[i-1])/(x[i]-x[i-1]);
                uxx[i] = (uxp-uxm)/(xip-xim);
                uip = 0.5*(u0[i+1]+u0[i]);
                uim=0.5*(u0[i]+u0[i-1]);
                ux=(uip-uim)/(xip-xim);
                xnuv=xnu+0.25*(xip-xim)*fabs(vit[i]);
                u1[i]=u0[i]+dt*(xnuv*uxx[i]-vit[i]*ux+f[i]-xlam*u0[i]);
            }
            u1[n]=u1[n-1]-ux*(xmax-x[n-1])+uxx[n-1]*(xmax-x[n-1])*(xmax-x[n-1])/2;

            som = norme(u1,u0,nn);
            double som0;
            som0 = 1;
            if(kt == 1 && som > 1.e-6)  som0 = som;
            som = som/som0;
            for(int i = 1; i <= n; i++) u0[i]=u1[i];
             if(som < 1.e-6)
             {
                 goto stop;
                 //cout << "kt = " << kt << " Residual = " << som << endl;
             }

             }while(t<time);
            stop:
                {
                    cout << "kt = " << kt << "residual = " << som << endl;
                }
				
            uxx[1] = uxx[2];
            uxx[n] = uxx[n-1];

            for(int i = 1; i <= n; i++) hloc[i]=min(max(hmin,sqrt(err/(fabs(uxx[i])+1.e-6))),hmax);
        }
            for(int i = 1; i <= n; i++)
            {
                    monFlux << x[i] << "       " << u1[i] << "      " << f[i] << "       " << vit[i] << "    " << uxx[i] << "      " << hloc[i] << "        " << u1new[i] << endl;
            }
        return 0;
        }