/Post/PViewVertexArrays.cpp

Large files files are truncated, but you can click here to view the full file

// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#include <string.h>
#include "GmshMessage.h"
#include "GmshDefines.h"
#include "onelab.h"
#include "Iso.h"
#include "MEdge.h"
#include "MFace.h"
#include "PView.h"
#include "PViewOptions.h"
#include "PViewData.h"
#include "PViewDataRemote.h"
#include "Numeric.h"
#include "VertexArray.h"
#include "SmoothData.h"
#include "Context.h"
#include "OpenFile.h"
#include "mathEvaluator.h"
#include "Options.h"
#include "StringUtils.h"

static void saturate(int nb, double **val, double vmin, double vmax,
                     int i0=0, int i1=1, int i2=2, int i3=3,
                     int i4=4, int i5=5, int i6=6, int i7=7)
{
  int id[8] = {i0, i1, i2, i3, i4, i5, i6, i7};
  for(int i = 0; i < nb; i++){
    if(val[id[i]][0] > vmax)      val[id[i]][0] = vmax;
    else if(val[id[i]][0] < vmin) val[id[i]][0] = vmin;
  }
}

static double saturateVector(double *val, int numComp2, double *val2,
                             double min, double max)
{
  double v = ComputeScalarRep(numComp2, val2);  // v >= 0
  if (v < min && v > 1e-15) {
    double f = min / v;
    for (int iComp = 0; iComp < numComp2; ++iComp)
      val2[iComp] *= f;
    val[0] *= f;
    val[1] *= f;
    val[2] *= f;
    return min;
  }
  if (v > max && v > 1e-15) {
    double f = max / v;
    for (int iComp = 0; iComp < numComp2; ++iComp)
      val2[iComp] *= f;
    val[0] *= f;
    val[1] *= f;
    val[2] *= f;
    return max;
  }
  return v;
}

static SVector3 normal3(double **xyz, int i0=0, int i1=1, int i2=2)
{
  SVector3 t1(xyz[i1][0] - xyz[i0][0],
              xyz[i1][1] - xyz[i0][1],
              xyz[i1][2] - xyz[i0][2]);
  SVector3 t2(xyz[i2][0] - xyz[i0][0],
              xyz[i2][1] - xyz[i0][1],
              xyz[i2][2] - xyz[i0][2]);
  SVector3 n = crossprod(t1, t2);
  n.normalize();
  return n;
}

static SVector3 getPointNormal(PView *p, double v)
{
  PViewOptions *opt = p->getOptions();
  SVector3 n(0., 0., 0.);
  if(opt->pointType > 0){
    // when we draw spheres, we use the normalized value (between 0
    // and 1) stored in the first component of the normal to modulate
    // the radius
    double d = opt->tmpMax - opt->tmpMin;
    n[0] = (v - opt->tmpMin) / (d ? d : 1.);
  }
  return n;
}

static void getLineNormal(PView *p, double x[2], double y[2], double z[2],
                          double *v, SVector3 n[2], bool computeNormal)
{
  PViewOptions *opt = p->getOptions();

  if(opt->lineType > 0){
    if(v){
      // when we draw tapered cylinders, we use the normalized values
      // (between 0 and 1) stored in the first component of the
      // normals to modulate the width
      double d = opt->tmpMax - opt->tmpMin;
      n[0][0] = (v[0] - opt->tmpMin) / (d ? d : 1.);
      n[1][0] = (v[1] - opt->tmpMin) / (d ? d : 1.);
    }
    else{
      // when we don't have values we use maximum width cylinders
      n[0][0] = n[1][0] = 1.;
    }
  }
  else if(computeNormal){
    SBoundingBox3d bb = p->getData()->getBoundingBox();
    if(bb.min().z() == bb.max().z())
      n[0] = n[1] = SVector3(0., 0., 1.);
    else if(bb.min().y() == bb.max().y())
      n[0] = n[1] = SVector3(0., 1., 0.);
    else if(bb.min().x() == bb.max().x())
      n[0] = n[1] = SVector3(1., 0., 0.);
    else{
      // we don't have any info about the normal, just pick one
      SVector3 t(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
      SVector3 ex(0., 0., 0.);
      if(t[0] == 0.)
        ex[0] = 1.;
      else if(t[1] == 0.)
        ex[1] = 1.;
      else
        ex[2] = 1.;
      n[0] = crossprod(t, ex);
      n[0].normalize();
      n[1] = n[0];
    }
  }
}

static bool getExternalValues(PView *p, int index, int ient, int iele,
                              int numNodes, int numComp, double **val,
                              int &numComp2, double **val2)
{
  PViewOptions *opt = p->getOptions();

  // use self by default
  numComp2 = numComp;
  for(int i = 0; i < numNodes; i++)
    for(int j = 0; j < numComp; j++)
      val2[i][j] = val[i][j];
  opt->externalMin = opt->tmpMin;
  opt->externalMax = opt->tmpMax;

  if(index < 0 || index >= (int)PView::list.size()) return false;

  PView *p2 = PView::list[index];
  PViewData *data2 = p2->getData(true); // use adaptive data if available

  if(iele >= data2->getNumElements(opt->timeStep, ient)) return false;

  if(!data2->skipElement(opt->timeStep, ient, iele) &&
     data2->getNumNodes(opt->timeStep, ient, iele) == numNodes){
    numComp2 = data2->getNumComponents(opt->timeStep, ient, iele);
    for(int i = 0; i < numNodes; i++)
      for(int j = 0; j < numComp2; j++)
        data2->getValue(opt->timeStep, ient, iele, i, j, val2[i][j]);
    if(opt->rangeType == PViewOptions::Custom){
      opt->externalMin = opt->customMin;
      opt->externalMax = opt->customMax;
    }
    else if(opt->rangeType == PViewOptions::PerTimeStep){
      opt->externalMin = data2->getMin(opt->timeStep);
      opt->externalMax = data2->getMax(opt->timeStep);
    }
    else{
      opt->externalMin = data2->getMin();
      opt->externalMax = data2->getMax();
    }
    return true;
  }
  return false;
}

static void applyGeneralRaise(PView *p, int numNodes, int numComp,
                              double **vals, double **xyz)
{
  PViewOptions *opt = p->getOptions();
  if(!opt->genRaiseEvaluator) return;

  std::vector<double> values(12, 0.), res(3);
  for(int k = 0; k < numNodes; k++) {
    for(int i = 0; i < 3; i++) values[i] = xyz[k][i];
    for(int i = 0; i < std::min(numComp, 9); i++) values[3 + i] = vals[k][i];
    if(opt->genRaiseEvaluator->eval(values, res))
      for(int i = 0; i < 3; i++)
        xyz[k][i] += opt->genRaiseFactor * res[i];
  }
}

void changeCoordinates(PView *p, int ient, int iele,
                       int numNodes, int type, int numComp,
                       double **xyz, double **val)
{
  PViewOptions *opt = p->getOptions();

  if(opt->explode != 1.) {
    double barycenter[3] = {0., 0., 0.};
    for(int i = 0; i < numNodes; i++)
      for(int j = 0; j < 3; j++)
        barycenter[j] += xyz[i][j];
    for(int j = 0; j < 3; j++)
      barycenter[j] /= (double)numNodes;
    for(int i = 0; i < numNodes; i++)
      for(int j = 0; j < 3; j++)
        xyz[i][j] = barycenter[j] + opt->explode * (xyz[i][j] - barycenter[j]);
  }

  if(opt->transform[0][0] != 1. || opt->transform[0][1] != 0. ||
     opt->transform[0][2] != 0. || opt->transform[1][0] != 0. ||
     opt->transform[1][1] != 1. || opt->transform[1][2] != 0. ||
     opt->transform[2][0] != 0. || opt->transform[2][1] != 0. ||
     opt->transform[2][2] != 1.){
    for(int i = 0; i < numNodes; i++) {
      double old[3] = {xyz[i][0], xyz[i][1], xyz[i][2]};
      for(int j = 0; j < 3; j++){
        xyz[i][j] = 0.;
        for(int k = 0; k < 3; k++)
          xyz[i][j] += opt->transform[j][k] * old[k];
      }
    }
  }

  if(opt->offset[0] || opt->offset[1] || opt->offset[2]){
    for(int i = 0; i < numNodes; i++)
      for(int j = 0; j < 3; j++)
        xyz[i][j] += opt->offset[j];
  }

  if(opt->raise[0] || opt->raise[1] || opt->raise[2]){
    for(int i = 0; i < numNodes; i++){
      double v = ComputeScalarRep(numComp, val[i]);
      for(int j = 0; j < 3; j++)
        xyz[i][j] += opt->raise[j] * v;
    }
  }

  if(opt->normalRaise && (type == TYPE_LIN || type == TYPE_TRI || type ==TYPE_QUA)){
    SVector3 n;
    if(type == TYPE_LIN){
      // assumes lines in z=const plane, and raises in that plane
      double x[2] = {xyz[0][0], xyz[1][0]};
      double y[2] = {xyz[0][1], xyz[1][1]};
      double z[2] = {xyz[0][2], xyz[1][2]};
      SVector3 p(0, 0, 1.);
      SVector3 t(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
      n = crossprod(t, p);
      n.normalize();
    }
    else
      n = normal3(xyz);
    for(int i = 0; i < numNodes; i++){
      double v = ComputeScalarRep(numComp, val[i]);
      for(int j = 0; j < 3; j++)
        xyz[i][j] += n[j] * opt->normalRaise * v;
    }
  }

  if(numComp == 3 && opt->vectorType == PViewOptions::Displacement){
    for(int i = 0; i < numNodes; i++){
      for(int j = 0; j < 3; j++)
        xyz[i][j] += opt->displacementFactor * val[i][j];
    }
  }

  if(opt->useGenRaise){
    int numComp2;
    double **val2 = new double*[numNodes];
    for(int i = 0; i < numNodes; i++)
      val2[i] = new double[9];
    getExternalValues(p, opt->viewIndexForGenRaise, ient, iele, numNodes,
                      numComp, val, numComp2, val2);
    applyGeneralRaise(p, numNodes, numComp2, val2, xyz);
    for(int i = 0; i < numNodes; i++)
      delete [] val2[i];
    delete [] val2;
  }
}

static double evalClipPlane(int clip, double x, double y, double z)
{
  return CTX::instance()->clipPlane[clip][0] * x +
    CTX::instance()->clipPlane[clip][1] * y +
    CTX::instance()->clipPlane[clip][2] * z +
    CTX::instance()->clipPlane[clip][3];
}

static double intersectClipPlane(int clip, int numNodes, double **xyz)
{
  double val = evalClipPlane(clip, xyz[0][0], xyz[0][1], xyz[0][2]);
  for(int i = 1; i < numNodes; i++){
    if(val * evalClipPlane(clip, xyz[i][0], xyz[i][1], xyz[i][2]) <= 0)
      return 0.; // the element intersects the cut plane
  }
  return val;
}

bool isElementVisible(PViewOptions *opt, int dim, int numNodes,
                      double **xyz)
{
  if(!CTX::instance()->clipWholeElements) return true;
  bool hidden = false;
  for(int clip = 0; clip < 6; clip++){
    if(opt->clip & (1 << clip)){
      if(dim < 3 && CTX::instance()->clipOnlyVolume){
      }
      else{
        double d = intersectClipPlane(clip, numNodes, xyz);
        if(dim == 3 && CTX::instance()->clipOnlyDrawIntersectingVolume && d){
          hidden = true;
          break;
        }
        else if(d < 0){
          hidden = true;
          break;
        }
      }
    }
  }
  return !hidden;
}

static void addOutlinePoint(PView *p, double **xyz, unsigned int color,
                            bool pre, int i0=0)
{
  if(pre) return;
  SVector3 n = getPointNormal(p, 1.);
  p->va_points->add(&xyz[i0][0], &xyz[i0][1], &xyz[i0][2], &n, &color, 0, true);
}

static void addScalarPoint(PView *p, double **xyz,
                           double **val, bool pre, int i0=0,
                           bool unique=false)
{
  if(pre) return;

  PViewOptions *opt = p->getOptions();

  double vmin = opt->tmpMin, vmax = opt->tmpMax;
  if(opt->saturateValues) saturate(1, val, vmin, vmax, i0);

  if(val[i0][0] >= vmin && val[i0][0] <= vmax){
    unsigned int col = opt->getColor(val[i0][0], vmin, vmax, false,
                                     (opt->intervalsType == PViewOptions::Discrete) ?
                                     opt->nbIso : -1);
    SVector3 n = getPointNormal(p, val[i0][0]);
    p->va_points->add(&xyz[i0][0], &xyz[i0][1], &xyz[i0][2], &n, &col, 0, unique);
  }
}

static void addOutlineLine(PView *p, double **xyz, unsigned int color,
                           bool pre, int i0=0, int i1=1)
{
  if(pre) return;

  const int in[2] = {i0, i1};
  unsigned int col[2];
  double x[2], y[2], z[2];
  for(int i = 0; i < 2; i++){
    x[i] = xyz[in[i]][0]; y[i] = xyz[in[i]][1]; z[i] = xyz[in[i]][2];
    col[i] = color;
  }
  SVector3 n[2];
  getLineNormal(p, x, y, z, 0, n, true);
  p->va_lines->add(x, y, z, n, col, 0, true);
}

static void addScalarLine(PView *p, double **xyz,
                          double **val, bool pre, int i0=0, int i1=1,
                          bool unique=false)
{
  if(pre) return;

  PViewOptions *opt = p->getOptions();

  if(opt->boundary > 0){
    opt->boundary--;
    addScalarPoint(p, xyz, val, pre, i0, true);
    addScalarPoint(p, xyz, val, pre, i1, true);
    opt->boundary++;
    return;
  }

  double vmin = opt->tmpMin, vmax = opt->tmpMax;
  if(opt->saturateValues) saturate(2, val, vmin, vmax, i0, i1);

  double x[2] = {xyz[i0][0], xyz[i1][0]};
  double y[2] = {xyz[i0][1], xyz[i1][1]};
  double z[2] = {xyz[i0][2], xyz[i1][2]};
  double v[2] = {val[i0][0], val[i1][0]};

  if(opt->intervalsType == PViewOptions::Continuous){
    SVector3 n[2];
    getLineNormal(p, x, y, z, v, n, true);
    if(val[i0][0] >= vmin && val[i0][0] <= vmax &&
       val[i1][0] >= vmin && val[i1][0] <= vmax){
      unsigned int col[2];
      for(int i = 0; i < 2; i++)
        col[i] = opt->getColor(v[i], vmin, vmax);
      p->va_lines->add(x, y, z, n, col, 0, unique);
    }
    else{
      double x2[2], y2[2], z2[2], v2[2];
      int nb = CutLine(x, y, z, v, vmin, vmax, x2, y2, z2, v2);
      if(nb == 2){
        unsigned int col[2];
        for(int i = 0; i < 2; i++)
          col[i] = opt->getColor(v2[i], vmin, vmax);
        p->va_lines->add(x2, y2, z2, n, col, 0, unique);
      }
    }
  }

  if(opt->intervalsType == PViewOptions::Discrete){
    for(int k = 0; k < opt->nbIso; k++){
      if(vmin == vmax) k = opt->nbIso / 2;
      double min = opt->getScaleValue(k, opt->nbIso + 1, vmin, vmax);
      double max = opt->getScaleValue(k + 1, opt->nbIso + 1, vmin, vmax);
      double x2[2], y2[2], z2[2], v2[2];
      int nb = CutLine(x, y, z, v, min, max, x2, y2, z2, v2);
      if(nb == 2){
        unsigned int color = opt->getColor(k, opt->nbIso);
        unsigned int col[2] = {color, color};
        SVector3 n[2];
        getLineNormal(p, x2, y2, z2, v2, n, true);
        p->va_lines->add(x2, y2, z2, n, col, 0, unique);
      }
      if(vmin == vmax) break;
    }
  }

  if(opt->intervalsType == PViewOptions::Iso){
    for(int k = 0; k < opt->nbIso; k++) {
      if(vmin == vmax) k = opt->nbIso / 2;
      double iso = opt->getScaleValue(k, opt->nbIso, vmin, vmax);
      double x2[1], y2[1], z2[1];
      int nb = IsoLine(x, y, z, v, iso, x2, y2, z2);
      if(nb == 1){
        unsigned int color = opt->getColor(k, opt->nbIso);
        SVector3 n = getPointNormal(p, iso);
        p->va_points->add(x2, y2, z2, &n, &color, 0, unique);
      }
      if(vmin == vmax) break;
    }
  }
}

static void addOutlineTriangle(PView *p, double **xyz,
                               unsigned int color, bool pre, int i0=0,
                               int i1=1, int i2=2)
{
  PViewOptions *opt = p->getOptions();

  const int il[3][2] = {{i0, i1}, {i1, i2}, {i2, i0}};

  SVector3 nfac = normal3(xyz, i0, i1, i2);

  for(int i = 0; i < 3; i++){
    double x[2] = {xyz[il[i][0]][0], xyz[il[i][1]][0]};
    double y[2] = {xyz[il[i][0]][1], xyz[il[i][1]][1]};
    double z[2] = {xyz[il[i][0]][2], xyz[il[i][1]][2]};
    SVector3 n[2] = {nfac, nfac};
    unsigned int col[2] = {color, color};
    if(opt->smoothNormals){
      for(int j = 0; j < 2; j++){
        if(pre) p->normals->add(x[j], y[j], z[j], n[j][0], n[j][1], n[j][2]);
        else p->normals->get(x[j], y[j], z[j], n[j][0], n[j][1], n[j][2]);
      }
    }
    getLineNormal(p, x, y, z, 0, n, false);
    if(!pre) p->va_lines->add(x, y, z, n, col, 0, true);
  }
}

static void addScalarTriangle(PView *p, double **xyz, double **val,
                              bool pre, int i0=0, int i1=1, int i2=2,
                              bool unique=false, bool skin=false)
{
  PViewOptions *opt = p->getOptions();

  const int il[3][2] = {{i0, i1}, {i1, i2}, {i2, i0}};

  if(opt->boundary > 0){
    opt->boundary--;
    for(int i = 0; i < 3; i++)
      addScalarLine(p, xyz, val, pre, il[i][0], il[i][1], true);
    opt->boundary++;
    return;
  }

  double vmin = opt->tmpMin, vmax = opt->tmpMax;
  if(opt->saturateValues) saturate(3, val, vmin, vmax, i0, i1, i2);

  double x[3] = {xyz[i0][0], xyz[i1][0], xyz[i2][0]};
  double y[3] = {xyz[i0][1], xyz[i1][1], xyz[i2][1]};
  double z[3] = {xyz[i0][2], xyz[i1][2], xyz[i2][2]};
  double v[3] = {val[i0][0], val[i1][0], val[i2][0]};

  SVector3 nfac = normal3(xyz, i0, i1, i2);

  if(opt->intervalsType == PViewOptions::Continuous){
    if(val[i0][0] >= vmin && val[i0][0] <= vmax &&
       val[i1][0] >= vmin && val[i1][0] <= vmax &&
       val[i2][0] >= vmin && val[i2][0] <= vmax){
      SVector3 n[3] = {nfac, nfac, nfac};
      unsigned int col[3];
      for(int i = 0; i < 3; i++){
        if(opt->smoothNormals){
          if(pre) p->normals->add(x[i], y[i], z[i], n[i][0], n[i][1], n[i][2]);
          else p->normals->get(x[i], y[i], z[i], n[i][0], n[i][1], n[i][2]);
        }
        col[i] = opt->getColor(v[i], vmin, vmax);
      }
      if(!pre) p->va_triangles->add(x, y, z, n, col, 0, unique, skin);
    }
    else{
      double x2[10], y2[10], z2[10], v2[10];
      int nb = CutTriangle(x, y, z, v, vmin, vmax, x2, y2, z2, v2);
      if(nb >= 3){
        for(int j = 2; j < nb; j++){
          double x3[3] = {x2[0], x2[j - 1], x2[j]};
          double y3[3] = {y2[0], y2[j - 1], y2[j]};
          double z3[3] = {z2[0], z2[j - 1], z2[j]};
          double v3[3] = {v2[0], v2[j - 1], v2[j]};
          SVector3 n[3] = {nfac, nfac, nfac};
          unsigned int col[3];
          for(int i = 0; i < 3; i++){
            if(opt->smoothNormals){
              if(pre) p->normals->add(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
              else p->normals->get(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
            }
            col[i] = opt->getColor(v3[i], vmin, vmax);
          }
          if(!pre) p->va_triangles->add(x3, y3, z3, n, col, 0, unique, skin);
        }
      }
    }
  }

  if(opt->intervalsType == PViewOptions::Discrete){
    for(int k = 0; k < opt->nbIso; k++){
      if(vmin == vmax) k = opt->nbIso / 2;
      double min = opt->getScaleValue(k, opt->nbIso + 1, vmin, vmax);
      double max = opt->getScaleValue(k + 1, opt->nbIso + 1, vmin, vmax);
      double x2[10], y2[10], z2[10], v2[10];
      int nb = CutTriangle(x, y, z, v, min, max, x2, y2, z2, v2);
      if(nb >= 3){
        unsigned int color = opt->getColor(k, opt->nbIso);
        unsigned int col[3] = {color, color, color};
        for(int j = 2; j < nb; j++){
          double x3[3] = {x2[0], x2[j - 1], x2[j]};
          double y3[3] = {y2[0], y2[j - 1], y2[j]};
          double z3[3] = {z2[0], z2[j - 1], z2[j]};
          SVector3 n[3] = {nfac, nfac, nfac};
          if(opt->smoothNormals){
            for(int i = 0; i < 3; i++){
              if(pre) p->normals->add(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
              else p->normals->get(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
            }
          }
          if(!pre) p->va_triangles->add(x3, y3, z3, n, col, 0, unique, skin);
        }
      }
      if(vmin == vmax) break;
    }
  }

  if(opt->intervalsType == PViewOptions::Iso){
    for(int k = 0; k < opt->nbIso; k++) {
      if(vmin == vmax) k = opt->nbIso / 2;
      double iso = opt->getScaleValue(k, opt->nbIso, vmin, vmax);
      double x2[3], y2[3], z2[3];
      int nb = IsoTriangle(x, y, z, v, iso, x2, y2, z2);
      if(nb == 2){
        unsigned int color = opt->getColor(k, opt->nbIso);
        unsigned int col[2] = {color, color};
        SVector3 n[2] = {nfac, nfac};
        if(opt->smoothNormals){
          for(int i = 0; i < 2; i++){
            if(pre) p->normals->add(x2[i], y2[i], z2[i], n[i][0], n[i][1], n[i][2]);
            else p->normals->get(x2[i], y2[i], z2[i], n[i][0], n[i][1], n[i][2]);
          }
        }
        double v[2] = {iso, iso};
        getLineNormal(p, x, y, z, v, n, false);
        if(!pre) p->va_lines->add(x2, y2, z2, n, col, 0, unique);
      }
      if(vmin == vmax) break;
    }
  }
}

static void addOutlineQuadrangle(PView *p, double **xyz,
                                 unsigned int color, bool pre, int i0=0, int i1=1,
                                 int i2=2, int i3=3)
{
  PViewOptions *opt = p->getOptions();

  const int il[4][2] = {{i0, i1}, {i1, i2}, {i2, i3}, {i3, i0}};

  SVector3 nfac = normal3(xyz, i0, i1, i2);

  for(int i = 0; i < 4; i++){
    double x[2] = {xyz[il[i][0]][0], xyz[il[i][1]][0]};
    double y[2] = {xyz[il[i][0]][1], xyz[il[i][1]][1]};
    double z[2] = {xyz[il[i][0]][2], xyz[il[i][1]][2]};
    SVector3 n[2] = {nfac, nfac};
    unsigned int col[2] = {color, color};
    if(opt->smoothNormals){
      for(int j = 0; j < 2; j++){
        if(pre) p->normals->add(x[j], y[j], z[j], n[j][0], n[j][1], n[j][2]);
        else p->normals->get(x[j], y[j], z[j], n[j][0], n[j][1], n[j][2]);
      }
    }
    getLineNormal(p, x, y, z, 0, n, false);
    if(!pre) p->va_lines->add(x, y, z, n, col, 0, true);
  }
}

static void addScalarQuadrangle(PView *p, double **xyz,
                                double **val, bool pre, int i0=0,
                                int i1=1, int i2=2, int i3=3, bool unique=false)
{
  PViewOptions *opt = p->getOptions();

  const int il[4][2] = {{i0, i1}, {i1, i2}, {i2, i3}, {i3, i0}};
  const int it[2][3] = {{i0, i1, i2}, {i0, i2, i3}};

  if(opt->boundary > 0){
    opt->boundary--;
    for(int i = 0; i < 4; i++)
      addScalarLine(p, xyz, val, pre, il[i][0], il[i][1], true);
    opt->boundary++;
    return;
  }

  for(int i = 0; i < 2; i++)
    addScalarTriangle(p, xyz, val, pre, it[i][0], it[i][1], it[i][2], unique);
}

static void addOutlinePolygon(PView *p, double **xyz,
                              unsigned int color, bool pre, int numNodes)
{
  for(int i = 0; i < numNodes / 3; i++)
    addOutlineTriangle(p, xyz, color, pre, 3*i, 3*i + 1, 3*i + 2);
}

static void addScalarPolygon(PView *p, double **xyz,
                             double **val, bool pre, int numNodes)
{
  PViewOptions *opt = p->getOptions();

  if(opt->boundary > 0){
    const int il[3][2] = {{0, 1}, {1, 2}, {2, 0}};
    std::map<MEdge, int, Less_Edge> edges;
    std::vector<MVertex *> verts;
    for(int i = 0; i < numNodes; i++)
      verts.push_back(new MVertex(xyz[i][0], xyz[i][1], xyz[i][2]));
    for(int i = 0; i < numNodes / 3; i++){
      for(int j = 0; j < 3; j++) {
        MEdge ed(verts[3*i+il[j][0]], verts[3*i+il[j][1]]);
        std::map<MEdge, int, Less_Edge>::iterator ite;
        for(ite = edges.begin(); ite != edges.end(); ite++)
          if((*ite).first == ed) break;
        if(ite == edges.end())
          edges[ed] = 100*i+j;
        else edges.erase(ite);
      }
    }

    opt->boundary--;
    for(std::map<MEdge, int, Less_Edge>::iterator ite = edges.begin();
        ite != edges.end(); ite++){
      int i = (int) (*ite).second / 100;
      int j = (*ite).second % 100;
      addScalarLine(p, xyz, val, pre, 3*i+il[j][0], 3*i+il[j][0], true);
    }
    opt->boundary++;

    for(int i = 0; i < numNodes; i++)
      delete verts[i];
    return;
  }

  for(int i = 0; i < numNodes / 3; i++)
    addScalarTriangle(p, xyz, val, pre, 3*i, 3*i+1, 3*i+2);
}

static void addOutlineTetrahedron(PView *p, double **xyz,
                                  unsigned int color, bool pre)
{
  const int it[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
  for(int i = 0; i < 4; i++)
    addOutlineTriangle(p, xyz, color, pre, it[i][0], it[i][1], it[i][2]);
}

static void addScalarTetrahedron(PView *p, double **xyz,
                                 double **val, bool pre, int i0=0,
                                 int i1=1, int i2=2, int i3=3)
{
  PViewOptions *opt = p->getOptions();

  const int it[4][3] = {{i0, i2, i1}, {i0, i1, i3}, {i0, i3, i2}, {i3, i1, i2}};

  if(opt->boundary > 0 ||
     opt->intervalsType == PViewOptions::Continuous ||
     opt->intervalsType == PViewOptions::Discrete){
    bool skin = (opt->boundary > 0) ? false : opt->drawSkinOnly;
    opt->boundary--;
    for(int i = 0; i < 4; i++)
      addScalarTriangle(p, xyz, val, pre, it[i][0], it[i][1], it[i][2], true, skin);
    opt->boundary++;
    return;
  }

  double vmin = opt->tmpMin, vmax = opt->tmpMax;
  if(opt->saturateValues) saturate(4, val, vmin, vmax, i0, i1, i2, i3);

  double x[4] = {xyz[i0][0], xyz[i1][0], xyz[i2][0], xyz[i3][0]};
  double y[4] = {xyz[i0][1], xyz[i1][1], xyz[i2][1], xyz[i3][1]};
  double z[4] = {xyz[i0][2], xyz[i1][2], xyz[i2][2], xyz[i3][2]};
  double v[4] = {val[i0][0], val[i1][0], val[i2][0], val[i3][0]};

  if(opt->intervalsType == PViewOptions::Iso){
    for(int k = 0; k < opt->nbIso; k++) {
      if(vmin == vmax) k = opt->nbIso / 2;
      double iso = opt->getScaleValue(k, opt->nbIso, vmin, vmax);
      double x2[6], y2[6], z2[6], nn[3];
      int nb = IsoSimplex(x, y, z, v, iso, x2, y2, z2, nn);
      if(nb >= 3){
        unsigned int color = opt->getColor(k, opt->nbIso);
        unsigned int col[3] = {color, color, color};
        for(int j = 2; j < nb; j++){
          double x3[3] = {x2[0], x2[j - 1], x2[j]};
          double y3[3] = {y2[0], y2[j - 1], y2[j]};
          double z3[3] = {z2[0], z2[j - 1], z2[j]};
          SVector3 n[3];
          for(int i = 0; i < 3; i++){
            n[i][0] = nn[0]; n[i][1] = nn[1]; n[i][2] = nn[2];
            if(opt->smoothNormals){
              if(pre) p->normals->add(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
              else p->normals->get(x3[i], y3[i], z3[i], n[i][0], n[i][1], n[i][2]);
            }
          }
          if(!pre) p->va_triangles->add(x3, y3, z3, n, col, 0, false, false);
        }
      }
      if(vmin == vmax) break;
    }
  }
}

static void addOutlineHexahedron(PView *p, double **xyz,
                                 unsigned int color, bool pre)
{
  const int iq[6][4] = {{0, 3, 2, 1}, {0, 1, 5, 4}, {0, 4, 7, 3},
                        {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 6, 7}};

  for(int i = 0; i < 6; i++)
    addOutlineQuadrangle(p, xyz, color, pre, iq[i][0], iq[i][1],
                         iq[i][2], iq[i][3]);
}

static void addScalarHexahedron(PView *p, double **xyz,
                                double **val, bool pre)
{
  PViewOptions *opt = p->getOptions();

  const int iq[6][4] = {{0, 3, 2, 1}, {0, 1, 5, 4}, {0, 4, 7, 3},
                        {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 6, 7}};
  const int is[6][4] = {{0, 1, 3, 4}, {1, 3, 4, 5}, {3, 4, 5, 7},
                        {1, 2, 3, 6}, {3, 1, 6, 5}, {6, 3, 5, 7}};

  if(opt->boundary > 0){
    opt->boundary--;
    for(int i = 0; i < 6; i++)
      addScalarQuadrangle(p, xyz, val, pre, iq[i][0], iq[i][1], iq[i][2], iq[i][3], true);
    opt->boundary++;
    return;
  }

  for(int i = 0; i < 6; i++)
    addScalarTetrahedron(p, xyz, val, pre, is[i][0], is[i][1], is[i][2], is[i][3]);
}

static void addOutlinePrism(PView *p, double **xyz, unsigned int color,
                            bool pre)
{
  const int iq[3][4] = {{0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
  const int it[2][3] = {{0, 2, 1}, {3, 4, 5}};

  for(int i = 0; i < 3; i++)
    addOutlineQuadrangle(p, xyz, color, pre, iq[i][0], iq[i][1], iq[i][2], iq[i][3]);
  for(int i = 0; i < 2; i++)
    addOutlineTriangle(p, xyz, color, pre, it[i][0], it[i][1], it[i][2]);
}

static void addScalarPrism(PView *p, double **xyz, double **val, bool pre)
{
  PViewOptions *opt = p->getOptions();
  const int iq[3][4] = {{0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
  const int it[2][3] = {{0, 2, 1}, {3, 4, 5}};
  const int is[3][4] = {{0, 1, 2, 3}, {3, 4, 5, 2}, {1, 2, 4, 3}};

  if(opt->boundary > 0){
    opt->boundary--;
    for(int i = 0; i < 3; i++)
      addScalarQuadrangle(p, xyz, val, pre, iq[i][0], iq[i][1], iq[i][2], iq[i][3], true);
    for(int i = 0; i < 2; i++)
      addScalarTriangle(p, xyz, val, pre, it[i][0], it[i][1], it[i][2], true);
    opt->boundary++;
    return;
  }

  for(int i = 0; i < 3; i++)
    addScalarTetrahedron(p, xyz, val, pre, is[i][0], is[i][1], is[i][2], is[i][3]);
}

static void addOutlinePyramid(PView *p, double **xyz,
                              unsigned int color, bool pre)
{
  const int it[4][3] = {{0, 1, 4}, {3, 0, 4}, {1, 2, 4}, {2, 3, 4}};

  addOutlineQuadrangle(p, xyz, color, pre, 0, 3, 2, 1);
  for(int i = 0; i < 4; i++)
    addOutlineTriangle(p, xyz, color, pre, it[i][0], it[i][1], it[i][2]);
}

static void addScalarPyramid(PView *p, double **xyz,
                             double **val, bool pre)
{
  PViewOptions *opt = p->getOptions();

  const int it[4][3] = {{0, 1, 4}, {3, 0, 4}, {1, 2, 4}, {2, 3, 4}};
  const int is[2][4] = {{0, 1, 2, 4}, {2, 3, 0, 4}};

  if(opt->boundary > 0){
    opt->boundary--;
    addScalarQuadrangle(p, xyz, val, pre, 0, 3, 2, 1, true);
    for(int i = 0; i < 4; i++)
      addScalarTriangle(p, xyz, val, pre, it[i][0], it[i][1], it[i][2], true);
    opt->boundary++;
    return;
  }

  for(int i = 0; i < 2; i++)
    addScalarTetrahedron(p, xyz, val, pre, is[i][0], is[i][1], is[i][2], is[i][3]);
}

static void addOutlinePolyhedron(PView *p, double **xyz,
                                 unsigned int color, bool pre, int numNodes)
{
  const int it[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
  std::map<MFace, int, Less_Face> triFaces;
  std::vector<MVertex *> verts;
  for(int i = 0; i < numNodes; i++)
    verts.push_back(new MVertex(xyz[i][0], xyz[i][1], xyz[i][2]));
  for(int i = 0; i < numNodes / 4; i++){
    for(int j = 0; j < 4; j++) {
      MFace f(verts[4*i+it[j][0]], verts[4*i+it[j][1]], verts[4*i+it[j][2]]);
      std::map<MFace, int, Less_Face>::iterator ite;
      for(ite = triFaces.begin(); ite != triFaces.end(); ite++)
        if((*ite).first == f) break;
      if(ite == triFaces.end())
        triFaces[f] = 100*i+j;
      else triFaces.erase(ite);
    }
  }
  for(std::map<MFace, int, Less_Face>::iterator ite = triFaces.begin();
      ite != triFaces.end(); ite++){
    int i = (int) (*ite).second / 100;
    int j = (*ite).second % 100;
    addOutlineTriangle(p, xyz, color, pre, 4*i+it[j][0], 4*i+it[j][1], 4*i+it[j][2]);
  }
  for(int i = 0; i < numNodes; i++)
    delete verts[i];
}

static void addScalarPolyhedron(PView *p, double **xyz,
                                double **val, bool pre, int numNodes)
{
  PViewOptions *opt = p->getOptions();

  if(opt->boundary > 0){
    return;
  }

  for(int i = 0; i < numNodes / 4; i++)
    addScalarTetrahedron(p, xyz, val, pre, 4*i, 4*i + 1, 4*i + 2, 4*i + 3);
}

static void addOutlineElement(PView *p, int type, double **xyz,
                              bool pre, int numNodes)
{
  PViewOptions *opt = p->getOptions();
  switch(type){
  case TYPE_PNT: addOutlinePoint(p, xyz, opt->color.point, pre); break;
  case TYPE_LIN: addOutlineLine(p, xyz, opt->color.line, pre); break;
  case TYPE_TRI: addOutlineTriangle(p, xyz, opt->color.triangle, pre); break;
  case TYPE_QUA: addOutlineQuadrangle(p, xyz, opt->color.quadrangle, pre); break;
  case TYPE_POLYG: addOutlinePolygon(p, xyz, opt->color.quadrangle, pre, numNodes); break;
  case TYPE_TET: addOutlineTetrahedron(p, xyz, opt->color.tetrahedron, pre); break;
  case TYPE_HEX: addOutlineHexahedron(p, xyz, opt->color.hexahedron, pre); break;
  case TYPE_PRI: addOutlinePrism(p, xyz, opt->color.prism, pre); break;
  case TYPE_PYR: addOutlinePyramid(p, xyz, opt->color.pyramid, pre); break;
  case TYPE_POLYH: addOutlinePolyhedron(p, xyz, opt->color.pyramid, pre, numNodes); break;
  }
}

static void addScalarElement(PView *p, int type, double **xyz,
                             double **val, bool pre, int numNodes)
{
  switch(type){
  case TYPE_PNT: addScalarPoint(p, xyz, val, pre); break;
  case TYPE_LIN: addScalarLine(p, xyz, val, pre); break;
  case TYPE_TRI: addScalarTriangle(p, xyz, val, pre); break;
  case TYPE_QUA: addScalarQuadrangle(p, xyz, val, pre); break;
  case TYPE_POLYG: addScalarPolygon(p, xyz, val, pre, numNodes); break;
  case TYPE_TET: addScalarTetrahedron(p, xyz, val, pre); break;
  case TYPE_HEX: addScalarHexahedron(p, xyz, val, pre); break;
  case TYPE_PRI: addScalarPrism(p, xyz, val, pre); break;
  case TYPE_PYR: addScalarPyramid(p, xyz, val, pre); break;
  case TYPE_POLYH: addScalarPolyhedron(p, xyz, val, pre, numNodes); break;
  }
}

static void addVectorElement(PView *p, int ient, int iele, int numNodes,
                             int type, double **xyz,
                             double **val, bool pre)
{
  // use adaptive data if available
  PViewData *data = p->getData(true);
  PViewOptions *opt = p->getOptions();

  int numComp2;
  double **val2 = new double*[numNodes];
  for(int i = 0; i < numNodes; i++)
    val2[i] = new double[9];
  getExternalValues(p, opt->externalViewIndex, ient, iele, numNodes,
                    3, val, numComp2, val2);

  if(opt->vectorType == PViewOptions::Displacement){
    for(int i = 0; i < numNodes; i++)
      val2[i][0] = ComputeScalarRep(numComp2, val2[i]);

    // add scalar element with correct min/max
    double min = opt->tmpMin, max = opt->tmpMax;
    opt->tmpMin = opt->externalMin;
    opt->tmpMax = opt->externalMax;
    addScalarElement(p, type, xyz, val2, pre, numNodes);
    opt->tmpMin = min;
    opt->tmpMax = max;

    // add point trajectories
    if(!pre && numNodes == 1 && opt->timeStep > 0 && opt->lineWidth){
      for(int ts = 0; ts < opt->timeStep; ts++){
        int numComp = data->getNumComponents(ts, ient, iele);
        double xyz0[3], dxyz[3][2] = {{0., 0.}, {0., 0.}, {0., 0.}};
        data->getNode(ts, ient, iele, 0, xyz0[0], xyz0[1], xyz0[2]);
        for(int j = 0; j < 3; j++){
          int comp = opt->forceNumComponents ? opt->componentMap[j] : j;
          if(comp >= 0 && comp < numComp){
            data->getValue(ts, ient, iele, 0, comp, dxyz[j][0]);
            data->getValue(ts + 1, ient, iele, 0, comp, dxyz[j][1]);
          }
        }
        unsigned int col[2];
        double norm[2];
        for(int i = 0; i < 2; i++){
          norm[i] = sqrt(dxyz[0][i] * dxyz[0][i] +
                         dxyz[1][i] * dxyz[1][i] +
                         dxyz[2][i] * dxyz[2][i]);
          col[i] = opt->getColor(norm[i], opt->tmpMin, opt->tmpMax);
        }
        for(int j = 0; j < 3; j++){
          dxyz[j][0] = xyz0[j] + dxyz[j][0] * opt->displacementFactor;
          dxyz[j][1] = xyz0[j] + dxyz[j][1] * opt->displacementFactor;
        }
        SVector3 n[2];
        getLineNormal(p, dxyz[0], dxyz[1], dxyz[2], norm, n, true);
        p->va_lines->add(dxyz[0], dxyz[1], dxyz[2], n, col, 0, false);
      }
    }
    for(int i = 0; i < numNodes; i++)
      delete [] val2[i];
    delete [] val2;
    return;
  }

  if(pre){
    for(int i = 0; i < numNodes; i++)
      delete [] val2[i];
    delete [] val2;
    return;
  }

  if(opt->glyphLocation == PViewOptions::Vertex){
    for(int i = 0; i < numNodes; i++){
      double v2 = opt->saturateValues ?
        saturateVector(val[i], numComp2, val2[i], opt->externalMin, opt->externalMax) :
        ComputeScalarRep(numComp2, val2[i]);
      if(v2 >= opt->externalMin && v2 <= opt->externalMax){
        unsigned int color = opt->getColor(v2, opt->externalMin, opt->externalMax, false,
                                           (opt->intervalsType == PViewOptions::Discrete) ?
                                           opt->nbIso : -1);
        unsigned int col[2] = {color, color};
        double dxyz[3][2];
        for(int j = 0; j < 3; j++){
          dxyz[j][0] = xyz[i][j];
          dxyz[j][1] = val[i][j];
        }
        p->va_vectors->add(dxyz[0], dxyz[1], dxyz[2], 0, col, 0, false);
      }
    }
  }

  if(opt->glyphLocation == PViewOptions::COG){
    SPoint3 pc(0., 0., 0.);
    double d[3] = {0., 0., 0.};
    double d2[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};

    for(int i = 0; i < numNodes; i++){
      pc += SPoint3(xyz[i][0], xyz[i][1], xyz[i][2]);
      for(int j = 0; j < 3; j++) d[j] += val[i][j];
      for(int j = 0; j < numComp2; j++) d2[j] += val2[i][j];
    }
    pc /= (double)numNodes;
    for(int j = 0; j < 3; j++) d[j] /= (double)numNodes;
    for(int j = 0; j < numComp2; j++) d2[j] /= (double)numNodes;

    // need tolerance since we compare computed results (the average)
    // instead of the raw data used to compute bounds
    double v2 = opt->saturateValues ?
      saturateVector(d, numComp2, d2, opt->externalMin, opt->externalMax) :
      ComputeScalarRep(numComp2, d2);
    if(v2 >= opt->externalMin * (1. - 1.e-15) &&
       v2 <= opt->externalMax * (1. + 1.e-15)){
      unsigned int color = opt->getColor(v2, opt->externalMin, opt->externalMax, false,
                                         (opt->intervalsType == PViewOptions::Discrete) ?
                                         opt->nbIso : -1);
      unsigned int col[2] = {color, color};
      double dxyz[3][2];
      for(int i = 0; i < 3; i++){
        dxyz[i][0] = pc[i];
        dxyz[i][1] = d[i];
      }
      p->va_vectors->add(dxyz[0], dxyz[1], dxyz[2], 0, col, 0, false);
    }
  }
  for(int i = 0; i < numNodes; i++)
    delete [] val2[i];
  delete [] val2;
}

static void addTensorElement(PView *p, int iEnt, int iEle, int numNodes, int type,
                             double **xyz, double **val,
                             bool pre)
{
  PViewOptions *opt = p->getOptions();
  fullMatrix <double> tensor(3,3);
  fullVector<double> S(3), imS (3);
  fullMatrix<double> V(3,3);
  fullMatrix <double> rightV(3,3);

  if(opt->tensorType == PViewOptions::VonMises){
    for(int i = 0; i < numNodes; i++)
      val[i][0] = ComputeVonMises(val[i]);
    addScalarElement(p, type, xyz, val, pre, numNodes);
  }
  else if (opt->tensorType == PViewOptions::Ellipse ||
           opt->tensorType == PViewOptions::Ellipsoid) {
    if(opt->glyphLocation == PViewOptions::Vertex){
      double vval[3][4]= {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}};
      for(int i = 0; i < numNodes; i++){
        for (int j = 0; j < 3; j++) {
          tensor(j,0) = val [i][0+j*3];
          tensor(j,1) = val [i][1+j*3];
          tensor(j,2) = val [i][2+j*3];
        }
        tensor.eig(S, imS, V, rightV, false);
        for (int k = 0; k < 3; k++) {
          vval[k][0] = xyz[i][k];
          for (int j = 0; j < 3; j++) {
            vval[k][j+1] = V(k,j)*S(j);
          }
        }
        double lmax = std::max(S(0), std::max(S(1), S(2)));
        unsigned int color = opt->getColor
          (lmax, opt->tmpMin, opt->tmpMax, false,
           (opt->intervalsType == PViewOptions::Discrete) ?  opt->nbIso : -1);
        unsigned int col[4] = {color, color, color, color};
        p->va_ellipses->add(vval[0], vval[1], vval[2], 0, col, 0, false);
      }
    }
    else if(opt->glyphLocation == PViewOptions::COG){
      double vval[3][4]= {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}};
      for(int i = 0; i < numNodes; i++) {
        for (int j = 0; j < 3; j++) {
          tensor(j,0) = val [i][0+j*3];
          tensor(j,1) = val [i][1+j*3];
          tensor(j,2) = val [i][2+j*3];
        }
        tensor.eig(S, imS, V, rightV, false);
        for (int j = 0; j < 3; j++) {
          vval[0][j+1] += V(0,j)*S(j)/numNodes;
          vval[1][j+1] += V(1,j)*S(j)/numNodes;
          vval[2][j+1] += V(2,j)*S(j)/numNodes;
        }
        vval[0][0] += xyz[i][0]/numNodes;
        vval[1][0] += xyz[i][1]/numNodes;
        vval[2][0] += xyz[i][2]/numNodes;
      }
      double lmax = std::max(S(0), std::max(S(1), S(2)));
      unsigned int color = opt->getColor
        (lmax, opt->tmpMin, opt->tmpMax, false,
         (opt->intervalsType == PViewOptions::Discrete) ?  opt->nbIso : -1);
      unsigned int col[4] = {color, color, color, color};
      p->va_ellipses->add(vval[0], vval[1], vval[2], 0, col, 0, false);
    }
  }
  else {
    double **vval[3] = {new double*[numNodes],
                        new double*[numNodes],
                        new double*[numNodes]};
    for(int i = 0; i < 3; i++)
      for(int j = 0; j < numNodes; j++)
        vval[i][j] = new double[3];
    for(int i = 0; i < numNodes; i++) {
      for (int j = 0; j < 3; j++) {
        tensor(j,0) = val [i][0+j*3];
        tensor(j,1) = val [i][1+j*3];
        tensor(j,2) = val [i][2+j*3];
      }
      tensor.eig(S, imS, V, rightV, opt->tensorType != PViewOptions::EigenVectors);
      if (PViewOptions::MinEigenValue == opt->tensorType)
        val[i][0] = S(0);
      else if (PViewOptions::MaxEigenValue == opt->tensorType)
        val[i][0] = S(2);
      else if (PViewOptions::EigenVectors == opt->tensorType) {
        for (int j = 0; j < 3; j++) {
          vval[0][i][j] = V(j,0)*S(0);
          vval[1][i][j] = V(j,1)*S(1);
          vval[2][i][j] = V(j,2)*S(2);
        }
      }
    }

    if (PViewOptions::EigenVectors == opt->tensorType) {
      addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[0], pre);
      addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[1], pre);
      addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[2], pre);
    }
    else
      addScalarElement(p, type, xyz, val, pre, numNodes);
    for(int i = 0; i < 3; i++){
      for(int j = 0; j < numNodes; j++)
        delete [] vval[i][j];
      delete [] vval[i];
    }
  }
}

static void addElementsInArrays(PView *p, bool preprocessNormalsOnly)
{
  static int numNodesError = 0, numCompError = 0;

  // use adaptive data if available
  PViewData *data = p->getData(true);
  PViewOptions *opt = p->getOptions();

  opt->tmpBBox.reset();

  int NMAX = PVIEW_NMAX;
  double **xyz = new double*[NMAX];
  double **val = new double*[NMAX];
  for(int i = 0; i < NMAX; i++){
    xyz[i] = new double[3];
    val[i] = new double[9];
  }
  for(int ent = 0; ent < data->getNumEntities(opt->timeStep); ent++){
    if(data->skipEntity(opt->timeStep, ent)) continue;
    for(int i = 0; i < data->getNumElements(opt->timeStep, ent); i++){
      if(data->skipElement(opt->timeStep, ent, i, true, opt->sampling)) continue;
      int type = data->getType(opt->timeStep, ent, i);
      if(opt->skipElement(type)) continue;
      int numComp = data->getNumComponents(opt->timeStep, ent, i);
      int numNodes = data->getNumNodes(opt->timeStep, ent, i);
      if(numNodes > PVIEW_NMAX){
        if(type == TYPE_POLYG || type == TYPE_POLYH){
          if(numNodes > NMAX){
            for(int j = 0; j < NMAX; j++){
              delete [] xyz[j];
              delete [] val[j];
            }
            delete [] xyz;
            delete [] val;
            NMAX = numNodes;
            xyz = new double*[NMAX];
            val = new double*[NMAX];
            for(int j = 0; j < NMAX; j++){
              xyz[j] = new double[3];
              val[j] = new double[9];
            }
          }
        }
        else {
          if(numNodesError != numNodes){
            numNodesError = numNodes;
            Msg::Error("You should never draw views with > %d nodes per element: use"
                       "'Adapt visualization grid' to view high-order datasets!",
                       PVIEW_NMAX);
          }
          continue;
        }
      }
      if((numComp > 9 && !opt->forceNumComponents) || opt->forceNumComponents > 9){
        if(numCompError != numComp) {
          numCompError = numComp;
          Msg::Error("You should never draw views with > 9 values per node: use"
                     "'Adapt visualization grid' to view high-order datasets!");
        }
        continue;
      }
      for(int j = 0; j < numNodes; j++){
        data->getNode(opt->timeStep, ent, i, j, xyz[j][0], xyz[j][1], xyz[j][2]);
        if(opt->forceNumComponents){
          for(int k = 0; k < opt->forceNumComponents; k++){
            int comp = opt->componentMap[k];
            if(comp >= 0 && comp < numComp)
              data->getValue(opt->timeStep, ent, i, j, comp, val[j][k]);
            else
              val[j][k] = 0.;
          }
        }
        else
          for(int k = 0; k < numComp; k++)
            data->getValue(opt->timeStep, ent, i, j, k, val[j][k]);
      }
      if(opt->forceNumComponents) numComp = opt->forceNumComponents;

      changeCoordinates(p, ent, i, numNodes, type, numComp, xyz, val);
      int dim = data->getDimension(opt->timeStep, ent, i);
      if(!isElementVisible(opt, dim, numNodes, xyz)) continue;

      for(int j = 0; j < numNodes; j++)
        opt->tmpBBox += SPoint3(xyz[j][0], xyz[j][1], xyz[j][2]);

      if(opt->showElement && !data->useGaussPoints())
        addOutlineElement(p, type, xyz, preprocessNormalsOnly, numNodes);

      if(opt->intervalsType != PViewOptions::Numeric){
        if(data->useGaussPoints()){
          for(int j = 0; j < numNodes; j++){
            double *x2 = new double[3]; double **xyz2 = &x2;
            double *v2 = new double[9]; double **val2 = &v2;
            xyz2[0][0] = xyz[j][0];
            xyz2[0][1] = xyz[j][1];
            xyz2[0][2] = xyz[j][2];
            for(int k = 0; k < numComp; k++)
              val2[0][k] = val[j][k];
            if(numComp == 1 && opt->drawScalars)
              addScalarElement(p, TYPE_PNT, xyz2, val2, preprocessNormalsOnly, numNodes);
            else if(numComp == 3 && opt->drawVectors)
              addVectorElement(p, ent, i, 1, TYPE_PNT, xyz2, val2, preprocessNormalsOnly);
            else if(numComp == 9 && opt->drawTensors)
              addTensorElement(p, ent, i, 1, TYPE_PNT, xyz2, val2, preprocessNormalsOnly);
            delete [] xyz2[0];
            delete [] val2[0];
          }
        }
        else if(numComp == 1 && opt->drawScalars)
          addScalarElement(p, type, xyz, val, preprocessNormalsOnly, numNodes);
        else if(numComp == 3 && opt->drawVectors)
          addVectorElement(p, ent, i, numNodes, type, xyz, val, preprocessNormalsOnly);
        else if(numComp == 9 && opt->drawTensors)
          addTensorElement(p, ent, i, numNodes, type, xyz, val, preprocessNormalsOnly);
      }
    }
  }
  for(int j = 0; j < NMAX; j++){
    delete [] xyz[j];
    delete [] val[j];
  }
  delete [] xyz;
  delete [] val;
}

class initPView {
 private:
  // we try to estimate how many primitives will end up in the vertex
  // arrays, since reallocating the arrays takes a huge amount of time
  // on Windows/Cygwin
  int _estimateIfClipped(PView *p, int num)
  {
    if(CTX::instance()->clipWholeElements &&
       CTX::instance()->clipOnlyDrawIntersectingVolume){
      PViewOptions *opt = p->getOptions();
      for(int clip = 0; clip < 6; clip++){
        if(opt->clip & (1 << clip))
          return (int)sqrt((double)num);
      }
    }
    return num;
  }
  int _estimateNumPoints(PView *p)
  {
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();
    int heuristic = data->getNumPoints(opt->timeStep);
    return heuristic + 10000;
  }
  int _estimateNumLines(PView *p)
  {
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();
    int heuristic = data->getNumLines(opt->timeStep);
    return heuristic + 10000;
  }
  int _estimateNumTriangles(PView *p)
  {
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();
    int tris = data->getNumTriangles(opt->timeStep);
    int quads = data->getNumQuadrangles(opt->timeStep);
    int polygs = data->getNumPolygons(opt->timeStep);
    int tets = data->getNumTetrahedra(opt->timeStep);
    int prisms = data->getNumPrisms(opt->timeStep);
    int pyrs = data->getNumPyramids(opt->timeStep);
    int hexas = data->getNumHexahedra(opt->timeStep);
    int polyhs = data->getNumPolyhedra(opt->timeStep);
    int heuristic = 0;
    if(opt->intervalsType == PViewOptions::Iso)
      heuristic = (tets + prisms + pyrs + hexas + polyhs) / 10;
    else if(opt->intervalsType == PViewOptions::Continuous)
      heuristic = (tris + 2 * quads + 3 * polygs + 6 * tets +
                   8 * prisms + 6 * pyrs + 12 * hexas + 10 * polyhs);
    else if(opt->intervalsType == PViewOptions::Discrete)
      heuristic = (tris + 2 * quads + 3 * polygs + 6 * tets +
                   8 * prisms + 6 * pyrs + 12 * hexas + 10 * polyhs) * 2;
    heuristic = _estimateIfClipped(p, heuristic);
    return heuristic + 10000;
  }
  int _estimateNumVectors(PView *p)
  {
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();
    int heuristic = data->getNumVectors(opt->timeStep);
    heuristic = _estimateIfClipped(p, heuristic);
    return heuristic + 1000;
  }
  int _estimateNumEllipses(PView *p)
  {
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();
    int heuristic = data->getNumTensors(opt->timeStep);
    heuristic = _estimateIfClipped(p, heuristic);
    return heuristic + 1000;
  }
 public:
  void operator () (PView *p)
  {
    // use adaptive data if available
    PViewData *data = p->getData(true);
    PViewOptions *opt = p->getOptions();

    if(data->getDirty() || !data->getNumTimeSteps() || !p->getChanged()) return;
    if(!opt->visible || opt->type != PViewOptions::Plot3D) return;

    p->deleteVertexArrays();

    if(data->isRemote()){
      // FIXME: need to rewrite option code and add nice serialization
      std::string fileName = CTX::instance()->homeDir + CTX::instance()->tmpFileName;
      PrintOptions(0, GMSH_FULLRC, 0, 0, fileName.c_str());
      std::string options = ConvertFileToString(fileName);
      data->fillRemoteVertexArrays(options);
      return;
    }

    if(opt->useGenRaise) opt->createGeneralRaise();

    if(opt->rangeType == PViewOptions::Custom){
      opt->tmpMin = opt->customMin;
      opt->tmpMax = opt->customMax;
    }
    else if(opt->rangeType == PViewOptions::PerTimeStep){
      opt->tmpMin = data->getMin(opt->timeStep);
      opt->tmpMax = data->getMax(opt->timeStep);
    }
    else{
      // FIXME: this is not perfect for multi-step adaptive views, as
      // we don't have the correct min/max info for the other steps
      opt->tmpMin = data->getMin();
      opt->tmpMax = data->getMax();
    }

    p->va_points = new VertexArray(1, _estimateNumPoints(p));
    p->va_lines = new VertexArray(2, _estimateNumLines(p));
    p->va_triangles = new VertexArray(3, _estimateNumTriangles(p));
    p->va_vectors = new VertexArray(2, _estimateNumVectors(p));
    p->va_ellipses = new VertexArray(4, _estimateNumEllipses(p));

    if(p->normals) delete p->normals;
    p->normals = new smooth_normals(opt->angleSmoothNormals);

    if(opt->smoothNormals) addElementsInArrays(p, true);
    addElementsInArrays(p, false);

    p->va_points->finalize();
    p->va_lines->finalize();
    p->va_triangles->finalize();
    p->va_vectors->finalize();
    p->va_ellipses->finalize();

    Msg::Debug("%d vertices in vertex arrays (%g Mb)", p->va_points->getNumVertices() +
               p->va_lines->getNumVertices() + p->va_triangles->getNumVertices…