/* Last edited on 2011-05-16 16:59:04 by stolfilocal */
/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *   Last edited on 2009-10-27 15:55:20 by stolfilocal
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include "arvore.h"
#include "jsmath.h"
#include "bool.h"
#include "affirm.h"
#include "definicoes.h"
#include "inicializa.h"
#include "interp_integ.h"
#include "gnuplot_i.h"
#include "imprime.h"
#include "timestep.h"
#include "plota.h"
#include "amr.h"
#include "bases.h"
#include "jstime.h"

#define DRAW_STEP 50
  /* Desenhar a arvore a cada {DRAW_STEP} iteracoes. */
  
#define PLOT_STEP 50
  /* Plotar a cada {PLOT_STEP} iteracoes. */

#define OUT_PREFIX "out"
#define OUT_PREFIX2 "outv"
#define OUT_PREFIX3 "outw"

void evoluiRK2(int d, bool_t op, bool_t ow, int N, int niv_tot, int tp, double eps, double Dt, double xmin[], double xmax[], VNo pac[], Funcao f, Fluxo flx, bool_t haar, bool_t podar);
void evoluiRK3(int d, bool_t op, bool_t ow, int N, int niv_tot, int tp, double eps, double Dt, double xmin[], double xmax[], VNo pac[], Funcao f, Fluxo flx, bool_t haar, bool_t podar);

double *read_vector(FILE *rd);

int main(int argc, char *argv[])
{
  bool_t op = FALSE; // ordem de previsao T <--> 4 , F <--> 2
  bool_t ow = FALSE; // ordem de previsao weno T <--> 3 , F <--> 2
  //int d = 3;
  //int d = 2;  /* Dimensao da grade. */
  int d = 1;  /* Dimensao da grade. */
  double xmin[d];  /* Coords inferiores da celula raiz. */
  double xmax[d];  /* Coords superiores da celula raiz. */
  int i;
  for (i = 0; i < d; i++) { xmin[i] = -1.0; xmax[i] = +1.0; }
  
  int niv_max = d*10;      /* Nivel maximo sem contar camada de seguranca. */
  int niv_tot = niv_max+1; /* Nivel maximo incluindo camada de seguranca. */

  int tp = op ? 9 : 5;   /* Tamanho do pacote em cada dimensao. */
  int tm = ipow(tp, d);  /* Tamanho total de um pacote. */
  int kc = tm/2;         /* Indice do elemento central do pacote. */
  //Funcao *u0 = &condicaoinicial;  /* Funcao de {x} que define o estado para {t=0}. */
  //Funcao *u0 = &gaussinicial; /* Funcao de {x} que define o estado para {t=0}. */
  //Funcao *u0 = &pulso;  /* Funcao de {x} que define o estado para {t=0}. */
  //Funcao *u0 = &seno;  /* Funcao de {x} que define o estado para {t=0}. */
  Funcao *u0 = &riemann;  /* Funcao de {x} que define o estado para {t=0}. */
  //Funcao *u0 = &riemann2;  /* Funcao de {x} que define o estado para {t=0}. */ 
  //Funcao *u0 = &riemann3;  /* Funcao de {x} que define o estado para {t=0}. */
  //Funcao *u0 = &bseno;  /* Funcao de {x} que define o estado para {t=0}. */
  
  Fluxo *flx = &adveccao; /* Funcao de {u} que define o fluxo continuo. ----- Adveccao */
  //Fluxo *flx = &burger;      /* Funcao de {u} que define o fluxo continuo. ----- Burger */
  //Fluxo *flx = &diferente;   /* Funcao de {u} que define o fluxo continuo. ----- Diferente */
  //Fluxo *flx = &buckleyLeverett;   /* Funcao de {u} que define o fluxo continuo. ----- Buckley-Leverett */

  double eps = 0.1;  /* Tolerancia para poda da arvore. */
  
  bool_t podar; /* elimina ou nao os nos superfluos? */
  bool_t haar;  // analise usando wavelets de haar 
  
  /* algoritmos de analise multiresolucao */
  //aproxima_funcao(d, xmin, xmax, f, niv_tot, eps, tp, pac, integra, interpola_quad, podar, haar);
  //bool_t sim =  TRUE;  // despreza coeficientes irrelevantes
  //bool_t recs = TRUE; // reconstroi ?
  //amr( d, pac, niv_tot, xmin, xmax, 0, interpola_quad, eps, sim, haar, recs);
  //apaga_arvore(pac[kc].p);

  //DoEPSTests(OUT_PREFIX);
  //fprintf(stderr,"-- Tchau --\n");
  //return EXIT_SUCCESS;

  //MakeTree(OUT_PREFIX, pac[kc].p, niv_tot);
  //DoMesh(OUT_PREFIX, pac[kc].p, 0);
  
  /* Desenha as tendas da base diadica */
  /*
  DrawTentMin(OUT_PREFIX, pac, 2, 0, interpola_quad, xmin, xmax);
  DrawTentMax(OUT_PREFIX, pac, 2, 0, interpola_quad, xmin, xmax);
  */
  /* Desenha um spline generico sobre uma malha (nos vertices o valor eh dado por {f} )*/
  //splineExemplo(OUT_PREFIX, pac, d, 0, interpola_quad, xmin, xmax, f); /* com problema */
  
  /* Determina menor dimensao {Dmin} do dominio: */
  double Dmin = +INF;
  for (i = 0; i < d; i++)
  { 
    double Dx = xmax[i] - xmin[i];
    if (Dx < Dmin) { Dmin = Dx; }
  }
  //double c = 1.0; /* por enquanto, depois tem que mudar a forma de calcular o passo no tempo */
  /* Determina o passo no tempo {Dt}: */
  double T = 0.3; //pow((2.0*sqrt(3.0)),d-1)*Dmin/c;                  /* Tempo para uma volta completa. */
  int tam1 = ipow(2,(niv_tot+d-1)/d); /* Numero de celulas por eixo (arred. p. cima). */
  double Dx = Dmin/tam1;              /* Tamanho minimo da celula. */
  double cfl = 0.6;
  double Dt = cfl*Dx;               /* Passo no tempo para atravessar meia celula. */ 
  int N = T/Dt + 1; //ifloor(T/Dt + 0.5, 1);      /* Numero de passos para volta completa. */
  
  fprintf(stderr, "T =    %.3f\n", T);
  fprintf(stderr, "Dmin = %.3f\n", Dmin);
  fprintf(stderr, "Dx =   %.6f\n", Dx);
  fprintf(stderr, "Dt =   %.6f\n", Dt);
  fprintf(stderr, "N  =   %d\n", N);

  /* gera arquivo analitic-solution.txt, com dados para plotar grafico da solucao analitica
     no tempo  t = N*Dt */
  //VNo pac[tm];
  //solanalitica_advec( d, xmin, xmax, u0, niv_tot, N-1, Dt, tp, pac, integraSA, interpola);

  //Funcao *sol = &sol_burg_riem; //confira o t dessa funcao para ver se eh o mesmo la de baixo.
  //solanalitica_burg( d, xmin, xmax, sol, niv_tot, N-1, Dt, tp, pac, integraSA, interpola);
  //for(i=0;i<tm;i++){ apaga_arvore(pac[i].p); }

fprintf(stderr, " T = %f \n", (N-1)*Dt);

  /* evolucao para malha uniforme */
  int nniveis = 8;
  int niv_ref = niv_tot+nniveis-1;
  podar = FALSE;
  haar = FALSE;
  double cpu_time_used_uniforme, cpu_time_used_adapt;
  VNo pacu[tm];  // Representacao do estado corrente. Malha uniforme 
  double inicio = user_cpu_time_usec();
  //evoluiRK2( d, op, ow, N, niv_tot, tp, eps, Dt, xmin, xmax, pacu, u0, flx, haar, podar);
  evoluiRK3( d, op, ow, N, niv_ref, tp, eps, Dt, xmin, xmax, pacu, u0, flx, haar, podar);
  double fim = user_cpu_time_usec();
  cpu_time_used_uniforme = (fim - inicio)/CLOCKS_PER_SEC;
  fprintf(stderr,"cpu - uniforme - = %f\n", cpu_time_used_uniforme);
  system("mv \"solucao.txt\" \"sol-ref.txt\" ");

  //DoMesh(OUT_PREFIX, pacu[kc].p, 0);

  //DoMeshND( d, OUT_PREFIX, pacu[kc].p, 0, xmin, xmax, xmin, xmax);
  
  //return EXIT_SUCCESS;
  // copia folhas para vetor 
  int ncf = 0;
  conta_folhas(pacu[kc].p, &ncf);
  int pos = 0;
  double vet[ncf]; // vetor com solucao na malha uniforme
  copia_folhas(pacu[kc].p, vet, &pos);
  int o;
  for(o=0;o<tm;o++){ apaga_arvore(pacu[o].p); }
  //fprintf(stderr,"passou\n");

  //pos = 0;
  //double soma = 0.0;
  //erro_norma_l1_v2( vet, pac[kc].p, &pos, &soma);
  //double cvol = vol_cel_folha_malha_uniforme(d, niv_tot, xmin, xmax);
  //double erroeu = soma*cvol; /* erro entre exata e uniforme */
  //double erroeu1 = soma/tam1;
  //fprintf(stderr,"%s = %f = %f\n", "erro entre exata e uniforme ", erroeu, erroeu1);
  //for(o=0;o<tm;o++){ apaga_arvore(pac[o].p); }
  //fprintf(stderr,"passou\n");

//return 0;

  /* evolucao para malha adaptativa */
  VNo paca[tm];  // Representacao do estado corrente. Malha adaptativa
  podar = FALSE;
  haar = FALSE;
  FILE *dado = fopen("dados.txt", "w");

  char EP[10], G[10], ER[10], U[10], A[10], CP[10]; 
  sprintf(EP,"log(eps) ");
  sprintf(G,"ganho ");
  sprintf(ER,"log(erroau) ");
  sprintf(U,"cpu-unif ");
  sprintf(A,"cpu-adap ");
  sprintf(CP,"Compress - tf "); // so no tempo final
  cpu_time_seq(dado, EP, G, ER, U, A, CP);

  int j;
  for(j=0;j<nniveis;j++)
  {   
    fprintf(stderr," iteracao numero = %d \n",j+1);
    double comeca = user_cpu_time_usec();
    //evoluiRK2( d, op, ow, N, niv_tot, tp, eps, Dt, xmin, xmax, paca, u0, flx, haar, podar);
    //evoluiRK3( d, op, ow, N, niv_tot, tp, eps, Dt, xmin, xmax, paca, u0, flx, haar, podar);
    evoluiRK3( d, op, ow, N, niv_tot+j, tp, eps, Dt, xmin, xmax, paca, u0, flx, haar, podar);
    double termina = user_cpu_time_usec();
    cpu_time_used_adapt = (termina - comeca) / CLOCKS_PER_SEC;
    
    int ncfa = 0;
    conta_folhas(paca[kc].p, &ncfa);
    double compressao = 100.0 - ncfa*100.0/ncf;
    
    /* mandar reconstruir adaptativo */
    reconstroi( d, op, paca, niv_ref, xmin, xmax, 0, interpola, haar);
    /* calcular erro norma1 niv_tot_1 */
    pos = 0;
    double soma = 0.0;
    erro_norma_l1_v2( vet, paca[kc].p, &pos, &soma);
    double cvol = vol_cel_folha_malha_uniforme(d, niv_ref, xmin, xmax);
    double erroau = soma*cvol; /* erro entre adaptiva e uniforme */
/*  // tem que arrumar a rotina que calcula solucao analitica
    soma = 0.0;
    pos = 0;
    erro_norma_l1_v2( vetex, paca[kc].p, &pos, &soma);
    double erroaa = soma*cvol; // erro entre adaptiva e analitica 
*/
    double ganho = cpu_time_used_uniforme/cpu_time_used_adapt;
    cpu_time(dado, log10(eps), ganho, log10(erroau), cpu_time_used_uniforme, cpu_time_used_adapt, compressao);
    fprintf(stderr,"eps = %f\n",eps);
    eps /= 10.0;
    
    /* limpa memoria */
    apaga_arvore(paca[kc].p);
  }
  fclose(dado);
  
  
//  MakeTree(OUT_PREFIX, pac[kc].p, niv_tot);
  /* calcular erros norma niv_tot_infty e niv_tot_1 */
//   fprintf(stderr, "Erro L_infty = %f\n", errmax( u, ua, niv_tot));
//   fprintf(stderr, "Erro L_1 = %f\n", errmid( u, ua, niv_tot));

  
  fprintf(stderr,"-- Tchau --\n");

  return EXIT_SUCCESS;
}

void evoluiRK3(int d, bool_t op,bool_t ow, int N, int niv_tot, int tp, double eps, double Dt, double xmin[], double xmax[], VNo pac[], Funcao f, Fluxo flx, bool_t haar, bool_t podar)
{

  aproxima_funcao(d, xmin, xmax, f, niv_tot, eps, tp, pac, integra, interpola, podar, haar,op);
  fprintf(stderr,"Criou a arvore inicial\n");
  
  /* Plota sol inicial */
  char str[100];
  //sprintf(str,"%s%08d%s%06f", "solucao-n-eps", 0, "-", eps);
  sprintf(str,"solucao.txt");

  FILE *disco = fopen(str, "w");
  escreve_em_arquivo( d, 0, disco, pac, xmin, xmax, interpola, 1, FALSE,op);
  fclose(disco);
  fprintf(stderr, "plotando: solucao inicial\n");

  /* Inicia gnuplot */
  gnuplot_ctrl *h1;
  h1 = gnuplot_init();
  //plotando( d, h1);
  //exportando( d, h1);
  
  int n;
  for(n=1;n<N;n++)
  {
    //fprintf(stderr, "============= n = %d ============\n", n);
    /* refina cada folha em um nivel */
    if(podar){refina( d, op, pac, niv_tot, xmin, xmax, 0, interpola, haar);}
    /* desenha a arvore refinada no tempo n */
    //if ((n <= 10) || (n == N-1) || (n % DRAW_STEP == 0))
    //{ MakeTree(OUT_PREFIX, pac[kc].p, niv_tot, n); }
    /* calcula fluxo numerico nas folhas */

    //calcula_fluxosT( d, 0, pac, xmin, xmax, interpola, flx, op);
    //calcula_fluxos_upwind( d, 0, pac, xmin, xmax, interpola, flx, op);
    calcula_fluxo_weno_v2( d, 0, pac, xmin, xmax, flx, interpola, ow, op);

    int kc = ipow(tp,d)/2;
    /* primeiro passo de Runge-Kutta 2 */
    avancaT13(d, 0, pac[kc].p, xmin, xmax, Dt, flx); 
    /* atualiza os nos internos da arvore */
    atualiza_folhas_v2( pac[kc].p, 0);
    atualiza_nos_internos(pac[kc].p);

    //calcula_fluxosT( d, 0, pac, xmin, xmax, interpola, flx, op); // fluxo de roe
    //calcula_fluxos_upwind( d, 0, pac, xmin, xmax, interpola, flx, op);
    calcula_fluxo_weno_v2( d, 0, pac, xmin, xmax, flx, interpola, ow, op);
      
    /* segundo passo de Runge-Kutta */
    avancaT23( d, 0, pac[kc].p, xmin, xmax, Dt, flx); 
    /* atualiza os nos internos da arvore */
    atualiza_folhas_v2( pac[kc].p, 1);
    atualiza_nos_internos(pac[kc].p);
    
    //calcula_fluxosT( d, 0, pac, xmin, xmax, interpola, flx, op); // fluxo de roe
    //calcula_fluxos_upwind( d, 0, pac, xmin, xmax, interpola, flx, op);
    calcula_fluxo_weno_v2( d, 0, pac, xmin, xmax, flx, interpola, ow, op);

    /* segundo passo de Runge-Kutta */
    avancaT33( d, 0, pac[kc].p, xmin, xmax, Dt, flx); 
    /* atualiza os nos internos da arvore */
    atualiza_folhas_v2( pac[kc].p, 2);
    atualiza_nos_internos(pac[kc].p);
    /* Remesh */
    /* poda */
    if(podar){poda( d, op, pac, niv_tot, xmin, xmax, 0, interpola, eps, podar, haar);}
    
    //grafica( d, n, N, PLOT_STEP, pac, interpola_quad, xmin, xmax, h1); // nao funciona, ainda 

    if(FALSE)// ((n <= 10) || (n == N-1) || (n % PLOT_STEP == 0))
    {
      /* mede erro entre solucao numerica e analitica na norma L1 e grava ni arquivo err */
      //double erro = 0.0;
      //erro_l1_analitica_numerica( d, 0, xmin, xmax, n*Dt, pac[kc].p, f, integra, &erro, c);
      //fprintf(err,"%f  %f\n", n*Dt, erro); 
      /* plot */
      //fprintf(stderr,"aqui\n");
      // so para mostrar na tela
      //sprintf(str,"solucao.txt");
      sprintf(str,"solucao-niv-%03d-n-%08d.txt", niv_tot, n);
  
      { FILE *risco = fopen(str, "w");
        escreve_em_arquivo( d, 0, risco, pac, xmin, xmax, interpola, 1, FALSE, op);
        fclose(risco);
      }

      // para que fique armazenado, caso eu precise
      //sprintf(str,"%s%08d%s%06f", "solucao-n-", n,"-eps-", eps);
      //FILE *cisco = fopen(str, "w");
      //escreve_em_arquivo( d, 0, cisco, pac, xmin, xmax, interpola, 1, FALSE, op);
      //fclose(cisco);

      //int nf = 0;
      //conta_folhas(pac[kc].p, &nf);
      //fprintf(stderr, "plotando: passo = %d folhas = %d\n",n,nf);
      
      //fprintf(comp,"%f  %f\n", n*Dt, 100.0-nf*100.0/pow(2,niv_tot)); 
      //exportando( d, h1);
      plotando(d,h1);
      //sleep(nf/100000 + 2);
      //MakeTree(OUT_PREFIX, pac[kc].p, niv_tot, n);
      //DoMesh(OUT_PREFIX, pac[kc].p, n);
      //DoMeshND( d, OUT_PREFIX, pac[kc].p, n, xmin, xmax, xmin, xmax);
    }
    
  }
  /* Grava solucao finaal: */
  { sprintf(str,"solucao-niv-%03d-final.txt", niv_tot);
    FILE *risco = fopen(str, "w");
    escreve_em_arquivo( d, 0, risco, pac, xmin, xmax, interpola, 1, FALSE, op);
    fclose(risco);
  }
  gnuplot_close(h1) ;
}


void evoluiRK2(int d, bool_t op,bool_t ow, int N, int niv_tot, int tp, double eps, double Dt, double xmin[], double xmax[], VNo pac[], Funcao f, Fluxo flx, bool_t haar, bool_t podar)
{

  aproxima_funcao(d, xmin, xmax, f, niv_tot, eps, tp, pac, integra, interpola, podar, haar,op);
  fprintf(stderr,"Criou a arvore inicial\n");
  
  /* Plota sol inicial */
  char str[100];
  //sprintf(str,"%s%08d%s%06f", "solucao-n-eps", 0, "-", eps);
  sprintf(str,"solucao.txt");

  FILE *disco = fopen(str, "w");
  escreve_em_arquivo( d, 0, disco, pac, xmin, xmax, interpola, 1, FALSE,op);
  fclose(disco);
  fprintf(stderr, "plotando: solucao inicial\n");

  /* Inicia gnuplot */
  gnuplot_ctrl *h1;
  h1 = gnuplot_init();
  //plotando( d, h1);
  //exportando( d, h1);
  
  int n;
  for(n=1;n<N;n++)
  {
    //fprintf(stderr, "============= n = %d ============\n", n);
    /* refina cada folha em um nivel */
    if(podar){refina( d, op, pac, niv_tot, xmin, xmax, 0, interpola, haar);}
    /* desenha a arvore refinada no tempo n */
    //if ((n <= 10) || (n == N-1) || (n % DRAW_STEP == 0))
    //{ MakeTree(OUT_PREFIX, pac[kc].p, niv_tot, n); }
    /* calcula fluxo numerico nas folhas */

      calcula_fluxosT( d, 0, pac, xmin, xmax, interpola, flx, op);
      //calcula_fluxos_upwind( d, 0, pac, xmin, xmax, interpola, flx, op);
      //calcula_fluxo_weno_v2( d, 0, pac, xmin, xmax, flx, interpola, ow, op);

    int kc = ipow(tp,d)/2;
    /* pacote de fluxos */
    VFl pc[d+1];
    int i;
    for(i=0;i<=d;i++)
    { 
      int k = i>1 ? i-1 : i;
      int pos = i!=1 ? kc - ipow(tp,k) : kc;
      pc[i].p = pac[pos].p;
      int dir;
      for(dir=0;dir<d;dir++)
      {
        pc[i].fvs[dir] = pac[pos].p->flxsup[dir];
        pc[i].fvi[dir] = pac[pos].p->flxinf[dir];
      }
    }
    /* primeiro passo de Runge-Kutta 2 */
    avancaTempoMeio(d, 0, pc, xmin, xmax, Dt, flx); /* aTM */
    /* atualiza as folhas da arvore (manda fv[0] <--> fval ) */
    atualiza_folhas(pac[kc].p);
    /* atualiza os nos internos da arvore */
    atualiza_arvore(pac[kc].p); /* com aTM */

      calcula_fluxosT( d, 0, pac, xmin, xmax, interpola, flx, op); // fluxo de roe
      //calcula_fluxos_upwind( d, 0, pac, xmin, xmax, interpola, flx, op);
      //calcula_fluxo_weno_v2( d, 0, pac, xmin, xmax, flx, interpola, ow, op);

    for(i=0;i<=d;i++)
    { 
      int k = i>1 ? i-1 : i;
      int pos = i!=1 ? kc - ipow(tp,k) : kc;
      pc[i].p = pac[pos].p;
      int dir;
      for(dir=0;dir<d;dir++)
      {
        pc[i].fvs[dir] = pac[pos].p->flxsup[dir];
        pc[i].fvi[dir] = pac[pos].p->flxinf[dir];
      }
    }
    /* segundo passo de Runge-Kutta */
    avancaTempoUm( d, 0, pc, xmin, xmax, Dt, flx);  /* aTU */
    /* Remesh */
    /* atualiza as folhas da arvore (manda fv[0] <--> fval ) */
    atualiza_folhas(pac[kc].p);
    /* atualisar os nos internos da arvore */
    atualiza_arvore(pac[kc].p);
    /* poda */
    if(podar){poda( d, op, pac, niv_tot, xmin, xmax, 0, interpola, eps, podar, haar);}
    
    //grafica( d, n, N, PLOT_STEP, pac, interpola_quad, xmin, xmax, h1); // nao funciona, ainda 

    if ((n <= 10) || (n == N-1) || (n % PLOT_STEP == 0))
    {
      /* mede erro entre solucao numerica e analitica na norma L1 e grava ni arquivo err */
      //double erro = 0.0;
      //erro_l1_analitica_numerica( d, 0, xmin, xmax, n*Dt, pac[kc].p, f, integra, &erro, c);
      //fprintf(err,"%f  %f\n", n*Dt, erro); 
      /* plot */
      //fprintf(stderr,"aqui\n");
      //char str[100];
      // so para mostrar na tela
      //sprintf(str,"solucao.txt");
      //FILE *risco = fopen(str, "w");
      //escreve_em_arquivo( d, 0, risco, pac, xmin, xmax, interpola, 1, FALSE, op);
      //fclose(risco);

      // para que fique armazenado, caso eu precise
      //sprintf(str,"%s%08d%s%06f", "solucao-n-", n,"-eps-", eps);
      //FILE *cisco = fopen(str, "w");
      //escreve_em_arquivo( d, 0, cisco, pac, xmin, xmax, interpola, 1, FALSE, op);
      //fclose(cisco);

      //int nf = 0;
      //conta_folhas(pac[kc].p, &nf);
      //fprintf(stderr, "plotando: passo = %d folhas = %d\n",n,nf);
      
      //fprintf(comp,"%f  %f\n", n*Dt, 100.0-nf*100.0/pow(2,niv_tot)); 
      //exportando( d, h1);
      //plotando(d,h1);
      //sleep(nf/100000 + 2);
      //MakeTree(OUT_PREFIX, pac[kc].p, niv_tot, n);
      //DoMesh(OUT_PREFIX, pac[kc].p, n);
      //DoMeshND( d, OUT_PREFIX, pac[kc].p, n, xmin, xmax, xmin, xmax);
    }
  }
  // so para mostrar na tela
  sprintf(str,"solucao.txt");
  FILE *risco = fopen(str, "w");
  escreve_em_arquivo( d, 0, risco, pac, xmin, xmax, interpola, 1, FALSE, op);
  fclose(risco);
  plotando(d,h1);
  gnuplot_close(h1) ;
}

double *read_vector(FILE *rd)
{
  int n;
  assert(fscanf(rd, "%d", &n) == 1);
  double *v = notnull(malloc(n*sizeof(double)), "no mem");
  int i;
  for (i = 0; i < n; i++) { assert(fscanf(rd, "%lg", &(v[i])) == 1); }
  return v;
}
  
      

/*
void solucao_analitica(int d, bool_t op, bool_t ow, int N, int niv_tot, int tp, double eps, double Dt, double xmin[], double xmax[], VNo pac[], Funcao f, Fluxo flx, bool_t haar, bool_t podar)
{ 

  auto double solanalitica(int d, int N, double Dt,  double x[], Funcao f);

  double solanalitica(int d, int N, double Dt,  double x[], Funcao f)
  {
    double xc = 0.0; // Supoe que dominio eh [-1,1]^d
    double dist = 0;
 
    double t = Dt*N;
  
    double a,b;
    int i;
    for (i = 0; i < d; i++)
    { 
      double px = x[i];
      double cx = xc + t;
      int lim = cx;
      if(lim%2==0){ a = lim - 1.0; b = lim + 1.0;}
      else{ a = lim; b = a + 2.0;}
      double n = -(b+a)/(b-a);
      double m = (1-n)/b;
      double d = m*px + n;
      double c = m*cx + n;
//     double dx = d - c;
      double dx = fabs(px + 2.0 - c) < fabs(px - c) ? fabs(px + 2.0 - c) : fabs(px - c);
//     double dx = px - c;
      dist += dx*dx;
    }
//   double r = 0.01; //Raio da gaussiana. //
//   double vf = 0.9*exp(-dist/r);
  // Gaussiana centrada em {xc,yc}  
      //dist = dx*dx + dy*dy;
      //r = 0.01;
      //vf = 3.0*exp(-dist/r);
  // caracteristica //
  double r = 0.5, vf;
  if( dist< r*r)
  { 
    vf = 1.0;
  }
  else
  { vf = 0.0; } 
  return vf;
  }
}
*/