Skip to content
Snippets Groups Projects
Select Git revision
  • 3f73249a8770896b03a093b3fe10bb7cc62f42e0
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

PViewDataGModel.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    PViewDataGModel.cpp 15.26 KiB
    // Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to <gmsh@geuz.org>.
    
    #include "PViewDataGModel.h"
    #include "MLine.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MTetrahedron.h"
    #include "MHexahedron.h"
    #include "MPrism.h"
    #include "MPyramid.h"
    #include "Numeric.h"
    #include "GmshMessage.h"
    
    PViewDataGModel::PViewDataGModel(DataType type) 
      : PViewData(), _min(VAL_INF), _max(-VAL_INF), _type(type)
    {
    }
    
    PViewDataGModel::~PViewDataGModel()
    {
      for(unsigned int i = 0; i < _steps.size(); i++) delete _steps[i];
    }
    
    bool PViewDataGModel::finalize()
    {
      _min = VAL_INF;
      _max = -VAL_INF;
      for(int step = 0; step < getNumTimeSteps(); step++){
        _steps[step]->setMin(VAL_INF);
        _steps[step]->setMax(-VAL_INF);
        if(_type == NodeData || _type == ElementData){
          // treat these 2 special cases separately for maximum efficiency
          int numComp = _steps[step]->getNumComponents();
          for(int i = 0; i < _steps[step]->getNumData(); i++){
            double *d = _steps[step]->getData(i);
            if(d){
              double val = ComputeScalarRep(numComp, d);
              _steps[step]->setMin(std::min(_steps[step]->getMin(), val));
              _steps[step]->setMax(std::max(_steps[step]->getMax(), val));
            }
          }
        }
        else{
          // general case (slower)
          for(int ent = 0; ent < getNumEntities(step); ent++){
            for(int ele = 0; ele < getNumElements(step, ent); ele++){
              if(skipElement(step, ent, ele)) continue;
              for(int nod = 0; nod < getNumNodes(step, ent, ele); nod++){
                double val;
                getScalarValue(step, ent, ele, nod, val);
                _steps[step]->setMin(std::min(_steps[step]->getMin(), val));
                _steps[step]->setMax(std::max(_steps[step]->getMax(), val));
              }
            }
          }
        }
        _min = std::min(_min, _steps[step]->getMin());
        _max = std::max(_max, _steps[step]->getMax());
      }
    
      // add interpolation data for known element types (this might be
      // overidden later)
      for(int step = 0; step < getNumTimeSteps(); step++){
        GModel *m = _steps[step]->getModel();
        for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
          if((*it)->lines.size()) 
            _addInterpolationMatricesForElement((*it)->lines[0]);
        }
        for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
          if((*it)->triangles.size()) 
            _addInterpolationMatricesForElement((*it)->triangles[0]);
          if((*it)->quadrangles.size()) 
            _addInterpolationMatricesForElement((*it)->quadrangles[0]);
        }
        for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
          if((*it)->tetrahedra.size()) 
            _addInterpolationMatricesForElement((*it)->tetrahedra[0]);
          if((*it)->hexahedra.size()) 
            _addInterpolationMatricesForElement((*it)->hexahedra[0]);
          if((*it)->prisms.size()) 
            _addInterpolationMatricesForElement((*it)->prisms[0]);
          if((*it)->pyramids.size()) 
            _addInterpolationMatricesForElement((*it)->pyramids[0]);
        }
      }
    
      return PViewData::finalize();
    }
    
    void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
    {
      int type = e->getType();
      const functionSpace *fs = e->getFunctionSpace();
      if(fs){
        if(e->getPolynomialOrder() > 1)
          setInterpolationMatrices(type, fs->coefficients, fs->monomials,
                                   fs->coefficients, fs->monomials);
        else
          setInterpolationMatrices(type, fs->coefficients, fs->monomials);
      }                               
    }
    
    MElement *PViewDataGModel::_getElement(int step, int ent, int ele)
    {
      static int lastStep = -1, lastEnt = -1, lastEle = -1;
      static MElement *curr = 0;
      if(step != lastStep || ent != lastEnt || ele != lastEle)
        curr = _steps[step]->getEntity(ent)->getMeshElement(ele);
      return curr;
    }
    
    int PViewDataGModel::getNumTimeSteps()
    {
      return _steps.size();
    }
    
    double PViewDataGModel::getTime(int step)
    {
      if(_steps.empty()) return 0.;
      return _steps[step]->getTime();
    }
    
    double PViewDataGModel::getMin(int step)
    {
      if(step < 0 || _steps.empty()) return _min;
      return _steps[step]->getMin();
    }
    
    double PViewDataGModel::getMax(int step)
    {
      if(step < 0 || _steps.empty()) return _max;
      return _steps[step]->getMax();
    }
    
    SBoundingBox3d PViewDataGModel::getBoundingBox(int step)
    { 
      if(step < 0 || _steps.empty()){
        SBoundingBox3d tmp;
        for(unsigned int i = 0; i < _steps.size(); i++)
          tmp += _steps[i]->getBoundingBox();
        return tmp;
      }
      return _steps[step]->getBoundingBox();
    }
    
    int PViewDataGModel::getNumScalars(int step)
    {
      if(_steps.empty()) return 0;
      // to generalize
      if(_steps[0]->getNumComponents() == 1) return getNumElements(0);
      return 0;
    }
    
    int PViewDataGModel::getNumVectors(int step)
    {
      if(_steps.empty()) return 0;
      // to generalize
      if(_steps[0]->getNumComponents() == 3) return getNumElements(0);
      return 0;
    }
    
    int PViewDataGModel::getNumTensors(int step)
    {
      if(_steps.empty()) return 0;
      // to generalize
      if(_steps[0]->getNumComponents() == 9) return getNumElements(0);
      return 0;
    }
    
    int PViewDataGModel::getNumPoints(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it)
        n += (*it)->points.size();
      return n;
    }
    
    int PViewDataGModel::getNumLines(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
        n += (*it)->lines.size();
      return n;
    }
    
    int PViewDataGModel::getNumTriangles(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
        n += (*it)->triangles.size();
      return n;
    }
    
    int PViewDataGModel::getNumQuadrangles(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
        n += (*it)->quadrangles.size();
      return n;
    }
    
    int PViewDataGModel::getNumTetrahedra(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        n += (*it)->tetrahedra.size();
      return n;
    }
    
    int PViewDataGModel::getNumHexahedra(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        n += (*it)->hexahedra.size();
      return n;
    }
    
    int PViewDataGModel::getNumPrisms(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        n += (*it)->prisms.size();
      return n;
    }
    
    int PViewDataGModel::getNumPyramids(int step)
    {
      if(_steps.empty()) return 0;
      GModel *m = _steps[0]->getModel(); // to generalize
      int n = 0;
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        n += (*it)->pyramids.size();
      return n;
    }
    
    int PViewDataGModel::getNumEntities(int step)
    {
      if(_steps.empty()) return 0;
      // to generalize
      if(step < 0) return _steps[0]->getNumEntities();
      return _steps[step]->getNumEntities();
    }
    
    int PViewDataGModel::getNumElements(int step, int ent)
    {
      if(_steps.empty()) return 0;
      // to generalize
      if(step < 0 && ent < 0) return _steps[0]->getModel()->getNumMeshElements();
      if(step < 0) return _steps[0]->getEntity(ent)->getNumMeshElements();
      if(ent < 0) return _steps[step]->getModel()->getNumMeshElements(); 
      return _steps[step]->getEntity(ent)->getNumMeshElements();
    }
    
    int PViewDataGModel::getDimension(int step, int ent, int ele)
    {
      return _getElement(step, ent, ele)->getDim();
    }
    
    int PViewDataGModel::getNumNodes(int step, int ent, int ele)
    {
      MElement *e = _getElement(step, ent, ele);
      if(_type == GaussPointData){
        return _steps[step]->getGaussPoints(e->getTypeForMSH()).size() / 3;
      }
      else{
        if(isAdaptive())
          return e->getNumVertices();
        else
          return e->getNumPrimaryVertices();
      }
    }
    
    int PViewDataGModel::getNode(int step, int ent, int ele, int nod, 
                                 double &x, double &y, double &z)
    {
      MElement *e = _getElement(step, ent, ele);
      if(_type == GaussPointData){ 
        std::vector<double> &p(_steps[step]->getGaussPoints(e->getTypeForMSH()));
        if(p[0] == 1.e22){
          // hack: the points are the element vertices
          x = e->getVertex(nod)->x();
          y = e->getVertex(nod)->y();
          z = e->getVertex(nod)->z();
        }
        else{
          double vx[8], vy[8], vz[8];
          for(int i = 0; i < e->getNumPrimaryVertices(); i++){
            vx[i] = e->getVertex(i)->x();
            vy[i] = e->getVertex(i)->y();
            vz[i] = e->getVertex(i)->z();
          }
          x = e->interpolate(vx, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2], 1, 1);
          y = e->interpolate(vy, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2], 1, 1);
          z = e->interpolate(vz, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2], 1, 1);
        }
        return 0;
      }
      else{
        MVertex *v = e->getVertex(nod);
        x = v->x();
        y = v->y();
        z = v->z();
        return v->getIndex();
      }
    }
    
    void PViewDataGModel::setNode(int step, int ent, int ele, int nod, 
                                  double x, double y, double z)
    {
      MVertex *v = _getElement(step, ent, ele)->getVertex(nod);
      v->x() = x;
      v->y() = y;
      v->z() = z;
    }
    
    void PViewDataGModel::tagNode(int step, int ent, int ele, int nod, int tag)
    {
      MVertex *v = _getElement(step, ent, ele)->getVertex(nod);
      v->setIndex(tag);
    }
    
    int PViewDataGModel::getNumComponents(int step, int ent, int ele)
    {
      return _steps[step]->getNumComponents();
    }
    
    int PViewDataGModel::getNumValues(int step, int ent, int ele)
    {
      if(_type == ElementNodeData || _type == NodeData){
        return getNumNodes(step, ent, ele) * getNumComponents(step, ent, ele);
      }
      else if(_type == ElementData){
        return getNumComponents(step, ent, ele);
      }
      else{
        Msg::Error("getNumValue() should not be used on this type of view");
        return getNumComponents(step, ent, ele);
      }
    }
    
    void PViewDataGModel::getValue(int step, int ent, int ele, int idx, double &val)
    {
      MElement *e = _getElement(step, ent, ele);
      if(_type == ElementNodeData || _type == ElementData){
        val = _steps[step]->getData(e->getNum())[idx];
      }
      else if(_type == NodeData){
        int numcomp = _steps[step]->getNumComponents();
        int nod = idx / numcomp;
        int comp = idx % numcomp;
        val = _steps[step]->getData(e->getVertex(nod)->getNum())[comp];
      }
      else{
        Msg::Error("getValue(index) should not be used on this type of view");
      }
    }
    
    void PViewDataGModel::getValue(int step, int ent, int ele, int nod, int comp, double &val)
    {
      MElement *e = _getElement(step, ent, ele);
      switch(_type){
      case NodeData: 
        val = _steps[step]->getData(e->getVertex(nod)->getNum())[comp];
        break;
      case ElementNodeData:
      case GaussPointData:
        val = _steps[step]->getData(e->getNum())[_steps[step]->getNumComponents() * nod + comp];
        break;
      case ElementData:
      default: 
        val = _steps[step]->getData(e->getNum())[comp];
        break;
      }
    }
    
    void PViewDataGModel::setValue(int step, int ent, int ele, int nod, int comp, double val)
    {
      MElement *e = _getElement(step, ent, ele);
      switch(_type){
      case NodeData: 
        _steps[step]->getData(e->getVertex(nod)->getNum())[comp] = val;
        break;
      case ElementNodeData:
      case GaussPointData: 
        _steps[step]->getData(e->getNum())[_steps[step]->getNumComponents() * nod + comp] = val;
        break;
      case ElementData: 
      default: 
        _steps[step]->getData(e->getNum())[comp] = val;
        break;
      }
    }
    
    int PViewDataGModel::getNumEdges(int step, int ent, int ele)
    { 
      return _getElement(step, ent, ele)->getNumEdges();
    }
    
    int PViewDataGModel::getType(int step, int ent, int ele)
    {
      return _getElement(step, ent, ele)->getType();
    }
    
    void PViewDataGModel::revertElement(int step, int ent, int ele)
    {
      if(!step) _getElement(step, ent, ele)->revert();
    }
    
    void PViewDataGModel::smooth()
    {
      if(_type == NodeData || _type == GaussPointData) return;
      std::vector<stepData<double>*> _steps2;
      for(unsigned int step = 0; step < _steps.size(); step++){
        GModel *m = _steps[step]->getModel();
        int numComp = _steps[step]->getNumComponents();
        _steps2.push_back(new stepData<double>(m, numComp, _steps[step]->getFileName(),
                                               _steps[step]->getFileIndex()));
        std::map<int, int> nodeConnect;
        for(int ent = 0; ent < getNumEntities(step); ent++){
          for(int ele = 0; ele < getNumElements(step, ent); ele++){
            MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
            double val;
            if(!getValueByIndex(step, e->getNum(), 0, 0, val)) continue;
            for(int nod = 0; nod < e->getNumVertices(); nod++){
              MVertex *v = e->getVertex(nod);
              if(nodeConnect.count(v->getNum()))
                nodeConnect[v->getNum()]++;
              else
                nodeConnect[v->getNum()] = 1;
              double *d = _steps2.back()->getData(v->getNum(), true);
              for(int j = 0; j < numComp; j++)
                if(getValueByIndex(step, e->getNum(), nod, j, val)) d[j] += val;
            }
          }
        }
        for(int i = 0; i < _steps2.back()->getNumData(); i++){
          double *d = _steps2.back()->getData(i);
          if(d){
            double f = nodeConnect[i];
            if(f) for(int j = 0; j < numComp; j++) d[j] /= f;
          }
        }
      }
      for(unsigned int i = 0; i < _steps.size(); i++) delete _steps[i];
      _steps = _steps2;
      _type = NodeData;
      finalize();
    }
    
    bool PViewDataGModel::skipEntity(int step, int ent)
    {
      if(step >= getNumTimeSteps()) return true;
      return !_steps[step]->getEntity(ent)->getVisibility();
    }
    
    bool PViewDataGModel::skipElement(int step, int ent, int ele, bool checkVisibility)
    {
      if(step >= getNumTimeSteps()) return true;
      stepData<double> *sd = _steps[step];
      if(!_steps[step]->getNumData()) return true;
      MElement *e = _getElement(step, ent, ele);
      if(checkVisibility && !e->getVisibility()) return true;
      if(_type == NodeData){
        for(int i = 0; i < e->getNumVertices(); i++)
          if(!sd->getData(e->getVertex(i)->getNum())) return true;
      }
      else{
        if(!sd->getData(e->getNum())) return true;
      }
      return false;
    }
    
    bool PViewDataGModel::hasTimeStep(int step)
    {
      if(step < getNumTimeSteps() && _steps[step]->getNumData()) return true;
      return false;
    }
    
    bool PViewDataGModel::hasPartition(int part)
    {
      return _partitions.find(part) != _partitions.end();
    }
    
    bool PViewDataGModel::hasMultipleMeshes()
    {
      if(_steps.size() <= 1) return false;
      GModel *m = _steps[0]->getModel();
      for(unsigned int i = 1; i < _steps.size(); i++)
        if(m != _steps[i]->getModel()) return true;
      return false;
    }
    
    bool PViewDataGModel::hasModel(GModel *model, int step)
    {
      if(step < 0){
        for(unsigned int i = 0; i < _steps.size(); i++)
          if(model == _steps[i]->getModel()) return true;
        return false;
      }
      return (model == _steps[step]->getModel());
    }
    
    GEntity *PViewDataGModel::getEntity(int step, int ent)
    {
      return _steps[step]->getEntity(ent);
    }
    
    bool PViewDataGModel::getValueByIndex(int step, int dataIndex, int nod, int comp, double &val)
    {
      double *d = _steps[step]->getData(dataIndex);
      if(!d) return false;
    
      if(_type == NodeData || _type == ElementData)
        val = d[comp];
      else
        val = d[_steps[step]->getNumComponents() * nod + comp];
      return true;
    }