// $Id: Post.cpp,v 1.52 2004-02-20 17:58:00 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
// 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.
// 
// Please report all bugs and problems to <gmsh@geuz.org>.

#include "Gmsh.h"
#include "GmshUI.h"
#include "Numeric.h"
#include "Geo.h"
#include "Mesh.h"
#include "Draw.h"
#include "Views.h"
#include "Context.h"
#include "gl2ps.h"

extern Context_T CTX;

static double Raise[3][8];
static double RaiseFactor[3];

// Give Value from Index

double GiveValueFromIndex_Lin(double ValMin, double ValMax, int NbIso,
                              int Iso)
{
  if(NbIso == 1)
    return (ValMax + ValMin) / 2.;
  return ValMin + Iso * (ValMax - ValMin) / (NbIso - 1.);
}

double GiveValueFromIndex_Log(double ValMin, double ValMax, int NbIso,
                              int Iso)
{
  if(NbIso == 1)
    return (ValMax + ValMin) / 2.;
  if(ValMin <= 0.)
    return 0.;
  return pow(10.,
             log10(ValMin) + Iso * (log10(ValMax) - log10(ValMin)) / (NbIso -
                                                                      1.));
}

double GiveValueFromIndex_DoubleLog(double ValMin, double ValMax, int NbIso,
                                    int Iso)
{
  if(NbIso == 1)
    return (ValMax + ValMin) / 2.;
  if(ValMin <= 0.)
    return 0.;

  double Iso2 = Iso / 2.;
  double NbIso2 = NbIso / 2.;

  return pow(10.,
             log10(ValMin) + Iso2 * (log10(ValMax) -
                                     log10(ValMin)) / (NbIso2 - 1.));
}

// Give Index From Value

int GiveIndexFromValue_Lin(double ValMin, double ValMax, int NbIso,
                           double Val)
{
  if(ValMin == ValMax)
    return NbIso / 2;
  return (int)((Val - ValMin) * (NbIso - 1) / (ValMax - ValMin));
}

int GiveIndexFromValue_Log(double ValMin, double ValMax, int NbIso,
                           double Val)
{
  if(ValMin == ValMax)
    return NbIso / 2;
  if(ValMin <= 0.)
    return 0;
  return (int)((log10(Val) - log10(ValMin)) * (NbIso - 1) / (log10(ValMax) -
                                                             log10(ValMin)));
}

int GiveIndexFromValue_DoubleLog(double ValMin, double ValMax, int NbIso,
                                 double Val)
{
  if(ValMin == ValMax)
    return NbIso / 2;
  if(ValMin <= 0.)
    return 0;
  return (int)((log10(Val) - log10(ValMin)) * (NbIso - 1) / (log10(ValMax) -
                                                             log10(ValMin)));
}


// Color Palette

void Palette1(Post_View * v, int nbi, int i)
{       /* i in [0,nbi-1] */
  int index;

  index = (nbi == 1) ?
    v->CT.size / 2 : (int)(i / (double)(nbi - 1) * (v->CT.size - 1) + 0.5);

  glColor4ubv((GLubyte *) & v->CT.table[index]);
}

void Palette2(Post_View * v, double min, double max, double val)
{       /* val in [min,max] */
  int index;

  index = (min == max) ?
    v->CT.size / 2 :
    (int)((val - min) / (max - min) * (v->CT.size - 1) + 0.5);

  glColor4ubv((GLubyte *) & v->CT.table[index]);
}

void RaiseFill(int i, double Val, double ValMin, double Raise[3][8])
{
  int j;
  for(j = 0; j < 3; j++)
    Raise[j][i] = (Val - ValMin) * RaiseFactor[j];
}


// Compute node coordinates taking Offset and Explode into account

void Get_Coords(double Explode, double *Offset, int nbnod,
                double *x1, double *y1, double *z1,
                double *x2, double *y2, double *z2)
{
  int i;
  double xc = 0., yc = 0., zc = 0.;

  if(Explode == 1) {
    for(i = 0; i < nbnod; i++) {
      x2[i] = x1[i] + Offset[0];
      y2[i] = y1[i] + Offset[1];
      z2[i] = z1[i] + Offset[2];
    }
  }
  else {
    for(i = 0; i < nbnod; i++) {
      xc += x1[i];
      yc += y1[i];
      zc += z1[i];
    }
    xc /= (double)nbnod;
    yc /= (double)nbnod;
    zc /= (double)nbnod;
    for(i = 0; i < nbnod; i++) {
      x2[i] = xc + Explode * (x1[i] - xc) + Offset[0];
      y2[i] = yc + Explode * (y1[i] - yc) + Offset[1];
      z2[i] = zc + Explode * (z1[i] - zc) + Offset[2];
    }
  }
}

// Compare barycenters with viewpoint (eye)

static double storedEye[3] = { 0., 0., 0. };

int changedEye()
{
  double zeye = 100 * CTX.lc, tmp[3];
  tmp[0] = CTX.rot[0][2] * zeye;
  tmp[1] = CTX.rot[1][2] * zeye;
  tmp[2] = CTX.rot[2][2] * zeye;
  if(fabs(tmp[0] - storedEye[0]) > 1.e-3 ||
     fabs(tmp[1] - storedEye[1]) > 1.e-3 ||
     fabs(tmp[2] - storedEye[2]) > 1.e-3) {
    storedEye[0] = tmp[0];
    storedEye[1] = tmp[1];
    storedEye[2] = tmp[2];
    Msg(DEBUG, "New eye = (%g %g %g)", tmp[0], tmp[1], tmp[2]);
    return 1;
  }
  return 0;
}

// to be rigorous, we should take Raise into account

int compareEye(double *q, double *w, int nbnodes)
{
  double d, dq, dw, cgq[3] = { 0., 0., 0. }, cgw[3] =
  {
  0., 0., 0.};
  for(int i = 0; i < nbnodes; i++) {
    cgq[0] += q[i];
    cgq[1] += q[i + nbnodes];
    cgq[2] += q[i + 2 * nbnodes];
    cgw[0] += w[i];
    cgw[1] += w[i + nbnodes];
    cgw[2] += w[i + 2 * nbnodes];
  }
  prosca(storedEye, cgq, &dq);
  prosca(storedEye, cgw, &dw);
  d = dq - dw;
  if(d > 0)
    return 1;
  if(d < 0)
    return -1;
  return 0;
}

int compareEye3Nodes(const void *a, const void *b)
{
  return compareEye((double *)a, (double *)b, 3);
}

int compareEye4Nodes(const void *a, const void *b)
{
  return compareEye((double *)a, (double *)b, 4);
}

int compareEye5Nodes(const void *a, const void *b)
{
  return compareEye((double *)a, (double *)b, 5);
}

int compareEye6Nodes(const void *a, const void *b)
{
  return compareEye((double *)a, (double *)b, 6);
}

int compareEye8Nodes(const void *a, const void *b)
{
  return compareEye((double *)a, (double *)b, 8);
}

// Draw_Post

void Draw_ScalarList(Post_View * v, double ValMin, double ValMax,
                     double Raise[3][8], List_T * list, int nbelm, int nbnod,
                     int smoothnormals, void (*draw) (Post_View *, int,
                                                      double, double,
                                                      double[3][8], double *,
                                                      double *, double *,
                                                      double *))
{
  int i, nb;
  double X[8], Y[8], Z[8];

  if(nbelm && v->DrawScalars) {
    nb = List_Nbr(list) / nbelm;
    if(smoothnormals && v->Light && v->SmoothNormals && v->Changed &&
       v->IntervalsType != DRAW_POST_ISO) {
      v->reset_normals(); 
      // we might save some normal stuff by checking if the change
      // actually changed the "geometry"... Should put e.g. Change=2
      // if timestep changed, etc.
      Msg(DEBUG, "Preprocessing of normals in view %d", v->Num);
      for(i = 0; i < List_Nbr(list); i += nb) {
        Get_Coords(v->Explode, v->Offset, nbnod,
                   (double *)List_Pointer_Fast(list, i),
                   (double *)List_Pointer_Fast(list, i + nbnod),
                   (double *)List_Pointer_Fast(list, i + 2 * nbnod), X, Y, Z);
        draw(v, 1, ValMin, ValMax, Raise, X, Y, Z,
             (double *)List_Pointer_Fast(list, i + 3 * nbnod));
      }
    }
    for(i = 0; i < List_Nbr(list); i += nb) {
      Get_Coords(v->Explode, v->Offset, nbnod,
                 (double *)List_Pointer_Fast(list, i),
                 (double *)List_Pointer_Fast(list, i + nbnod),
                 (double *)List_Pointer_Fast(list, i + 2 * nbnod), X, Y, Z);
      draw(v, 0, ValMin, ValMax, Raise, X, Y, Z,
           (double *)List_Pointer_Fast(list, i + 3 * nbnod));
    }
  }
}

void Draw_VectorList(Post_View * v, double ValMin, double ValMax,
                     double Raise[3][8], List_T * list, int nbelm, int nbnod,
                     void (*draw) (Post_View *, double, double, double[3][8],
                                   double *, double *, double *, double *))
{
  int i, nb;
  double X[8], Y[8], Z[8];

  if(nbelm && v->DrawVectors) {
    nb = List_Nbr(list) / nbelm;
    for(i = 0; i < List_Nbr(list); i += nb) {
      Get_Coords(v->Explode, v->Offset, nbnod,
                 (double *)List_Pointer_Fast(list, i),
                 (double *)List_Pointer_Fast(list, i + nbnod),
                 (double *)List_Pointer_Fast(list, i + 2 * nbnod), X, Y, Z);
      draw(v, ValMin, ValMax, Raise, X, Y, Z,
           (double *)List_Pointer_Fast(list, i + 3 * nbnod));
    }
  }
}

void Draw_TensorList(Post_View * v, double ValMin, double ValMax,
                     double Raise[3][8], List_T * list, int nbelm, int nbnod,
                     void (*draw) (Post_View *, double, double, double[3][8],
                                   double *, double *, double *, double *))
{
  int i, nb;
  double X[8], Y[8], Z[8];

  if(nbelm && v->DrawTensors) {
    nb = List_Nbr(list) / nbelm;
    for(i = 0; i < List_Nbr(list); i += nb) {
      Get_Coords(v->Explode, v->Offset, nbnod,
                 (double *)List_Pointer_Fast(list, i),
                 (double *)List_Pointer_Fast(list, i + nbnod),
                 (double *)List_Pointer_Fast(list, i + 2 * nbnod), X, Y, Z);
      draw(v, ValMin, ValMax, Raise, X, Y, Z,
           (double *)List_Pointer_Fast(list, i + 3 * nbnod));
    }
  }
}

void Draw_Post(void)
{
  int iView, j, k, nb;
  double ValMin = 0., ValMax = 0., AbsMax;
  Post_View *v;

  if(!CTX.post.list)
    return;

  if(!CTX.post.draw) {  // draw only the bbox of the visible views
    for(iView = 0; iView < List_Nbr(CTX.post.list); iView++) {
      v = (Post_View *) List_Pointer(CTX.post.list, iView);
      if(v->Visible && v->Type == DRAW_POST_3D) {
        glColor4ubv((GLubyte *) & CTX.color.fg);
        glBegin(GL_LINE_LOOP);
        glVertex3d(v->BBox[0], v->BBox[2], v->BBox[4]);
        glVertex3d(v->BBox[1], v->BBox[2], v->BBox[4]);
        glVertex3d(v->BBox[1], v->BBox[3], v->BBox[4]);
        glVertex3d(v->BBox[0], v->BBox[3], v->BBox[4]);
        glEnd();
        glBegin(GL_LINE_LOOP);
        glVertex3d(v->BBox[0], v->BBox[2], v->BBox[5]);
        glVertex3d(v->BBox[1], v->BBox[2], v->BBox[5]);
        glVertex3d(v->BBox[1], v->BBox[3], v->BBox[5]);
        glVertex3d(v->BBox[0], v->BBox[3], v->BBox[5]);
        glEnd();
        glBegin(GL_LINES);
        glVertex3d(v->BBox[0], v->BBox[2], v->BBox[4]);
        glVertex3d(v->BBox[0], v->BBox[2], v->BBox[5]);
        glVertex3d(v->BBox[1], v->BBox[2], v->BBox[4]);
        glVertex3d(v->BBox[1], v->BBox[2], v->BBox[5]);
        glVertex3d(v->BBox[1], v->BBox[3], v->BBox[4]);
        glVertex3d(v->BBox[1], v->BBox[3], v->BBox[5]);
        glVertex3d(v->BBox[0], v->BBox[3], v->BBox[4]);
        glVertex3d(v->BBox[0], v->BBox[3], v->BBox[5]);
        glEnd();
      }
    }
    return;
  }

  for(iView = 0; iView < List_Nbr(CTX.post.list); iView++) {

    v = (Post_View *) List_Pointer(CTX.post.list, iView);

    if(v->Visible && !v->Dirty) {

      // Sort the data % eye for transparency. Hybrid views (e.g. tri
      // + qua) or multiple views will be sorted incorrectly... You
      // can use Plugin(DecomposeInSimplex) and/or View->Combine for
      // this.

      if(CTX.alpha && ColorTable_IsAlpha(&v->CT) && 
	 (changedEye() || v->Changed)) {
        Msg(DEBUG, "Sorting view %d", v->Num);

        if(v->DrawScalars) {

          if(v->IntervalsType != DRAW_POST_ISO) {

            if(v->NbST && v->DrawTriangles) {
              nb = List_Nbr(v->ST) / v->NbST;
              qsort(v->ST->array, v->NbST, nb * sizeof(double),
                    compareEye3Nodes);
              v->Changed = 1;
            }
            if(v->NbSQ && v->DrawQuadrangles) {
              nb = List_Nbr(v->SQ) / v->NbSQ;
              qsort(v->SQ->array, v->NbSQ, nb * sizeof(double),
                    compareEye4Nodes);
              v->Changed = 1;
            }

          }

          // The following is of course not rigorous: we should store
          // the triangles generated during the iso computation, and
          // sort these... But this is better than doing nothing. If
          // you want a rigorous sorting of the iso-surfaces, just use
          // Plugin(CutMap).

          if(v->NbSS && v->DrawTetrahedra) {
            nb = List_Nbr(v->SS) / v->NbSS;
            qsort(v->SS->array, v->NbSS, nb * sizeof(double),
                  compareEye4Nodes);
            v->Changed = 1;
          }
          if(v->NbSH && v->DrawHexahedra) {
            nb = List_Nbr(v->SH) / v->NbSH;
            qsort(v->SH->array, v->NbSH, nb * sizeof(double),
                  compareEye8Nodes);
            v->Changed = 1;
          }
          if(v->NbSI && v->DrawPrisms) {
            nb = List_Nbr(v->SI) / v->NbSI;
            qsort(v->SI->array, v->NbSI, nb * sizeof(double),
                  compareEye6Nodes);
            v->Changed = 1;
          }
          if(v->NbSY && v->DrawPyramids) {
            nb = List_Nbr(v->SY) / v->NbSY;
            qsort(v->SY->array, v->NbSY, nb * sizeof(double),
                  compareEye5Nodes);
            v->Changed = 1;
          }

        }

      }

      if(CTX.post.display_lists && !v->Changed && v->DisplayListNum > 0) {

        Msg(DEBUG, "Call display List %d", v->DisplayListNum);
        glCallList(v->DisplayListNum);

      }
      else {

        if(CTX.post.display_lists) {
          if(v->DisplayListNum > 0) {
            Msg(DEBUG, "Delete display List %d", v->DisplayListNum);
            glDeleteLists(v->DisplayListNum, 1);
          }
          else {
            v->DisplayListNum = glGenLists(1);
            Msg(DEBUG, "Gen display list -> %d", v->DisplayListNum);
          }

          if(v->DisplayListNum > 0) {
            Msg(DEBUG, "New display List %d", v->DisplayListNum);
            glNewList(v->DisplayListNum, GL_COMPILE_AND_EXECUTE);
          }
          else
            Msg(GERROR, "Unable to create display list");
        }

        glPointSize(v->PointSize);
        gl2psPointSize(v->PointSize * CTX.print.eps_point_size_factor);

        glLineWidth(v->LineWidth);
        gl2psLineWidth(v->LineWidth * CTX.print.eps_line_width_factor);

        if(v->ShowElement)
          glEnable(GL_POLYGON_OFFSET_FILL);

        // force this
        if(v->IntervalsType == DRAW_POST_CONTINUOUS) {
          glShadeModel(GL_SMOOTH);
          glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
        }
        else {  // there is a bug in CutTriangle2D!! See Iso.cpp
          glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
        }

        switch (v->RangeType) {
        case DRAW_POST_RANGE_DEFAULT:
          ValMin = v->Min;
          ValMax = v->Max;
          break;
        case DRAW_POST_RANGE_CUSTOM:
          ValMin = v->CustomMin;
          ValMax = v->CustomMax;
          break;
        case DRAW_POST_RANGE_PER_STEP:
          ValMin = v->TimeStepMin[v->TimeStep];
          ValMax = v->TimeStepMax[v->TimeStep];
          break;
        }

        switch (v->ScaleType) {
        case DRAW_POST_LINEAR:
          v->GIFV = GiveIndexFromValue_Lin;
          v->GVFI = GiveValueFromIndex_Lin;
          break;
        case DRAW_POST_LOGARITHMIC:
          v->GIFV = GiveIndexFromValue_Log;
          v->GVFI = GiveValueFromIndex_Log;
          break;
        case DRAW_POST_DOUBLELOGARITHMIC:
          v->GIFV = GiveIndexFromValue_DoubleLog;
          v->GVFI = GiveValueFromIndex_DoubleLog;
          break;
        }

        AbsMax = DMAX(fabs(ValMin), fabs(ValMax));
        AbsMax = (AbsMax == 0.) ? 1. : AbsMax;

        for(j = 0; j < 3; j++) {
          RaiseFactor[j] = v->Raise[j] / AbsMax;
          for(k = 0; k < 5; k++)
            Raise[j][k] = 0.;
        }

        // Points
        if(v->DrawPoints) {
          if(v->Type == DRAW_POST_3D)
            Draw_ScalarList(v, ValMin, ValMax, Raise, v->SP, v->NbSP, 1, 0,
                            Draw_ScalarPoint);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VP, v->NbVP, 1,
                          Draw_VectorPoint);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TP, v->NbTP, 1,
                          Draw_TensorPoint);
        }

        // Lines
        if(v->DrawLines) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SL, v->NbSL, 2, 0,
                          Draw_ScalarLine);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VL, v->NbVL, 2,
                          Draw_VectorLine);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TL, v->NbTL, 2,
                          Draw_TensorLine);
        }

        // Triangles
        if(v->DrawTriangles) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->ST, v->NbST, 3, 1,
                          Draw_ScalarTriangle);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VT, v->NbVT, 3,
                          Draw_VectorTriangle);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TT, v->NbTT, 3,
                          Draw_TensorTriangle);
        }

        // Quadrangles
        if(v->DrawQuadrangles) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SQ, v->NbSQ, 4, 1,
                          Draw_ScalarQuadrangle);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VQ, v->NbVQ, 4,
                          Draw_VectorQuadrangle);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TQ, v->NbTQ, 4,
                          Draw_TensorQuadrangle);
        }

        // Tetrahedra
        if(v->DrawTetrahedra) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SS, v->NbSS, 4, 1,
                          Draw_ScalarTetrahedron);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VS, v->NbVS, 4,
                          Draw_VectorTetrahedron);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TS, v->NbTS, 4,
                          Draw_TensorTetrahedron);
        }

        // Hexahedra
        if(v->DrawHexahedra) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SH, v->NbSH, 8, 1,
                          Draw_ScalarHexahedron);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VH, v->NbVH, 8,
                          Draw_VectorHexahedron);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TH, v->NbTH, 8,
                          Draw_TensorHexahedron);
        }

        // Prisms
        if(v->DrawPrisms) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SI, v->NbSI, 6, 1,
                          Draw_ScalarPrism);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VI, v->NbVI, 6,
                          Draw_VectorPrism);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TI, v->NbTI, 6,
                          Draw_TensorPrism);
        }

        // Pyramids
        if(v->DrawPyramids) {
          Draw_ScalarList(v, ValMin, ValMax, Raise, v->SY, v->NbSY, 5, 1,
                          Draw_ScalarPyramid);
          Draw_VectorList(v, ValMin, ValMax, Raise, v->VY, v->NbVY, 5,
                          Draw_VectorPyramid);
          Draw_TensorList(v, ValMin, ValMax, Raise, v->TY, v->NbTY, 5,
                          Draw_TensorPyramid);
        }

        // Strings
        if(v->DrawStrings) {
          glColor4ubv((GLubyte *) & CTX.color.text);
          Draw_Text2D3D(3, v->TimeStep, v->NbT3, v->T3D, v->T3C);
        }

        if(CTX.post.display_lists)
          glEndList();

        v->Changed = 0;

        if(v->ShowElement || v->VectorType == DRAW_POST_DISPLACEMENT)
          glDisable(GL_POLYGON_OFFSET_FILL);

      }

    }

  }
}