V3VisualGrid.h

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding, Matthias Kraatz
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_V3VISUALGRID_H
#define AMROC_V3VISUALGRID_H

#include "VisualGrid.h" 
#include "HDFToV3/V3VisualGridBase.h"
#include "HDFToV3/V3VisualGridData.h"

#define NSURF 6

/* strcpn copies str2 into str1 up to, but not including the first null */

int strcpn(char* str1, char* str2) {
  register int i = 0;
  while (str2[i]!='\0') {
    str1[i] = str2[i];
    i++;
  }  
  return i;
}

/* strcpb copies str2 into str1 up to, but not including the first null and
   then pads with blanks to the len character */

int strcpb( char* str1, char* str2, int len ) {
  register int i = 0, j;
  while (str2[i] != '\0') {
    str1[i] = str2[i];
    i++;
  }
  for (j=i; j < len; j++) str1[j] = ' ';  
  return i;
}


template <class VectorType>
class V3VisualGrid : public V3VisualGridBase, public VisualGrid<VectorType> {
  typedef typename VectorType::InternalDataType DataType;
  typedef VisualGrid<VectorType> base;

public:
  typedef typename base::data_block_type database_block_type;
  typedef V3VisualGridData<VectorType> data_block_type;
  typedef typename base::block_list_type blockbase_list_type;
  typedef std::vector<data_block_type*> block_list_type;
  typedef typename base::evaluator_type evaluator_type;
  typedef typename base::hierarchy_storage_type hierarchy_storage_type;
  typedef typename base::hdata_block_type hdata_block_type;

  typedef std::pair<int, int> node_pair_type;
  typedef std::multimap<int, int> node_pair_list_type;

public:
  V3VisualGrid(evaluator_type* ev) : V3VisualGridBase(ev), base(ev), 
                                     _Equiv(0), _PlotCrossSection(0) {
    for (int s=0; s<NSURF; s++) 
      SurfaceDisp[s] = 1;
    base::_BlockList = (blockbase_list_type *) new block_list_type;
    base::Make3d = true;
  }

  virtual ~V3VisualGrid() {
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++)
      delete (*bit);
    BlockList().clear();      
    delete base::_BlockList;
  }

  virtual void update() {    
    std::string iname;
    if (UpdateName()) iname = *UpdateName();
    else iname = "display.in";

    // Generating display.in if file is missing; 
    std::ifstream infile;
    std::ofstream outfile;
    infile.open(iname.c_str(), std::ios::in);
    if (!infile) {
      if (!UpdateName()) {
        outfile.open(iname.c_str(), std::ios::out);
        outfile << "Type  6\nUseEquiv 1\n\nDisplayMinLevel    0\nDisplayMaxLevel    10\n\n";
        outfile << "ScriptBegin {\n        movewindow 7,2000,0\n";
        outfile << "        Edgesofftransparenton\n}\n";
        outfile.close();
      }
    }
    else{
      infile.close();
    }

    ControlDevice DisplayCtrl(GetFileControlDevice(iname.c_str(),""));
    RegisterAt(DisplayCtrl,"Type",base::_FKeys); 
    RegisterAt(DisplayCtrl,"UseEquiv",_Equiv);
    RegisterAt(DisplayCtrl,"DisplayMinLevel",base::_MinLevel); 
    RegisterAt(DisplayCtrl,"DisplayMaxLevel",base::_MaxLevel);
    RegisterAt(DisplayCtrl,"CutOut",base::_Eliminate);
    RegisterAt(DisplayCtrl,"CutOutValue",base::_EliminateValue);
    delete base::PlotLevel;
    base::PlotLevel = new int[NLevels()];
    for (int lev=0; lev<NLevels(); lev++) {
      char LevelName[16];      
      std::sprintf(LevelName,"PlotGrid(%d)",lev);
      RegisterAt(DisplayCtrl,LevelName,base::PlotLevel[lev]);
      base::PlotLevel[lev] = 0;
    }    
    int d;
    for (d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"ShowMin(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::show[2*d]);
      std::sprintf(VariableName,"ShowMax(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::show[2*d+1]);     
      std::sprintf(VariableName,"ShowCut(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::cut[d]);     
      std::sprintf(VariableName,"Symmetry(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::Symmetry[d]);     
      std::sprintf(VariableName,"Periodic(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::Periodic[d]);     
    }   
    RegisterAt(DisplayCtrl,"PlotCrossSection",_PlotCrossSection);
    for (d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"CrossPlaneP1(%d)",d);
      RegisterAt(DisplayCtrl,VariableName,CrossPlane[d]);
    }
    for (d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"CrossPlaneP2(%d)",d);
      RegisterAt(DisplayCtrl,VariableName,CrossPlane[3+d]);
    }
    for (d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"CrossPlaneP3(%d)",d);
      RegisterAt(DisplayCtrl,VariableName,CrossPlane[6+d]);
    }
    TheEvaluator().register_at(DisplayCtrl,"");
    DisplayCtrl.update();  

    TheEvaluator().update();
  }

  virtual int SetUp(int step, std::string** filenames, int Nfilenames) {
    if (base::SetUp(step,filenames,Nfilenames) >= 0) 
      for (typename block_list_type::iterator bit=BlockList().begin();
           bit!=BlockList().end(); bit++)
        (*bit)->SetCrossPlane(CrossPlane);
    return 0;
  }

  virtual void SetGrid(float xyz[][3]) {
    SetGridNodes(xyz);
  }

  virtual void SetScal(int *key,float v[]) {
    if (FKeys() == 1)
      SetScalNodes(key, v);
    else
      SetScalCells(key, v);
    if (TheEvaluator().mvals[2*(*key-1)] ==
        TheEvaluator().mvals[2*(*key-1)+1]) {
      if (TheEvaluator().mvals[2*(*key-1)] == 0.0) {
        TheEvaluator().mvals[2*(*key-1)]   = -0.1;
        TheEvaluator().mvals[2*(*key-1)+1] =  0.1; 
      }
      else {
        TheEvaluator().mvals[2*(*key-1)]   -=
          0.1*TheEvaluator().mvals[2*(*key-1)];
        TheEvaluator().mvals[2*(*key-1)+1] +=
          0.1*TheEvaluator().mvals[2*(*key-1)+1];
      }
    }    
    
    int state;
    if (FKeys() == 1) state = 5;
    else state = 9;        
    long int kk=(long int)(*key); 
    V3_SETSTATE(&state, &kk, base::TheEvaluator().Mvals(*key-1), 
                base::TheEvaluator().TitleKeys(*key-1), 
                std::strlen(TheEvaluator().TitleKeys(*key-1)));
  }

  virtual void Vect(int *key,float v[][3]) {
    if (FKeys() == 1)
      VectNodes(key, v);
    else
      VectCells(key, v);
  }

  virtual void FindNodePairs() {
    register int i,j,k;
    int kn,knmatch,ijk[3];
    std::cerr << std::endl << "Finding matching node pairs";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++){
       
      BBox bbnodes(growupper((*bit)->DB().bbox(),1));
      
      std::cerr << ".";
      for (typename block_list_type::iterator bitb=bit+1;
           bitb!=BlockList().end(); bitb++) { 
        
        BBox bbnodes_cur(growupper((*bitb)->DB().bbox(),1));
        BBox bbint = bbnodes * bbnodes_cur;
        if (bbint.empty())
          continue;
        bbint.lower() = (bbint.lower()/bbint.stepsize())*bbint.stepsize();
        bbint.upper() = (bbint.upper()/bbint.stepsize())*bbint.stepsize();

        //starting with front surface....;
        if (bbint.lower(2) == bbnodes.lower(2)) {
          ijk[2]=bbint.lower(2);
          for(j=bbint.lower(1); j<=bbint.upper(1); j+=bbint.stepsize(1)) {
            ijk[1]=j;
            for(i=bbint.lower(0); i<=bbint.upper(0); i+=bbint.stepsize(0)) {
              ijk[0]=i;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end()) 
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }
          }
        }

        //doing surface back ...;
        if (bbint.upper(2) == bbnodes.upper(2)) {
          ijk[2]=bbint.upper(2);
          for(j=bbint.lower(1); j<=bbint.upper(1); j+=bbint.stepsize(1)) {
            ijk[1]=j;
            for(i=bbint.lower(0); i<=bbint.upper(0); i+=bbint.stepsize(0)) {
              ijk[0]=i;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end())
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }
          }
        }

        //doing surface left ...;
        if (bbint.lower(0) == bbnodes.lower(0)) {
          ijk[0]=bbint.lower(0);
          for(k=bbint.lower(2); k<=bbint.upper(2); k+=bbint.stepsize(2)) {
            ijk[2]=k;
            for(j=bbint.lower(1); j<=bbint.upper(1); j+=bbint.stepsize(1)) {
              ijk[1]=j;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end())
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }
          }
        }

        //doing surface right ...;
        if (bbint.upper(0) == bbnodes.upper(0)) {
          ijk[0]=bbint.upper(0);
          for(k=bbint.lower(2); k<=bbint.upper(2); k+=bbint.stepsize(2)) {
            ijk[2]=k;
            for(j=bbint.lower(1); j<=bbint.upper(1); j+=bbint.stepsize(1)) {
              ijk[1]=j;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end())
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }     
          }
        }
        
        //doing surface down ...;
        if (bbint.lower(1) == bbnodes.lower(1)) {
          ijk[1]=bbint.lower(1);
          for(k=bbint.lower(2); k<=bbint.upper(2); k+=bbint.stepsize(2)) {
            ijk[2]=k;
            for(i=bbint.lower(0); i<=bbint.upper(0); i+=bbint.stepsize(0)) {
              ijk[0]=i;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end())
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }
          }
        }

        //doing surface up ...;
        if (bbint.upper(1) == bbnodes.upper(1)) {
          ijk[1]=bbint.upper(1);
          for(k=bbint.lower(2); k<=bbint.upper(2); k+=bbint.stepsize(2)) {
            ijk[2]=k;
            for(i=bbint.lower(0); i<=bbint.upper(0); i+=bbint.stepsize(0)) {
              ijk[0]=i;
              kn = (*bit)->FindMatchingNode(ijk);          
              knmatch = (*bitb)->FindMatchingNode(ijk);
              if(kn<knmatch) 
                if (NodePairs.empty() ||
                    NodePairs.find(knmatch)==NodePairs.end())
                  NodePairs.insert(node_pair_type(knmatch, kn));
            }
          }       
        }
      }
    }
    std::cerr << std::endl;
  }
    
  virtual void SetNodePairs(int listequiv[][2]){
    int i = 0;
    for (node_pair_list_type::iterator npit=NodePairs.begin();
         npit!=NodePairs.end(); npit++, i++) {
      listequiv[i][0]=(*npit).second;
      listequiv[i][1]=(*npit).first;
    }
    NodePairs.clear();
  }
       
  virtual int NFacetsOnSurface() {    
    int dummy[1][4];
    int ksurf=0;
    for (int kn=0; kn<6; kn++) {
      int d=kn/2;
      BBoxList surflist;
      typename block_list_type::iterator bit;
      for (bit=BlockList().begin();bit!=BlockList().end(); bit++) {
        BBox bbnodes(growupper((*bit)->DB().bbox(),1));
        if (kn%2) bbnodes.upper(d) = bbnodes.lower(d);
        else bbnodes.lower(d) = bbnodes.upper(d); 
        surflist.add(bbnodes);
      }
      for (bit=BlockList().begin();bit!=BlockList().end(); bit++)
        (*bit)->SetSurfaces(kn, surflist, ksurf, dummy, 0);
    }
    return ksurf;
  }
  
  virtual void SetSurfaces(int nsurf[][2], int surf[], int scel[][4],
                           char tsurf[][20]) {
    int ksurf = 0;
    for (int kn=0; kn<6; kn++) {
      int d=kn/2;
      BBoxList surflist;
      typename block_list_type::iterator bit;
      for (bit=BlockList().begin();bit!=BlockList().end(); bit++) {
        BBox bbnodes(growupper((*bit)->DB().bbox(),1));
        if (kn%2) bbnodes.upper(d) = bbnodes.lower(d);
        else bbnodes.lower(d) = bbnodes.upper(d); 
        surflist.add(bbnodes);
      }
      for (bit=BlockList().begin();bit!=BlockList().end(); bit++)
        (*bit)->SetSurfaces(kn, surflist, ksurf, scel, 1);
      nsurf[kn][0] = ksurf;
      nsurf[kn][1] = SurfaceDisp[kn];
    }
    strcpb(tsurf[0],"West",20);
    strcpb(tsurf[1],"East",20);   
    strcpb(tsurf[2],"Bottom",20);
    strcpb(tsurf[3],"Top",20);
    strcpb(tsurf[4],"Front",20);
    strcpb(tsurf[5],"Back",20);

    std::cerr << "Actual ksurf: "  << ksurf << std::endl;
  }

  virtual void Draw3D() {
    DrawFrameAroundGrid(); 
    Draw3DGrids();
    Draw3dCrossSections();
    //    DrawOrigin(); 
    //    DrawXYZ();      
  }      

  void Draw3DGrids() {
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      if (base::PlotLevel[(*bit)->Level()] > 0)
        (*bit)->Draw3DGrids();
  }

  void Draw3dCrossSections() {
    if(_PlotCrossSection>0){
      for (typename block_list_type::iterator bit=BlockList().begin();
           bit!=BlockList().end(); bit++) 
        if (base::PlotLevel[(*bit)->Level()] != 0)
          (*bit)->Draw3DCrossSections();
    }
  }
  
  void DrawOrigin() {
    int itype = 3;
    int icolor = 1;
    int np = 1*2;
    float radii = 0.;
    int intzero = 0;
    float linesx[2*3] = { 5., 0., 0., 0., 0., 0.};
    GeomMapField(linesx, 2*3);
    float colorx[2] = { 0., 0.,};
    V3_OBJECT3D(&itype, &icolor, linesx, &radii, colorx, &intzero, &np);

    float linesy[2*3] = { 0., 5., 0., 0., 0., 0.};
    GeomMapField(linesy, 2*3);
    float colory[2] = {0.5, 0.5};
    V3_OBJECT3D(&itype, &icolor, linesy, &radii, colory, &intzero, &np);
    
    float linesz[2*3] = { 0., 0., 5., 0., 0., 0.};
    GeomMapField(linesz, 2*3);
    float colorz[2] = { 1., 1.};
    V3_OBJECT3D(&itype, &icolor, linesz, &radii, colorz, &intzero, &np);
    
    float floatzero = 0.;
    itype = 1;
    np = 2*3;
    float trianglesx[2*3*3] = { 6., 0., 0., 5., .5, 0., 5.,-.5, 0.,
                                  6., 0., 0., 5., 0., .5, 5., 0.,-.5};
    GeomMapField(trianglesx, 2*3*3);
    float tricolorx[2*3] = { 0., 0., 0., 0., 0., 0.};
    V3_OBJECT3D(&itype, &icolor, trianglesx, &radii, tricolorx, &intzero, &np);
    
    float trianglesy[2*3*3] = { 0., 6., 0., .5, 5., 0.,-.5, 5., 0.,
                                  0., 6., 0., 0., 5., .5, 0., 5.,-.5};
    GeomMapField(trianglesy, 2*3*3);
    float tricolory[2*3] = { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
    V3_OBJECT3D(&itype, &icolor, trianglesy, &radii, tricolory, &intzero, &np);
      
    float trianglesz[2*3*3] = { 0., 0., 6., 0., .5, 5., 0.,-.5, 5.,
                                  0., 0., 6., .5, 0., 5.,-.5, 0., 5.};
    GeomMapField(trianglesz, 2*3*3);
    float tricolorz[2*3] = { 1., 1., 1., 1., 1., 1.};     
    V3_OBJECT3D(&itype, &icolor, trianglesz, &radii, tricolorz, &intzero, &np);
      
    icolor = 0;
    itype = 5;
    np = 1;
    float midpoint[1*3] = { 0., 0., 0.};
    radii = .5*base::dx[0];
    V3_OBJECT3D(&itype, &icolor, midpoint, &radii, &floatzero, &intzero, &np);
  }

  void DrawXYZ(){
    int itype = 0;
    int icolor = 1;
    int np = 2*4;
    float radii = 0.;
    int intzero = 0;
    float slashx[2*4*3] = { 6,-1,-1,6.5,-1,-1,8,-4,-1,7.5,-4,-1,
                           7.5,-1,-1,8,-1,-1,6.5,-4,-1,6,-4,-1};
    GeomMapField(slashx, 2*4*3);
    float colorx[2*4] = { 0,0,0,0,0,0,0,0};
    V3_OBJECT3D(&itype, &icolor, slashx, &radii, colorx, &intzero, &np);
      
    float slashy[2*4*3] = { -1.5,8,-1,-1,8,-1,-2.5,5,-1,-3,5,-1,
                           -3,8,-1,-2.5,8,-1,-2,6.5,-1,-2.42,6,-1};
      
    GeomMapField(slashy, 2*4*3);
    float colory[2*4] = { 0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5};
    V3_OBJECT3D(&itype, &icolor, slashy, &radii, colory, &intzero, &np);

    np=3*4;
    float slashz[3*4*3] = { -3,-1,7,-1,-1,7,-1,-1.5,7,-3,-1.5,7,
                           -1.5,-1,7,-1,-1,7,-2.5,-4,7,-3,-4,7,
                           -3,-3.5,7,-1,-3.5,7,-1,-4,7,-3,-4,7};
      
    GeomMapField(slashz, 3*4*3);
    float colorz[3*4] = {1,1,1,1,1,1,1,1,1,1,1,1};
    V3_OBJECT3D(&itype, &icolor, slashz, &radii, colorz, &intzero, &np);
  }

  void DrawFrameAroundGrid() {
    float radii = 0.;
    float floatzero = 0.;
    int intzero = 0;
    int itype = 4;
    int icolor = 0;
    int np = 5;
    float line1[5*3] = { base::geom[0], base::geom[2], base::geom[4], 
                         base::geom[0], base::geom[3], base::geom[4],
                         base::geom[1], base::geom[3], base::geom[4], 
                         base::geom[1], base::geom[2], base::geom[4],
                         base::geom[0], base::geom[2], base::geom[4] };      
    V3_OBJECT3D(&itype, &icolor, line1, &radii, &floatzero, &intzero, &np);

    float line2[5*3] = { base::geom[0], base::geom[2], base::geom[5], 
                         base::geom[0], base::geom[3], base::geom[5],
                         base::geom[1], base::geom[3], base::geom[5], 
                         base::geom[1], base::geom[2], base::geom[5],
                         base::geom[0], base::geom[2], base::geom[5] };      
    V3_OBJECT3D(&itype, &icolor, line2, &radii, &floatzero, &intzero, &np);

    itype = 3;
    np = 2*4;
    float lines[2*4*3] = { base::geom[0], base::geom[2], base::geom[4], 
                           base::geom[0], base::geom[2], base::geom[5],
                           base::geom[0], base::geom[3], base::geom[4], 
                           base::geom[0], base::geom[3], base::geom[5],
                           base::geom[1], base::geom[2], base::geom[4], 
                           base::geom[1], base::geom[2], base::geom[5],
                           base::geom[1], base::geom[3], base::geom[4], 
                           base::geom[1], base::geom[3], base::geom[5]};        
    V3_OBJECT3D(&itype, &icolor, lines, &radii, &floatzero, &intzero, &np);
  }    
  
  virtual const int& Equiv() const { return _Equiv; }
  virtual int Kequiv() const { return NodePairs.size(); }

  virtual void BlockListAppend(hdata_block_type& workdata, BBox& bb, int& comp, int& Level,
                               int& proc, float* dx, float* geom) {
    BlockList().push_back(new data_block_type(*this, workdata, bb, comp,
                                              TheEvaluator(), Level, proc, dx, geom));
  }

  block_list_type& BlockList() { return (block_list_type &)*base::_BlockList; }
  block_list_type& BlockList() const { return (block_list_type &)*base::_BlockList; }
 
  virtual void register_at(ControlDevice& Ctrl,const std::string& prefix) {
    base::register_at(Ctrl,prefix); }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, ""); }
  virtual int NCells() { return base::NCells(); }
  virtual int NNodes() { return base::NNodes(); }
  virtual void SetBlocks(int blocks[]) { base::SetBlocks(blocks); }
  virtual void SetGridCells(float xyz[][3]) {
    base::SetGridCells(xyz); }
  virtual void SetGridNodes(float xyz[][3]) {
    base::SetGridNodes(xyz); }
  virtual void SetGridConns(int con[][8]) {
    base::SetGridConns(con); }
  virtual void SetScalCells(int *key,float v[]) {
    base::SetScalCells(key,v); }
  virtual void SetScalNodes(int *key,float v[]) {
    base::SetScalNodes(key,v); }
  virtual void VectNodes(int *key,float v[][3]) {
    base::VectNodes(key,v); }
  virtual void VectCells(int *key,float v[][3]) {
    base::VectCells(key,v); }
  virtual int Size() const { return base::Size(); }
  virtual const int& MinLevel() const { return base::MinLevel(); }
  virtual const int& MaxLevel() const { return base::MaxLevel(); }
  virtual const int& NLevels() const { return base::NLevels(); }
  virtual const int& RefineFactor(const int l) const {
    return base::RefineFactor(l); }
  virtual const BBox& GlobalBBox() const {
    return base::GlobalBBox(); }
  virtual const int& FKeys() const { return base::FKeys(); }
  virtual evaluator_type& TheEvaluator() const {
    return base::TheEvaluator(); }
  virtual evaluator_type& TheEvaluator() {
    return base::TheEvaluator(); }
  virtual float64 RealTime() const { return base::RealTime(); }
  virtual void SetUpdateName(std::string* iname) {
    base::SetUpdateName(iname); }
  virtual std::string* UpdateName() { return base::UpdateName(); }
  virtual std::string* UpdateName() const {
    return base::UpdateName(); }

private:
  void GeomMapField(float field[], int size) {
    for (int i=0; i<size; i+=3) {
      field[i] *= base::dx[0]; field[i+1] *= base::dx[1]; field[i+2] *= base::dx[2];
      field[i] += base::geom[0]; field[i+1] += base::geom[2]; field[i+2] += base::geom[4]; 
    }
  }

private:
  int _Equiv;
  node_pair_list_type NodePairs;
  int SurfaceDisp[NSURF];
  float CrossPlane[9];
  int _PlotCrossSection;
};

#endif

---