/* SOHeatTransf -- Simulates the process of heat transference. */
/* Last edited on 2003-09-15 17:28:31 by ra014852 */


#include <SOGrid.h>
#include <SOMatrix.h>
#include <SOFunction.h>
#include <SOLinCombFunction.h>
#include <SOProcFunction.h>
#include <SOApprox.h>
#include <SOFuncMap.h> 
#include <SOParams.h>
#include <SOBasic.h>
#include <SOPlotParams.h>
#include <SOPlot.h>
#include <SO2DPlot.h>

#include <rn.h>
#include <js.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <values.h>
#include <limits.h>


typedef struct Options
  { char *basisName;      /* Prefix for basis filenames (minus extension). */
    char *matName;        /* Prefix for matrix filenames (default {basisName}). */
    char *outName;        /* Prefix for output filenames. */
    char *initfName;      /* Prefix for the initial temperature function filename. */
    char *finalfName;     /* Prefix for the final temperature function filename. */
    double tIni;       /* Initial simulation time (in seconds). */
    double tFin;       /* Total simulation time (in seconds). */
    double tStep;      /* Time interval between global iterations (in seconds). */
    nat gaussOrder;    /* Gaussian quadrature order for dot products. */

    /* Method used to solve the linear system: */
    bool gaussElim;    /* TRUE uses the generic Gauss elimination method. */
    bool gaussSeidel;  /* TRUE uses Gauss-Seidel iteration. */
    bool conjGrad;     /* TRUE uses the conjugate gradient method. */
    double omega;      /* Relaxation factor for Gauss-Seidel. */
    /* Plotting options: */
    bool plot;         /* TRUE to plot results. */
    bool plotsub;      /* TRUE to add subtitles while plotting results. */
    PlotOptions plt;   /* Plot format and style options. */
  } Options;

Options GetOptions(int argn, char **argc);

void open_plot
  (
    SOPlot_PSStream **temp_ps,
    SOPlot_PSStream **dtemp_ps,
    Options *o,
    int *meshDepth, 
    double *minR, double *maxR,
    int *totCount
  );

int main(int argn, char **argc)
  { Options o = GetOptions(argn, argc);
    int i;
    double time, K = 1.2; /* K: heat diffusivity coefficient. */
    char *init_time, *end_time;

    printf("DBG-> o.basisName: %s \n", o.basisName);

    /* Gets basis and domain dimension */
    char *basName = txtcat(o.basisName, ".bas");
    Basis bas = SOApprox_ReadBasis(basName);
    unsigned int dim = bas.nel;
    SOGrid_Dim pDim = bas.el[0]->d->pDim;
    double maxP[pDim];

    /* Gets the {SOGrid} tree */
    char treeFile[100];
    strcpy(treeFile, o.basisName); i = strlen(treeFile);

    treeFile[i - 5] = treeFile[i - 1]; // max rank number.
    if(treeFile[i - 1] == '0' || treeFile[i - 1] == '2')
      { treeFile[i - 6] = treeFile[i - 2];  treeFile[i - 5] = treeFile[i - 1]; } // max rank number.

    treeFile[i - 4] = '.'; treeFile[i - 3] = 't'; treeFile[i - 2] = 'r';
    treeFile[i - 1] = 'e'; treeFile[i] = 'e'; treeFile[i + 1] = '\0';
    FILE *rd = open_read(treeFile);
    SOGrid_Tree *tree = SOGrid_Tree_read(rd);

    double_vec T = double_vec_new(dim);        /* Coeffs for the temperature approximation function. */
    double_vec dT = double_vec_new(dim);       /* Coeffs for the differential temperature approximation function. */
    double_vec b = double_vec_new(dim);        /* Right-hand side vector. */

    double TotH[(dim + 1) * (dim)];
    SOMatrix H = SOMatrix_Null(dim,dim);  /* Left hand side: Dot product matrix * area size. */
    SOMatrix Hb = SOMatrix_Null(dim,dim); /* Gradient dot product matrix. */
    
    SOFunction *temperature;   /* Current temperature approximation function. */
    SOFunction *dtemperature;  /* Current differential temperature approximation function. */
 
    int iter = 0, meshDepth, plot_count_dtemp = 0, plot_count_temp = 0, totCount;
    double minR[2], maxR[2];   /* Rectangle to plot. */
    double tempfRange = o.plt.fRange, tempfStep = (double)o.plt.fRange / 10;
    double dtempfRange = o.plt.fRange, dtempfStep = (double)o.plt.fRange / 10;
    double tempMin, tempMax = 0.0, plottempMax = 0.0;
    double dtempMin, dtempMax = 0.0;
    SOPlot_PSStream *temp_ps, *dtemp_ps;
    FuncMap NoMap = IdentityMap(1, pDim);

    open_plot(&temp_ps, &dtemp_ps, &o, &meshDepth, minR, maxR, &totCount);

    auto void plot_state
      ( SOFunction *temp, SOFunction *dtemp, 
        double t, char *tag, int iter, double *plottempMax
      );

    void plot_state
      ( SOFunction *temp, SOFunction *dtemp, 
        double t, char *tag, int iter, double *plottempMax
      )
      {
        char tmpstr[200]; 
        if (temp != NULL)
          { 
            double tPlotMin = 0.0, tPlotMax = 0.0;    
            sprintf(tmpstr,"t: %.6f - max: %.4f - min: %.4f", t, tempMax, tempMin);

            if(iter == 0) tempfStep = (double)tempfRange / 8.0; // 10; 
            else {
                   tempfStep = (double)tempMax / 8.0; // 10; 
                   if(tempfStep < 1.0)tempfStep = 1.0;
                 }
  
            /* Plot extra do erro -- manter padrao imagens. */ 
            if(*tag == 'e'){tempfStep =  tempfStep / 10.0; printf("** tempfStep ERRO: %f \n\n", tempfStep);}

            printf("** tempfRange: %f, tempfStep: %f \n\n", tempfRange, tempfStep);

            o.plt.fRange = tempfRange;  o.plt.fStep = tempfStep;
            SO2DPlot_AddFunctionPicture ( temp_ps, temp, &NoMap, tree, minR, maxR, -tempfRange, +tempfRange,
                                          tempfStep, o.plt.isolineWidth, meshDepth - 6, 2, meshDepth,
                                          o.plt.gridWidth, (! o.plt.noGrid), 20,
                                          (! o.plt.noBands), (! o.plt.noIsolines), &tPlotMin, &tPlotMax,
                                          o.plotsub?tmpstr:(char*)NULL
                                        );
            *plottempMax = tPlotMax;

/*             if(iter == 1) */
/*               { tempfRange = fabs(tPlotMax) > fabs(tPlotMin) ? fabs(tPlotMax) : fabs(tPlotMin); */
/*                 tempfStep = (double)tempfRange / 15;  */
/*               } */
            plot_count_temp++;
            if(!(plot_count_temp % totCount))
              { sprintf(tmpstr,"%08d", iter);
                SOPlot_AddPage(temp_ps, tmpstr);
              }

            printf("observed TEMPERATURE function extrema, at TIME %f :", t);
            printf(" min = %16.12f max = %16.12f\n\n", tPlotMin, tPlotMax);
          }
        if (dtemp != NULL)
          {
            double dtPlotMin = 0.0, dtPlotMax = 0.0;
            sprintf(tmpstr,"t: %.6f - max: %.4f - min: %.4f", t, dtempMax, dtempMin);

            if(iter == 0) dtempfStep = (double)dtempfRange / 10; 
            else dtempfStep = (double)dtempMax / 10;
            if(iter == 1) dtempfRange = fabs(dtempMax);

            printf("** iter: %d, dtempfRange: %f, dtempfStep: %f \n\n", iter, dtempfRange, dtempfStep);

/*             SO2DPlot_AddFunctionPicture ( dtemp_ps, dtemp, &NoMap, tree, minR, maxR, -dtempfRange, +dtempfRange, */
/*                                           dtempfStep, o.plt.isolineWidth, meshDepth - 6, 2, meshDepth, */
/*                                           o.plt.gridWidth, (! o.plt.noGrid), 20, */
/*                                           (! o.plt.noBands), (! o.plt.noIsolines), &dtPlotMin, &dtPlotMax, */
/*                                           o.plotsub?tmpstr:(char*)NULL */
/*                                         ); */
            if(iter == 1)
              { dtempfRange = fabs(dtPlotMax) > fabs(dtPlotMin) ? fabs(dtPlotMax) : fabs(dtPlotMin);
                dtempfStep = (double)dtempfRange / 15;
              }
            plot_count_dtemp++;
            if(!(plot_count_dtemp % totCount))
              { sprintf(tmpstr,"%08d", iter);
                SOPlot_AddPage(dtemp_ps, tmpstr);
              }

            printf( "observed DIFFERENTIAL TEMPERATURE function extrema, at TIME %f :", t);
            printf( " min = %16.12f max = %16.12f\n", dtPlotMin, dtPlotMax);
            printf( "\n");
          }
      }   

    
    /* Generate the initial temperature fields {T}: */
    if(!strcmp(o.initfName, o.basisName))
      { int sourceind; /* Index of the heat source. */
        double sourcepos[2] = {0.5, 0.5};                           
        SOApprox_GetBasisIndeces(bas, pDim, 1, &sourceind, sourcepos);
        for(i = 0; i < dim; i++) { T.el[i] = 0.0; dT.el[i] = 0.0; }
        T.el[sourceind] = 1000.0; 
      }
    else 
      { char *initfName = txtcat(o.initfName, ".fun");
        FILE *initf = open_read(initfName);
        SOLinCombFunction *initTemp = SOLinCombFunction_Cast(SOFunction_Read(initf));
	//        for(i = 0; i < dim; i++) { T.el[i] = (double)1000 * initTemp->d->coef.el[i]; dT.el[i] = 0.0; }
        for(i = 0; i < dim; i++){ T.el[i] = initTemp->d->coef.el[i]; dT.el[i] = 0.0; }
      }

    temperature = SOApprox_BuildIterFunction(bas, T, basName, pDim, 1);
    dtemperature = SOApprox_BuildIterFunction(bas, dT, basName, pDim, 1);

    /* Obtain the system's matrix (actually, the regular eval dot product matrix): */
    H = SOApprox_ReadMatrix(txtcat(o.matName, "-ev.mat"));
    Hb = SOApprox_ReadMatrix(txtcat(o.matName, "-gr.mat"));
    //    Hc = SOApprox_ReadMatrix(txtcat(o.matName, "-gr1.mat"));

    time = o.tIni;
    if (o.plot) plot_state(temperature, dtemperature, time, "i", 0, &plottempMax);

    /* Gets the initial experiment time. */
    init_time = today();
 
    while (time <= o.tFin)
      { printf( "=== time = %f ===============================\n", time);

        /* Given the temperature field {T} for this instant, computes its change {dT}. */

        /* Obtain the system's right-hand side: */
        for(i = 0; i < dim; i++) b.el[i] = 0.0;
        for(i = 0; i < Hb.ents.nel; i++) 
          b.el[Hb.ents.el[i].row] += (double)Hb.ents.el[i].va * T.el[Hb.ents.el[i].col]; 
        for(i = 0; i < dim; i++) b.el[i] = - K * b.el[i];

/*         for(i = 0; i < b.nel; i++) printf(" DEBUG: b[%d]->  %f \n", i, b.el[i]); */
 
        /* Solve the system obtaining the new guesses {AdBnew}: */
        if (o.gaussElim)
          { SOApprox_GaussElim(TotH, H, b, dT, FALSE); }
        else
          { assert(FALSE , "no solution method?"); }

/*         for(i = 0; i < dT.nel; i++) printf(" DEBUG: dT[%d]->  %f \n", i, dT.el[i]); */
/*         for(i = 0; i < T.nel; i++) printf(" DEBUG : T[%d]->  %f \n", i, T.el[i]); */
                
        /* The new temperature and differential temperature functions are generated. */
        for(i = 0; i < dim; i++) T.el[i] = T.el[i] + dT.el[i] * o.tStep;  
        
        SOApprox_GetLimitValues(temperature, maxP, &tempMax, &tempMin);
        SOApprox_GetLimitValues(dtemperature, maxP, &dtempMax, &dtempMin);
        printf( "Maximum TEMPERATURE, at TIME %f :", time);
        printf( " %16.12f with DIFFERENCE: %16.12f \n", tempMax, dtempMax);

        iter++;

        /* 2D Plotting */
        /* Plots the new temperature and differential temperature functions. */
        if (o.plot && (pDim == 2) && ((iter % 50) == 0 || iter < 10 || time > 0.0056)) 
          plot_state(temperature, dtemperature, time, "f", iter, &plottempMax);
          
        time += o.tStep;
      }    

    /* Gets the final experiment time. */
    end_time = today();
    printf( "INITIAL experiment TIME: %s \n", init_time);
    printf( "FINAL experiment TIME: %s \n\n", end_time);

    /* Calculates the distance between the simulation result and the analytical solution. */
    if(o.finalfName != NULL)
      {
        double distance = 0.0, ePlotMin = 0.0, ePlotMax = 0.0;
        char *finalfName = txtcat(o.finalfName, ".fun");
        FILE *finalf = open_read(finalfName);
        FILE *rd2 = open_read("../SOFindTentBasis/regular10.tree");
        SOLinCombFunction *finalTemp = SOLinCombFunction_Cast(SOFunction_Read(finalf));
        SOGrid_Tree *tree10 = SOGrid_Tree_read(rd2);
	//        SOFunction *finalTempf = SOFunction_Cast(finalTemp);

        /* Builds error function: */
        Basis error_bas = SOFunctionRef_vec_new(2);
        double_vec errc = double_vec_new(2); errc.el[0] = 1; errc.el[1] = -1;
        error_bas.el[0] = (SOFunction *)temperature; error_bas.el[1] = (SOFunction *)finalTemp;
        SOFunction *e = (SOFunction *)SOLinCombFunction_Make(pDim, 1,"Errorf.bas", error_bas, errc);

        o.plt.fRange = tempfRange / 10.0;  o.plt.fStep = tempfStep / 10.0;
        SOApprox_PlotSingleFunction(e, &NoMap, tree, txtcat(o.outName, "-err"), &(o.plt), &ePlotMin, &ePlotMax);

        /* Plot extra do erro -- manter padrao imagens. */ 
        tempMax = fabs(ePlotMax) > fabs(ePlotMin) ? fabs(ePlotMax) : fabs(ePlotMin);
        iter++; plot_state(e, dtemperature, time, "e", iter, &ePlotMax);
        iter++; plot_state(e, dtemperature, time, "e", iter, &ePlotMax);
        /* Plot extra do erro -- manter padrao imagens. */ 

        printf( "observed ERROR function extrema:");
        printf( " min = %16.12f max = %16.12f\n", ePlotMin, ePlotMax);
        printf( "\n");

        SOFunction_Dot(e, &NoMap, e, &NoMap, (SOFunction *)NULL, tree10, &distance, TRUE);
        distance = sqrt(distance);
        printf( "observed DISTANCE ( |f-g| = sqrt(<f-g, f-g>) ): %.16g \n", distance);
      }

    SOPlot_CloseStream(temp_ps);
    SOPlot_CloseStream(dtemp_ps);
    return 0;
  }
  

void open_plot
  (
    SOPlot_PSStream **temp_ps,
    SOPlot_PSStream **dtemp_ps,
    Options *o,
    int *meshDepth, 
    double *minR, double *maxR,
    int *totCount
  )
  {
    char paperSize[50];
    int pos = 0, hCount = 1, vCount = 1;

    if(!o->plt.eps)
      { 
        strcpy(paperSize,"a3"); 
        pos = 2; 
        hCount = 2; vCount = 3;
      }
    
    *totCount = hCount * vCount; 
    paperSize[pos] = '\0';    

    *temp_ps = SOPlot_OpenStream(o->plt.eps, txtcat("temp_", o->outName), 
                                 paperSize, 350.0, 350.0 * SQRTHALF);
    *dtemp_ps = SOPlot_OpenStream(o->plt.eps, txtcat("dtemp_", o->outName), 
                                  paperSize, 350.0, 350.0 * SQRTHALF);   

    *meshDepth = (int)ceil(2*log(o->plt.figSize/o->plt.meshSize)/log(2));
    minR[X] = 0.0;  maxR[X] = 1.0;
    minR[Y] = 0.0;  maxR[Y] = SQRTHALF; //1.0;

    SOPlot_ChangeLayout(*temp_ps, paperSize, 300.0, 300.0 * SQRTHALF,(o->plt.eps?0:1), hCount, vCount);
    SOPlot_ChangeLayout(*dtemp_ps, paperSize, 300.0, 300.0 * SQRTHALF,(o->plt.eps?0:1), hCount, vCount);
  }


#define PPUSAGE SOParams_SetUsage

Options GetOptions(int argn, char **argc)
  {
    Options o;
    SOParams_T *pp = SOParams_NewT(stderr, argn, argc);

    PPUSAGE(pp, "SOUNiSolve \\\n");
    PPUSAGE(pp, "  -basisName NAME [ -matName NAME ] \\\n");
    PPUSAGE(pp, "    -tIni TOTAL_TIME \\\n");
    PPUSAGE(pp, "    -tFin TOTAL_TIME \\\n");
    PPUSAGE(pp, "    -tStep TIME_INTERVAL \\\n");
    PPUSAGE(pp, "  [ -gaussSeidel [-omega NUM] | \\\n");
    PPUSAGE(pp, "    -conjGrad | \\\n");
    PPUSAGE(pp, "    -gaussElim ] \\\n");
    PPUSAGE(pp, "  [ -gaussOrder NUM ] \\\n");
    PPUSAGE(pp, "  [ -initfName NAME ] \\\n");
    PPUSAGE(pp, "  [ -finalfName NAME ] \\\n");
    PPUSAGE(pp, "  -outName NAME \\\n");
    PPUSAGE(pp, "  [ -plot ] \\\n");
    PPUSAGE(pp, "  [ -plotsub ] \\\n");
    PPUSAGE(pp, SOPlotParams_FunctionHelp " \n");

    SOParams_GetKeyword(pp, "-basisName");                               
    o.basisName = SOParams_GetNext(pp);  
       
    if (SOParams_KeywordPresent(pp, "-matName"))
      { o.matName = SOParams_GetNext(pp); }
    else
      { o.matName = o.basisName; }

    if (SOParams_KeywordPresent(pp, "-initfName"))
      { o.initfName = SOParams_GetNext(pp); }
    else
      { o.initfName = o.basisName; }

    if (SOParams_KeywordPresent(pp, "-finalfName"))
      { o.finalfName = SOParams_GetNext(pp); }
    else
      { o.finalfName = NULL; }

    SOParams_GetKeyword(pp, "-outName");                               
    o.outName = SOParams_GetNext(pp);  
                                                 
    o.gaussSeidel = SOParams_KeywordPresent(pp, "-gaussSeidel");
    o.conjGrad = SOParams_KeywordPresent(pp, "-conjGrad");
    o.gaussElim = SOParams_KeywordPresent(pp, "-gaussElim");
    if (! (o.gaussSeidel || o.conjGrad || o.gaussElim))
      { SOParams_Error(pp, "must specify a linear solution method"); }
    else if ((o.gaussSeidel + o.conjGrad + o.gaussElim) > 1)
      { SOParams_Error(pp, "please specify only one linear solution method"); }

    if (o.gaussSeidel)
      { if (SOParams_KeywordPresent(pp, "-omega"))
          { o.omega = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); }
        else
          { o.omega = 1.0; }
      }

    SOParams_GetKeyword(pp, "-gaussOrder");
    o.gaussOrder = SOParams_GetNextInt(pp, 1, INT_MAX);

    SOParams_GetKeyword(pp, "-tIni");
    o.tIni = SOParams_GetNextDouble(pp, 0.0, 1.0e8);  

    SOParams_GetKeyword(pp, "-tFin");                               
    o.tFin = SOParams_GetNextDouble(pp, 0.0, 1.0e8);  

    SOParams_GetKeyword(pp, "-tStep");                               
    o.tStep = SOParams_GetNextDouble(pp, 0.0, 1.0e8);

    /* Plotting options: */

    o.plot = SOParams_KeywordPresent(pp, "-plot");

    o.plotsub = SOParams_KeywordPresent(pp, "-plotsub");

    o.plt = SOPlotParams_FunctionParse(pp);
    
    SOParams_Finish(pp);
    return o;
  }