OctreePost.cpp

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

#include "Octree.h"
#include "OctreePost.h"
#include "PView.h"
#include "PViewData.h"
#include "PViewDataList.h"
#include "PViewDataGModel.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "shapeFunctions.h"
#include "GModel.h"
#include "MElement.h"
#include "Context.h"

// helper routines for list-based views

static void minmax(int n, double *X, double *Y, double *Z,
                   double *min, double *max)
{
  min[0] = X[0];
  min[1] = Y[0];
  min[2] = Z[0];
  max[0] = X[0];
  max[1] = Y[0];
  max[2] = Z[0];
  for(int i = 1; i < n; i++) {
    min[0] = (X[i] < min[0]) ? X[i] : min[0];
    min[1] = (Y[i] < min[1]) ? Y[i] : min[1];
    min[2] = (Z[i] < min[2]) ? Z[i] : min[2];
    max[0] = (X[i] > max[0]) ? X[i] : max[0];
    max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
    max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
  }

  // make bounding boxes larger up to (absolute) geometrical tolerance
  double eps = CTX::instance()->geom.tolerance;
  for(int i = 0; i < 3; i++){
    min[i] -= eps;
    max[i] += eps;
  }
}

static void centroid(int n, double *X, double *Y, double *Z, double *c)
{
  const double oc = 1./(double)n;
  c[0] = X[0];
  c[1] = Y[0];
  c[2] = Z[0];
  for(int i = 1; i < n; i++) {
    c[0] += X[i];
    c[1] += Y[i];
    c[2] += Z[i];
  }
  c[0] *= oc;
  c[1] *= oc;
  c[2] *= oc;
}

void linBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[2], *Z = &X[4];
  minmax(2, X, Y, Z, min, max);
}

void triBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[3], *Z = &X[6];
  minmax(3, X, Y, Z, min, max);
}

void quaBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8];
  minmax(4, X, Y, Z, min, max);
}

void tetBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8];
  minmax(4, X, Y, Z, min, max);
}

void hexBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[8], *Z = &X[16];
  minmax(8, X, Y, Z, min, max);
}

void priBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[6], *Z = &X[12];
  minmax(6, X, Y, Z, min, max);
}

void pyrBB(void *a, double *min, double *max)
{
  double *X = (double*) a, *Y = &X[5], *Z = &X[10];
  minmax(5, X, Y, Z, min, max);
}

int linInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[2], *Z = &X[4], uvw[3];
  line lin(X, Y, Z);
  lin.xyz2uvw(x, uvw);
  return lin.isInside(uvw[0], uvw[1], uvw[2]);
}

int triInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[3], *Z = &X[6], uvw[3];
  triangle tri(X, Y, Z);
  tri.xyz2uvw(x, uvw);
  return tri.isInside(uvw[0], uvw[1], uvw[2]);
}

int quaInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
  quadrangle qua(X, Y, Z);
  qua.xyz2uvw(x, uvw);
  return qua.isInside(uvw[0], uvw[1], uvw[2]);
}

int tetInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
  tetrahedron tet(X, Y, Z);
  tet.xyz2uvw(x, uvw);
  return tet.isInside(uvw[0], uvw[1], uvw[2]);
}

int hexInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[8], *Z = &X[16], uvw[3];
  hexahedron hex(X, Y, Z);
  hex.xyz2uvw(x, uvw);
  return hex.isInside(uvw[0], uvw[1], uvw[2]);
}

int priInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[6], *Z = &X[12], uvw[3];
  prism pri(X, Y, Z);
  pri.xyz2uvw(x, uvw);
  return pri.isInside(uvw[0], uvw[1], uvw[2]);
}

int pyrInEle(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[5], *Z = &X[10], uvw[3];
  pyramid pyr(X, Y, Z);
  pyr.xyz2uvw(x, uvw);
  return pyr.isInside(uvw[0], uvw[1], uvw[2]);
}

void linCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[2], *Z = &X[4];
  centroid(2, X, Y, Z, x);
}

void triCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[3], *Z = &X[6];
  centroid(3, X, Y, Z, x);
}

void quaCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8];
  centroid(4, X, Y, Z, x);
}

void tetCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[4], *Z = &X[8];
  centroid(4, X, Y, Z, x);
}

void hexCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[8], *Z = &X[16];
  centroid(8, X, Y, Z, x);
}

void priCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[6], *Z = &X[12];
  centroid(6, X, Y, Z, x);
}

void pyrCentroid(void *a, double *x)
{
  double *X = (double*) a, *Y = &X[5], *Z = &X[10];
  centroid(5, X, Y, Z, x);
}

static void addListOfStuff(Octree *o, std::vector<double> &l, int nbelm)
{
  for(unsigned int i = 0; i < l.size(); i += nbelm)
    Octree_Insert(&l[i], o);
}

// OctreePost implementation

OctreePost::~OctreePost() 
{
  Octree_Delete(_SL); Octree_Delete(_VL); Octree_Delete(_TL);
  Octree_Delete(_ST); Octree_Delete(_VT); Octree_Delete(_TT);
  Octree_Delete(_SQ); Octree_Delete(_VQ); Octree_Delete(_TQ);
  Octree_Delete(_SS); Octree_Delete(_VS); Octree_Delete(_TS);
  Octree_Delete(_SH); Octree_Delete(_VH); Octree_Delete(_TH);
  Octree_Delete(_SI); Octree_Delete(_VI); Octree_Delete(_TI);
  Octree_Delete(_SY); Octree_Delete(_VY); Octree_Delete(_TY);
}

OctreePost::OctreePost(PView *v) 
{
  _create(v->getData(true)); // use adaptive data if available
}

OctreePost::OctreePost(PViewData *data) 
{
  _create(data);
}

void OctreePost::_create(PViewData *data)
{
  _SL = _VL = _TL = _ST = _VT = _TT = _SQ = _VQ = _TQ = 0; 
  _SS = _VS = _TS = _SH = _VH = _TH = _SI = _VI = _TI = 0;
  _SY = _VY = _TY = 0;
  _theViewDataList = 0;
  _theViewDataGModel = 0;

  _theViewDataGModel = dynamic_cast<PViewDataGModel*>(data);

  if(_theViewDataGModel) return; // the octree is already available in the model

  _theViewDataList = dynamic_cast<PViewDataList*>(data);

  if(_theViewDataList){
    PViewDataList *l = _theViewDataList;

    if(l->haveInterpolationMatrices() && !l->isAdapted()){
      Msg::Error("Cannot create octree for non-adapted high-order list-based view: you need");
      Msg::Error("to select 'Adapt visualization grid' first");
      return;
    }

    SBoundingBox3d bb = l->getBoundingBox();
    // make bounding box larger up to (absolute) geometrical tolerance
    double eps = CTX::instance()->geom.tolerance;
    SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
    bbmin -= bbeps;
    bbmax += bbeps;
    double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
    double size[3] = {bbmax.x() - bbmin.x(),
                      bbmax.y() - bbmin.y(),
                      bbmax.z() - bbmin.z()};
    const int maxElePerBucket = 100; // memory vs. speed trade-off
    
    _SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
    addListOfStuff(_SL, l->SL, 6 + 2 * l->getNumTimeSteps());
    Octree_Arrange(_SL);
    _VL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
    addListOfStuff(_VL, l->VL, 6 + 6 * l->getNumTimeSteps());
    Octree_Arrange(_VL);
    _TL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
    addListOfStuff(_TL, l->TL, 6 + 18 * l->getNumTimeSteps());
    Octree_Arrange(_TL);
    
    _ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
    addListOfStuff(_ST, l->ST, 9 + 3 * l->getNumTimeSteps());
    Octree_Arrange(_ST);
    _VT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
    addListOfStuff(_VT, l->VT, 9 + 9 * l->getNumTimeSteps());
    Octree_Arrange(_VT);
    _TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
    addListOfStuff(_TT, l->TT, 9 + 27 * l->getNumTimeSteps());
    Octree_Arrange(_TT);
    
    _SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
    addListOfStuff(_SQ, l->SQ, 12 + 4 * l->getNumTimeSteps());
    Octree_Arrange(_SQ);
    _VQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
    addListOfStuff(_VQ, l->VQ, 12 + 12 * l->getNumTimeSteps());
    Octree_Arrange(_VQ);
    _TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
    addListOfStuff(_TQ, l->TQ, 12 + 36 * l->getNumTimeSteps());
    Octree_Arrange(_TQ);
    
    _SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
    addListOfStuff(_SS, l->SS, 12 + 4 * l->getNumTimeSteps());
    Octree_Arrange(_SS);
    _VS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
    addListOfStuff(_VS, l->VS, 12 + 12 * l->getNumTimeSteps());
    Octree_Arrange(_VS);
    _TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
    addListOfStuff(_TS, l->TS, 12 + 36 * l->getNumTimeSteps());
    Octree_Arrange(_TS);
    
    _SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
    addListOfStuff(_SH, l->SH, 24 + 8 * l->getNumTimeSteps());
    Octree_Arrange(_SH);
    _VH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
    addListOfStuff(_VH, l->VH, 24 + 24 * l->getNumTimeSteps());
    Octree_Arrange(_VH);
    _TH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
    addListOfStuff(_TH, l->TH, 24 + 72 * l->getNumTimeSteps());
    Octree_Arrange(_TH);
    
    _SI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
    addListOfStuff(_SI, l->SI, 18 + 6 * l->getNumTimeSteps());
    Octree_Arrange(_SI);
    _VI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
    addListOfStuff(_VI, l->VI, 18 + 18 * l->getNumTimeSteps());
    Octree_Arrange(_VI);
    _TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
    addListOfStuff(_TI, l->TI, 18 + 54 * l->getNumTimeSteps());
    Octree_Arrange(_TI);
    
    _SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
    addListOfStuff(_SY, l->SY, 15 + 5 * l->getNumTimeSteps());
    Octree_Arrange(_SY);
    _VY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
    addListOfStuff(_VY, l->VY, 15 + 15 * l->getNumTimeSteps());
    Octree_Arrange(_VY);
    _TY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
    addListOfStuff(_TY, l->TY, 15 + 45 * l->getNumTimeSteps());
    Octree_Arrange(_TY);
  }
}

bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp, 
                           double P[3], int step, double *values,
                           double *elementSize)
{
  if(!in) return false;

  double *X = (double*)in, *Y = &X[nbNod], *Z = &X[2*nbNod], *V = &X[3*nbNod], U[3];

  elementFactory factory;
  element *e = factory.create(nbNod, dim, X, Y, Z);
  if(!e) return false;

  e->xyz2uvw(P, U);
  if(step < 0){
    for(int i = 0; i < _theViewDataList->getNumTimeSteps(); i++)
      for(int j = 0; j < nbComp; j++)
        values[nbComp * i + j] = e->interpolate(&V[nbNod * nbComp * i + j], 
                                                U[0], U[1], U[2], nbComp);
  }
  else{
    for(int j = 0; j < nbComp; j++)
      values[j] = e->interpolate(&V[nbNod * nbComp * step + j], 
                                 U[0], U[1], U[2], nbComp);
  }

  if(elementSize) *elementSize = e->maxEdgeLength();

  delete e;
  return true;
} 

bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep,
                           double *values, double *elementSize)
{
  if(!in) return false;

  if(_theViewDataGModel->getNumComponents(0, 0, 0) != nbComp) return false;

  MElement *e = (MElement*)in;

  std::vector<int> dataIndex(e->getNumVertices());
  if(_theViewDataGModel->getType() == PViewDataGModel::NodeData)
    for(int i = 0; i < e->getNumVertices(); i++)
      dataIndex[i] = e->getVertex(i)->getNum();
  else
    for(int i = 0; i < e->getNumVertices(); i++)
      dataIndex[i] = e->getNum();
    
  double U[3];
  e->xyz2uvw(P, U);

  std::vector<double> nodeval(e->getNumVertices() * 9);
  for(int step = 0; step < _theViewDataGModel->getNumTimeSteps(); step++){
    if(!_theViewDataGModel->hasTimeStep(step)) continue;
    if(timestep < 0 || step == timestep){
      for(int nod = 0; nod < e->getNumVertices(); nod++){
        for(int comp = 0; comp < nbComp; comp++)
          _theViewDataGModel->getValueByIndex(step, dataIndex[nod], nod, comp, 
                                              nodeval[nod * nbComp + comp]);
      }
      for(int comp = 0; comp < nbComp; comp++){
        double val = e->interpolate(&nodeval[comp], U[0], U[1], U[2], nbComp);
        if(timestep < 0)
          values[nbComp * step + comp] = val;
        else
          values[comp] = val;
      }
    }
  }
  
  if(elementSize) *elementSize = e->maxEdge();
  return true;
} 

bool OctreePost::searchScalar(double x, double y, double z, double *values, 
                              int step, double *size)
{
  double P[3] = {x, y, z};

  if(step < 0){
    int numSteps = 1;
    if(_theViewDataList) numSteps = _theViewDataList->getNumTimeSteps();
    else if(_theViewDataGModel) numSteps = _theViewDataGModel->getNumTimeSteps();
    for(int i = 0; i < numSteps; i++){
      values[i] = 0.0; 
    }
  }
  else
    values[0] = 0.0;

  if(_theViewDataList){
    if(_getValue(Octree_Search(P, _SS), 3, 4, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _SH), 3, 8, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _SI), 3, 6, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _SY), 3, 5, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _ST), 2, 3, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _SQ), 2, 4, 1, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _SL), 1, 2, 1, P, step, values, size)) return true;
  }
  else if(_theViewDataGModel){
    GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
    if(m){
      SPoint3 pt(P);
      if(_getValue(m->getMeshElementByCoord(pt), 1, P, step, values, size)) return true;
    }
  }
  
  return false;
}

bool OctreePost::searchScalarWithTol(double x, double y, double z, double *values, 
                                     int step, double *size, double tol)
{
  bool a = searchScalar(x, y, z, values, step, size);
  if(!a && tol != 0.){
    double oldtol1 = element::getTolerance();
    double oldtol2 = MElement::getTolerance();
    element::setTolerance(tol);
    MElement::setTolerance(tol);
    a = searchScalar(x, y, z, values, step, size);
    element::setTolerance(oldtol1);
    MElement::setTolerance(oldtol2);
  }    
  return a;
}

bool OctreePost::searchVector(double x, double y, double z, double *values, 
                              int step, double *size)
{
  double P[3] = {x, y, z};

  if(step < 0){
    int numSteps = 1;
    if(_theViewDataList) numSteps = _theViewDataList->getNumTimeSteps();
    else if(_theViewDataGModel) numSteps = _theViewDataGModel->getNumTimeSteps();
    for(int i = 0; i < 3 * numSteps; i++)
      values[i] = 0.0; 
  }
  else
    for(int i = 0; i < 3; i++)
      values[i] = 0.0;

  if(_theViewDataList){
    if(_getValue(Octree_Search(P, _VS), 3, 4, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VH), 3, 8, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VI), 3, 6, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VY), 3, 5, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VT), 2, 3, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VQ), 2, 4, 3, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _VL), 1, 2, 3, P, step, values, size)) return true;
  }
  else if(_theViewDataGModel){
    GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
    if(m){
      SPoint3 pt(P);
      if(_getValue(m->getMeshElementByCoord(pt), 3, P, step, values, size)) return true;
    }
  }
  
  return false;
}

bool OctreePost::searchTensor(double x, double y, double z, double *values, 
                              int step, double *size)
{
  double P[3] = {x, y, z};

  if(step < 0){
    int numSteps = 1;
    if(_theViewDataList) numSteps = _theViewDataList->getNumTimeSteps();
    else if(_theViewDataGModel) numSteps = _theViewDataGModel->getNumTimeSteps();
    for(int i = 0; i < 9 * numSteps; i++)
      values[i] = 0.0; 
  }
  else
    for(int i = 0; i < 9; i++)
      values[i] = 0.0;

  if(_theViewDataList){
    if(_getValue(Octree_Search(P, _TS), 3, 4, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TH), 3, 8, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TI), 3, 6, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TY), 3, 5, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TT), 2, 3, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TQ), 2, 4, 9, P, step, values, size)) return true;
    if(_getValue(Octree_Search(P, _TL), 1, 2, 9, P, step, values, size)) return true;
  }
  else if(_theViewDataGModel){
    GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
    if(m){
      SPoint3 pt(P);
      if(_getValue(m->getMeshElementByCoord(pt), 9, P, step, values, size)) return true;
    }
  }
  
  return false;
}