/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *                                                                         *
 *   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.             *
 ***************************************************************************/

#include "fluxo.h"

double flux(double u1, double u2, double u3, double c)
{
  double result;
  if(fabs(u3-u2) <= fabs(u2-u1)) { result = u3 - u2;}
  else { result = u2 - u1; }
  
  result = u2 + 0.5*result;
  result = c*result;
  return result;
}

/* para um caso mais geral tem que passar a funcao de fluxo 
   a constante c ta fazendo papel da funcao de fluxo f (edp) 
     u_t + f(u)_x = 0, f(u) = c*u */
double fluxo_roe(double u1, double u2, double u3, double u4, double div, Fluxo f)
{
  double ul = u2 + 0.5*minmod(u3 - u2, u2 - u1);
  double ur = u3 - 0.5*minmod(u4 - u3, u3 - u2);

  /* a constante c ta fazendo papel da funcao de fluxo f (edp) 
     u_t + f(u)_x = 0, f(u) = c*u */
  //  double ca = ul!=ur ? (f(ur) - f(ul))/(ur - ul) : div;
  double ca = fabs(ul-ur)>0.0000000000001 ? (f(ur) - f(ul))/(ur - ul) : div;
  
  double absA = fabs(ca);
  /* caso geral return 0.5*(f(ur)+f(ul)-fabs(A(ul,ur))*(ur - ul)); */
  return 0.5*( f(ur) + f(ul) - absA*(ur - ul) );
}

double fluxo_upwind(double u1, double u2, double df, Fluxo f)
{
  /* caso geral return 0.5*(f(ur)+f(ul)-fabs(A(ul,ur))*(ur - ul)); */
  double fu;
  if(df>=0){ fu = f(u1); }
  else{ fu = f(u2); }
  return  fu;
}

double minmod(double a, double b)
{
  double mm = (fabs(a)<=fabs(b)) ? a : b;
  return mm;
}

bool_t folhasounulos(int d, int ord, VNo pac[])
{
  int ctr = ipow(ord,d)/2;
  bool_t efon = TRUE;
  int i;
  for(i=0;i<d;i++)
  {
    int corr = ipow(ord,i);
    int t;
    for(t=0;t<2;t++)
    {
      int s = ipow(-1,t);
      //pode ser qqr filho, estou interessado em saber se os vizinhos sao folhas ou nulos
      //assert((pac[ctr+s*corr].p->fil[0] == NULL) == (pac[ctr+s*corr].p->fil[1] == NULL));
      No *pf = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
      if(pf != NULL) 
      {
        efon = FALSE;
      }
    }
  }
  return efon;
}

void calcula_fluxosT(int d, int prof, VNo pac[], double xmin[], double xmax[], Preditor pred, Fluxo f, int lim, bool_t op)
{
  int ord = op ? 9 : 5;
  int tam = ipow(ord,d);
  int ctr = tam/2;
  
  if(pac[ctr].p == NULL)
  {
    return;
  }

  No *pf0 = pac[ctr].p->fil[0];
  No *pf1 = pac[ctr].p->fil[1];
  assert((pf0 == NULL) == (pf1 == NULL));

  if((pf0 == NULL) && (prof < lim )){ return;} /* porque isso aqui mesmo ? */
  if( prof > lim ){ return;}  
  if((pf0 == NULL) &&  (prof == lim ))
  { //fprintf(stderr,"prof = %d lim = %d\n", prof, lim);
    int dir;
    for(dir=0;dir<d;dir++)
    {
      int corr = ipow(ord,dir);
      double u1 = pac[ctr-2*corr].p == NULL ? pac[ctr-2*corr].fv : pac[ctr - 2*corr].p->fval;
      double u2 = pac[ctr-1*corr].p == NULL ? pac[ctr-1*corr].fv : pac[ctr - 1*corr].p->fval;
      double u3 = pac[ctr+0*corr].p == NULL ? pac[ctr+0*corr].fv : pac[ctr + 0*corr].p->fval;
      double u4 = pac[ctr+1*corr].p == NULL ? pac[ctr+1*corr].fv : pac[ctr + 1*corr].p->fval;
      double u5 = pac[ctr+2*corr].p == NULL ? pac[ctr+2*corr].fv : pac[ctr + 2*corr].p->fval;
//      double df = pac[ctr].p->divf[dir];
//      double df = (f(u[ord/2+corr])-f(u[ord/2]))/(u[ord/2+1]-u[ord/2]);
      double h = fabs(xmax[dir]-xmin[dir]);
      double df = (f(u4)-f(u3))/(h);
      /* se todos os vizinhos sao folhas ou nulos */
      bool_t tvsfon = folhasounulos( d, ord, pac);
      if(tvsfon)
      {
        //fprintf(stderr,"tudo folha\n");
        pac[ctr].p->flxsup[dir] = fluxo_roe(u2, u3, u4, u5, df, f);
        pac[ctr].p->flxinf[dir] = fluxo_roe(u1, u2, u3, u4, df, f);
	double tmp1 = u3;
        double tmp2 = tmp1;
	double tmp3 = u4;
        double tmp4 = tmp3;
        if(tmp1!=tmp2)
        {
          if(pac[ctr-2*corr].p == NULL){fprintf(stderr,"NULO\n");}
          fprintf(stderr,"fs1 = %f\tfs2 = %f\n",tmp1,tmp2);
          fprintf(stderr,"fs1 = %f\tfs2 = %f\n",tmp1,tmp2);
          fprintf(stderr,"u3 = %f\tu3 = %f\n", pac[ctr+0*corr].fv , pac[ctr + 0*corr].p->fval);
          fprintf(stderr,"u4 = %f\tu4 = %f\n", pac[ctr+1*corr].fv , pac[ctr + 1*corr].p->fval);
        }
        if(tmp3!=tmp4)
        {
          fprintf(stderr,"fi1 = %f\tfi2 = %f\n",tmp3,tmp4);  
        }
      }
      else
      { //fprintf(stderr,"lascou\n");
        int t;
        for(t=0;t<2;t++)
        {
          int s = ipow(-1,t);
          /* (TRUE) se eh vizinho da esquerda e (FALSE) para vizinho da direita */
          bool_t vizesq = s<0 ? TRUE : FALSE;
          if(vizesq)
          {
            No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
            No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
            assert((pf0==NULL)==(pf1==NULL));
            /* o vizinho eh refinado? */
            bool_t vizref = pf0==NULL ? FALSE : TRUE;
            if(vizref)
            {
              int j;
              double xml[d], xMl[d];
              for(j=0;j<d;j++)
              {
                xml[j] = j == dir ? xmin[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmin[j];
                xMl[j] = j == dir ? xmax[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmax[j];
              }
              /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
              {pac[ctr+s*corr].p} e some os fluxos superiores desses descendentes. */
              double somadosfluxos = 0.0;
              encontra_vizinhos_inferiores_esq_refinada( d, pac[ctr-corr].p, prof, dir, xml, xMl, xmin, &somadosfluxos);
              pac[ctr].p->flxinf[dir] = somadosfluxos;
            }
            else
            { 
              pac[ctr].p->flxinf[dir] = fluxo_roe(u1, u2, u3, u4, df, f);
            }
          }
          else /* direita */
          {
            No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
            No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
            assert((pf0==NULL)==(pf1==NULL));
            /* o vizinho eh refinado? */
            bool_t vizref = pf0==NULL ? FALSE : TRUE;
            if(vizref)
            {
              int j;
              double xml[d], xMl[d];
              for(j=0;j<d;j++)
              {
                xml[j] = j == dir ? xmin[j] + s*(xmax[j]-xmin[j]) : xmin[j];
                xMl[j] = j == dir ? xmax[j] + s*(xmax[j]-xmin[j]) : xmax[j];
              }
              /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
              {pac[ctr+s*corr].p} e some os fluxos inferiores desses descendentes. */
              double somadosfluxos = 0.0;
              encontra_vizinhos_superiores_dir_refinada( d, pac[ctr+corr].p, prof, dir, xml, xMl, xmax, &somadosfluxos);
              pac[ctr].p->flxsup[dir] = somadosfluxos;
            }
            else
            { 
              pac[ctr].p->flxsup[dir] = fluxo_roe(u2, u3, u4, u5, df, f);
            }
          }
        }
      }
      //return;
    }//fim for(dir)
  }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    int tam = ipow(ord,d);
    int ex =  prof % d;
    double xminf[d], xmaxf[d];
    VNo pacf[tam];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    divide_pacote_geral( d, pac,ex, pred, qual, pacf, op);
    calcula_fluxosT( d, prof+1, pacf, xminf, xmaxf, pred,f, lim,op);
  }
}

void calcula_fluxo_weno_v2(int d, int prof, VNo pac[], double xmin[], double xmax[], Fluxo f, Preditor pred, bool_t rd, bool_t op, int lim)
{
  
  /* ord==9 eh compativel com r==3?? */
  int ord = op ? 9 : 5;
  int r = rd ? 3 : 2;
  int ctr = ipow(ord,d)/2;
  if(pac[ctr].p == NULL) {return;}
  
  No *pf0 = pac[ctr].p->fil[0];
  No *pf1 = pac[ctr].p->fil[1];
  assert((pf0 == NULL) == (pf1 == NULL));

//  if((pf0 == NULL) && (prof < lim )){ return;} /* porque isso aqui mesmo ? */
  if( prof > lim ){ return;}  
  if((pf0 == NULL) /*&&  (prof == lim )*/)
  { 
    /* Reconstrucao em pac[ctr].p*/
    int dir;
    for(dir=0;dir<d;dir++)
    {
      int corr = ipow(ord,dir);
      double u[ord];
      int i, linf = ord/2, lsup = (ord/2+1);
      for(i=-linf;i<lsup;i++)
      {
        u[i+linf] = pac[ctr+i*corr].p == NULL ? pac[ctr+i*corr].fv : pac[ctr+i*corr].p->fval;
      }

      double varf = f(u[ord/2+corr])-f(u[ord/2-corr]);
        // na verdade, o que importa aqui eh o sinal da derivada (sdfu)
      int sdfu = varf <= 0.0 ? -1 : 1;

      double xmaxx = xmax[dir];
      double xminn = xmin[dir];

      if(r==2)
      {

        /* se todos os vizinhos sao folhas ou nulos */
        bool_t tvsfon = folhasounulos( d, ord, pac);
        if(tvsfon)
        { 
          //double Reconstrucao(int j, double x, int sdfu, double xmaxx, double xminn)
          int j = ord/2-1;
          double A = ReconstrucaoR2(j, xmax[dir], sdfu, xmaxx, xminn, u);
          double B = ReconstrucaoR2(j+1, xmax[dir], sdfu, xmaxx, xminn, u);
          pac[ctr].p->flxsup[dir] = H(A, B, f);

          A = ReconstrucaoR2(j-1, xmin[dir], sdfu, xmaxx, xminn, u);
          B = ReconstrucaoR2(j, xmin[dir], sdfu, xmaxx, xminn, u);
          pac[ctr].p->flxinf[dir] = H( A, B, f);
        }
        else
        { 
          int t;
          for(t=0;t<2;t++)
          {
            int s = ipow(-1,t);
            /* (TRUE) se eh vizinho da esquerda e (FALSE) para vizinho da direita */
            bool_t vizesq = s<0 ? TRUE : FALSE;
            if(vizesq)
            {
              No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
              No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
              assert((pf0==NULL)==(pf1==NULL));
              /* o vizinho eh refinado? */
              bool_t vizref = pf0==NULL ? FALSE : TRUE;
              if(vizref)
              {
                int j;
                double xml[d], xMl[d];
                for(j=0;j<d;j++)
                {
                  double hl = (xmax[j]-xmin[j]);
                  xml[j] = j == dir ? xmin[j] /*+ s**/ - hl : xmin[j];
                  xMl[j] = j == dir ? xmax[j] /*+ s**/ - hl : xmax[j];
                }
                /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
                {pac[ctr+s*corr].p} e some os fluxos superiores desses descendentes. */
                double somadosfluxos = 0.0;
                encontra_vizinhos_inferiores_esq_refinada( d, pac[ctr-corr].p, prof, dir, xml, xMl, xmin, &somadosfluxos);
                pac[ctr].p->flxinf[dir] = somadosfluxos;
              }
              else
              { 
                double A = ReconstrucaoR2(ord/2-1, xmin[dir], sdfu, xmaxx, xminn, u);
                double B = ReconstrucaoR2(ord/2, xmin[dir], sdfu, xmaxx, xminn, u);
                pac[ctr].p->flxinf[dir] = H( A, B, f);

              }
            }
            else /* direita */
            {
              No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
              No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
              assert((pf0==NULL)==(pf1==NULL));
              /* o vizinho eh refinado? */
              bool_t vizref = pf0==NULL ? FALSE : TRUE;
              if(vizref)
              {
                int j;
                double xml[d], xMl[d];
                for(j=0;j<d;j++)
                {
                  xml[j] = j == dir ? xmin[j] + s*(xmax[j]-xmin[j]) : xmin[j];
                  xMl[j] = j == dir ? xmax[j] + s*(xmax[j]-xmin[j]) : xmax[j];
                }
                /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
                {pac[ctr+s*corr].p} e some os fluxos inferiores desses descendentes. */
                double somadosfluxos = 0.0;
                encontra_vizinhos_superiores_dir_refinada( d, pac[ctr+corr].p, prof, dir, xml, xMl, xmax, &somadosfluxos);
                pac[ctr].p->flxsup[dir] = somadosfluxos;
              }
              else
              { 
                double A = ReconstrucaoR2(ord/2, xmax[dir], sdfu, xmaxx, xminn,u);
                double B = ReconstrucaoR2(ord/2+1, xmax[dir], sdfu, xmaxx, xminn,u);  
                pac[ctr].p->flxsup[dir] = H( A, B, f);
              }
            }
          }
        }
      }
      else //r==3
      {

        /* se todos os vizinhos sao folhas ou nulos */
        bool_t tvsfon = folhasounulos( d, ord, pac);
        if(tvsfon)
        { 
          double A = ReconstrucaoR3(ord/2, xmax[dir], sdfu, xmax[dir], xmin[dir],u);
          double B = ReconstrucaoR3(ord/2+1, xmax[dir], sdfu, xmax[dir], xmin[dir],u);
          pac[ctr].p->flxsup[dir] = H( A, B, f);
          A = ReconstrucaoR3(ord/2-1, xmin[dir], sdfu, xmax[dir], xmin[dir], u);
          B = ReconstrucaoR3(ord/2, xmin[dir], sdfu, xmax[dir], xmin[dir], u);
          pac[ctr].p->flxinf[dir] = H( A, B, f);
        }
        else
        { 
          int t;
          for(t=0;t<2;t++)
          {
            int s = ipow(-1,t);
            /* (TRUE) se eh vizinho da esquerda e (FALSE) para vizinho da direita */
            bool_t vizesq = s<0 ? TRUE : FALSE;
            if(vizesq)
            {
              No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
              No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
              assert((pf0==NULL)==(pf1==NULL));
              /* o vizinho eh refinado? */
              bool_t vizref = pf0==NULL ? FALSE : TRUE;
              if(vizref)
              {
                int j;
                double xml[d], xMl[d];
                for(j=0;j<d;j++)
                {
                  xml[j] = j == dir ? xmin[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmin[j];
                  xMl[j] = j == dir ? xmax[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmax[j];
                }
                /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
                {pac[ctr+s*corr].p} e some os fluxos superiores desses descendentes. */
                double somadosfluxos = 0.0;
                encontra_vizinhos_inferiores_esq_refinada( d, pac[ctr-corr].p, prof, dir, xml, xMl, xmin, &somadosfluxos);
                pac[ctr].p->flxinf[dir] = somadosfluxos;
              }
              else
              {
               double A = ReconstrucaoR3(ord/2-1, xmin[dir], sdfu, xmax[dir], xmin[dir], u);
               double B = ReconstrucaoR3(ord/2, xmin[dir], sdfu, xmax[dir], xmin[dir], u);
               pac[ctr].p->flxinf[dir] = H( A, B, f);
              }
            }
            else /* direita */
            {
              No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
              No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
              assert((pf0==NULL)==(pf1==NULL));
              /* o vizinho eh refinado? */
              bool_t vizref = pf0==NULL ? FALSE : TRUE;
              if(vizref)
              {
                int j;
                double xml[d], xMl[d];
                for(j=0;j<d;j++)
                {
                  xml[j] = j == dir ? xmin[j] + s*(xmax[j]-xmin[j]) : xmin[j];
                  xMl[j] = j == dir ? xmax[j] + s*(xmax[j]-xmin[j]) : xmax[j];
                }
                /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
                {pac[ctr+s*corr].p} e some os fluxos inferiores desses descendentes. */
                double somadosfluxos = 0.0;
                encontra_vizinhos_superiores_dir_refinada( d, pac[ctr+corr].p, prof, dir, xml, xMl, xmax, &somadosfluxos);
                pac[ctr].p->flxsup[dir] = somadosfluxos;
              }
              else
              { // tem que ser derivada em R(x)
                double A = ReconstrucaoR3(ord/2, xmax[dir], sdfu, xmax[dir], xmin[dir], u);
                double B =  ReconstrucaoR3(ord/2+1, xmax[dir], sdfu, xmax[dir], xmin[dir], u);
                pac[ctr].p->flxsup[dir] = H( A, B, f);
              }
            }
          }
        } 
      }
    }
  }
  int qual;
  for(qual=0;qual<2;qual++)
  {
    int tam = ipow(ord,d);
    int ex =  prof % d;
    double xminf[d], xmaxf[d];
    VNo pacf[tam];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    divide_pacote_geral( d, pac, ex, pred, qual, pacf, op);
    calcula_fluxo_weno_v2( d, prof+1, pacf, xminf, xmaxf, f, pred, rd, op, lim);
  }
}

void calcula_fluxos_upwind(int d, int prof, VNo pac[], double xmin[], double xmax[], Preditor pred, Fluxo f, int lim, bool_t op)
{
  int ord = 5;
  int tam = ipow(ord,d);
  int ctr = tam/2;

  if(pac[ctr].p == NULL)
  {
    return;
  }

  No *pf0 = /*pac[ctr].p == NULL ? NULL :*/ pac[ctr].p->fil[0];
  No *pf1 = /*pac[ctr].p == NULL ? NULL :*/ pac[ctr].p->fil[1];
  assert((pf0 == NULL) == (pf1 == NULL));

  if((pf0 == NULL) && (prof < lim )){ return;}
  if( prof > lim ){ return;}  
  if((pf0 == NULL) &&  (prof == lim ))
  { //fprintf(stderr,"prof = %d lim = %d\n", prof, lim);
    int dir;
    for(dir=0;dir<d;dir++)
    {
      int corr = ipow(ord,dir);
      //double u1 = pac[ctr-2*corr].p == NULL ? pac[ctr-2*corr].fv : pac[ctr - 2*corr].p->fval;
      double u2 = pac[ctr-1*corr].p == NULL ? pac[ctr-1*corr].fv : pac[ctr - 1*corr].p->fval;
      double u3 = pac[ctr+0*corr].p == NULL ? pac[ctr+0*corr].fv : pac[ctr + 0*corr].p->fval;
      double u4 = pac[ctr+1*corr].p == NULL ? pac[ctr+1*corr].fv : pac[ctr + 1*corr].p->fval;
      //double u5 = pac[ctr+2*corr].p == NULL ? pac[ctr+2*corr].fv : pac[ctr + 2*corr].p->fval;
      double df = (f(u4)-f(u3))/(u4-u3);
      /* se todos os vizinhos sao folhas ou nulos */
      bool_t tvsfon = folhasounulos( d, ord, pac);
      if(tvsfon)
      {
        //fprintf(stderr,"tudo folha\n");
        pac[ctr].p->flxsup[dir] = fluxo_upwind(u3, u4, df, f);
        pac[ctr].p->flxinf[dir] = fluxo_upwind(u2, u3, df, f);
      }
      else
      {
        int t;
        for(t=0;t<2;t++)
        {
          int s = ipow(-1,t);
          /* (TRUE) se eh vizinho da esquerda e (FALSE) para vizinho da direita */
          bool_t vizesq = s<0 ? TRUE : FALSE;
          if(vizesq)
          {
            No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
            No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
            assert((pf0==NULL)==(pf1==NULL));
            /* o vizinho eh refinado? */
            bool_t vizref = pf0==NULL ? FALSE : TRUE;
            if(vizref)
            {
              int j;
              double xml[d], xMl[d];
              for(j=0;j<d;j++)
              {
                xml[j] = j == dir ? xmin[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmin[j];
                xMl[j] = j == dir ? xmax[j] /*-*/ + s*(xmax[j]-xmin[j]) : xmax[j];
              }
              /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
              {pac[ctr+s*corr].p} e some os fluxos superiores desses descendentes. */
              double somadosfluxos = 0.0;
              encontra_vizinhos_inferiores_esq_refinada( d, pac[ctr-corr].p, prof, dir, xml, xMl, xmin, &somadosfluxos);
              pac[ctr].p->flxinf[dir] = somadosfluxos;
            }
            else
            { 
              pac[ctr].p->flxinf[dir] = fluxo_upwind(u2, u3, df, f);
            }
          }
          else /* direita */
          {
            No *pf0 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[0];
            No *pf1 = pac[ctr+s*corr].p == NULL ? NULL : pac[ctr+s*corr].p->fil[1];
            assert((pf0==NULL)==(pf1==NULL));
            /* o vizinho eh refinado? */
            bool_t vizref = pf0==NULL ? FALSE : TRUE;
            if(vizref)
            {
              int j;
              double xml[d], xMl[d];
              for(j=0;j<d;j++)
              {
                xml[j] = j == dir ? xmin[j] + s*(xmax[j]-xmin[j]) : xmin[j];
                xMl[j] = j == dir ? xmax[j] + s*(xmax[j]-xmin[j]) : xmax[j];
              }
              /* procure os descendentes de {pac[ctr+s*corr].p} que fazem fronteira com
              {pac[ctr+s*corr].p} e some os fluxos inferiores desses descendentes. */
              double somadosfluxos = 0.0;
              encontra_vizinhos_superiores_dir_refinada( d, pac[ctr+corr].p, prof, dir, xml, xMl, xmax, &somadosfluxos);
              pac[ctr].p->flxsup[dir] = somadosfluxos;
            }
            else
            { 
              pac[ctr].p->flxsup[dir] = fluxo_upwind(u3, u4, df, f);
            }
          }
        }
      }
    }//fim for(dir)
  }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    int tam = ipow(ord,d);
    int ex =  prof % d;
    double xminf[d], xmaxf[d];
    VNo pacf[tam];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    divide_pacote_geral( d, pac,ex, pred, qual, pacf,op);
    calcula_fluxos_upwind( d, prof+1, pacf, xminf, xmaxf, pred,f, lim,op);
  }
}


double ReconstrucaoR2(int j, double x, int sdfu, double xmaxx, double xminn, double u[])
{
  double dx = fabs(xmaxx - xminn); // h
  double pj, pj1, ISj, ISj1;
  double eps = 0.00001;
  double Alpha[2];
  double xj = 0.5*(xmaxx + xminn); // rever

/*
  pj  = u[j]+((u[j]-u[j-1])/h)*(x-xj-h); 
  pj1 = u[j]+((u[j+1]-u[j])/h)*(x-xj);
  Sj  = (u[j]-u[j-1])*(u[j]-u[j-1]);
  Sj1 = (u[j+1]-u[j])*(u[j+1]-u[j]);
*/
  pj  = u[j]  +  ((x-xj+dx) * (u[j]-u[j-1]) )/dx; 
  pj1 = u[j]  +  ( (x-xj)   * (u[j+1]-u[j]) )/dx;
  ISj  = (u[j]-u[j-1])*(u[j]-u[j-1]);
  ISj1 = (u[j+1]-u[j])*(u[j+1]-u[j]);

  if(sdfu>0)
  {
    Alpha[0] = 1.0/(   2.0*(eps+ISj)*(eps+ISj)   ); 
    Alpha[1] = 1.0/(   (eps+ISj1)*(eps+ISj1)     ); 
  }
  else
  {
    Alpha[0] = 1.0/(   (eps+ISj)*(eps+ISj)       ); 
    Alpha[1] = 1.0/(   2.0*(eps+ISj1)*(eps+ISj1) );
  }

  double S = Alpha[0] + Alpha[1];
fprintf(stderr,"S = %f\n",S);
  double R = (Alpha[0]/S)*pj + (Alpha[1]/S)*pj1; 
  double aaa = R;
  assert(aaa==R);
  return R;
}


double ReconstrucaoR3(int j, double x, int sdfu, double xmaxx, double xminn, double u[])
{
  double pj, pj1, pj2, Sj, Sj1, Sj2;
  double Alpha[3];
  double eps = 0.00001;
  //int dif = ord/2 - j;
  double xj = 0.5*(xmaxx + xminn); // + dif*h;

  double h = fabs(xmaxx - xminn); 

  pj  = ((u[j]-2.0*u[j-1]+u[j-2])/(2.0*h*h))*(x-xj-h)*(x-xj-h)+((u[j]-u[j-2])/(2.0*h))*(x-xj-h)+u[j-1]-(u[j]-2.0*u[j-1]+u[j-2])/(24.0);
  pj1 = ((u[j+1]-2.0*u[j]+u[j-1])/(2.0*h*h))*(x-xj)*(x-xj)+((u[j+1]-u[j-1])/(2.0*h))*(x-xj)+u[j]-(u[j+1]-2.0*u[j]+u[j-1])/(24.0);
  pj2 = ((u[j+2]-2.0*u[j+1]+u[j])/(2.0*h*h))*(x-(xj+h))*(x-(xj+h))+((u[j+2]-u[j])/(2.0*h))*(x-(xj+h))+u[j+1]-(u[j+2]-2.0*u[j+1]+u[j])/(24.0);
  Sj =  ((u[j-1]-u[j-2])*(u[j-1]-u[j-2])+(u[j]-u[j-1])*(u[j]-u[j-1]))/2.0 + (u[j]-2.0*u[j-1]+u[j-2])*(u[j]-2.0*u[j-1]+u[j-2]);
  Sj1 = ((u[j]-u[j-1])*(u[j]-u[j-1])+(u[j+1]-u[j])*(u[j+1]-u[j]))/2.0 + (u[j+1]-2.0*u[j]+u[j-1])*(u[j+1]-2.0*u[j]+u[j-1]);
  Sj2 = ((u[j+1]-u[j])*(u[j+1]-u[j])+(u[j+2]-u[j+1])*(u[j+2]-u[j+1]))/2.0 + (u[j+2]-2.0*u[j+1]+u[j])*(u[j+2]-2.0*u[j+1]+u[j]);

  if(sdfu>0)
  {
    Alpha[0] = 1.0/(12.0*(eps+Sj)*(eps+Sj)*(eps+Sj));
    Alpha[1] = 1.0/(2.0*(eps+Sj1)*(eps+Sj1)*(eps+Sj1));
    Alpha[2] = 1.0/(4.0*(eps+Sj2)*(eps+Sj2)*(eps+Sj2));
  }
  else
  {
    Alpha[0] = 1.0/(4.0*(eps+Sj)*(eps+Sj)*(eps+Sj));
    Alpha[1] = 1.0/(2.0*(eps+Sj1)*(eps+Sj1)*(eps+Sj1));
    Alpha[2] = 1.0/(12.0*(eps+Sj2)*(eps+Sj2)*(eps+Sj2));
  }

  double S = Alpha[0]+Alpha[1]+Alpha[2];
  double R = (Alpha[0]/S)*pj + (Alpha[1]/S)*pj1+ (Alpha[2]/S)*pj2;
  return R;
}


      // implementar aqui o fluxo H(.,.), que eh Lipschitz
      // monotona nao decrescente no primeiro argumento e
      // monotona nao crescente no segundo argumento 
      // Liu, Osher, Chan, Journal of Computational Physics 115, 200-212  1994.
double H(double a, double b, Fluxo f)
      {
        // usando de Godunov, mais facil de implementar
        double V; 
        int i;
        double fmax = -INF;
        double fmin = +INF;
        for(i=0;i<21;i++)
        {
          double linf = MIN(a,b);
          double lsup = MAX(a,b);
          double h = (lsup-linf)/20.0;
          double v = f(linf+i*h);
          double temp = fmin;
          //fmin = temp > v ? v : temp;
          fmin = MIN(temp, v);
          double tmp = fmax;
          //fmax = tmp < v ? v : tmp;
          fmax = MAX(tmp, v );
        }
        if(a<=b){ V = fmin;}
        else{ V = fmax;}

        return V;
      }

/*
      double H(double a, double b, Fluxo f)
      {
        // usando de Roe com entropy fix
        double V; 
        int i;
        double linf = MIN(a,b);
        double lsup = MAX(a,b);
        double beta = -INF;
        bool_t cresc = FALSE, decre = FALSE;
        double h = fabs( (lsup-linf)/20.0 );
        for(i=0;i<20;i++)
        {
          double df = (  f(linf+(i+1)*h) - f(linf+(i-1)*h)   )/(2.0*h);
          double v = fabs(df);
          double temp = beta;
          //beta = temp < v ? v : temp;
          beta = MAX(temp,v);
          if(df>=0.0){ cresc = TRUE; }
          else{ decre = TRUE; }
        }
        if(cresc && decre)
        { V = 0.5*( f(a) + f(b) -beta*(b-a) ); }
        else if(cresc && (!decre))
        { V = f(a); }
        else
        { V = f(b); }

        return V;
      }

*/