浏览代码

Fix MPI_Tree::Write2File(...) for particles

Dhairya Malhotra 10 年之前
父节点
当前提交
6157c25c2d
共有 9 个文件被更改,包括 364 次插入128 次删除
  1. 2 1
      include/cheb_node.hpp
  2. 69 49
      include/cheb_node.txx
  3. 6 0
      include/fmm_node.hpp
  4. 110 0
      include/fmm_node.txx
  5. 2 1
      include/mpi_node.hpp
  6. 122 42
      include/mpi_node.txx
  7. 16 0
      include/mpi_tree.hpp
  8. 36 34
      include/mpi_tree.txx
  9. 1 1
      include/pvfmm.hpp

+ 2 - 1
include/cheb_node.hpp

@@ -176,7 +176,8 @@ class Cheb_Node: public MPI_Node<Real_t>{
   /**
    * \brief Append node VTU data to vectors.
    */
-  virtual void VTU_Data(std::vector<Real_t>& coord, std::vector<Real_t>& value, std::vector<int32_t>& connect, std::vector<int32_t>& offset, std::vector<uint8_t>& types, int lod=-1);
+  template <class VTUData_t, class Node_t>
+  static void VTU_Data(VTUData_t& vtu_data, std::vector<Node_t*>& nodes, int lod);
 
   /**
    * \brief Compute gradient of the data.

+ 69 - 49
include/cheb_node.txx

@@ -190,59 +190,79 @@ void Cheb_Node<Real_t>::Unpack(PackedData p0, bool own_data) {
 }
 
 template <class Real_t>
-void Cheb_Node<Real_t>::VTU_Data(std::vector<Real_t>& coord, std::vector<Real_t>& value, std::vector<int32_t>& connect, std::vector<int32_t>& offset, std::vector<uint8_t>& types, int lod){
-  //for(size_t i=0;i<this->pt_coord.Dim();i++) coord.push_back(this->pt_coord[i]);
-  //for(size_t i=0;i<this->pt_value.Dim();i++) value.push_back(this->pt_value[i]);
+template <class VTUData_t, class Node_t>
+void Cheb_Node<Real_t>::VTU_Data(VTUData_t& vtu_data, std::vector<Node_t*>& nodes, int lod){
+  typedef typename VTUData_t::VTKReal_t VTKReal_t;
+  //MPI_Node<Real_t>::VTU_Data(vtu_data, nodes, lod);
+
+  VTUData_t new_data;
+  { // Set new data
+    new_data.value.resize(1);
+    new_data.name.push_back("cheb_value");
+    std::vector<VTKReal_t>& coord=new_data.coord;
+    std::vector<VTKReal_t>& value=new_data.value[0];
+
+    std::vector<int32_t>& connect=new_data.connect;
+    std::vector<int32_t>& offset =new_data.offset;
+    std::vector<uint8_t>& types  =new_data.types;
+
+    int gridpt_cnt;
+    std::vector<Real_t> grid_pts;
+    {
+      int cheb_deg=lod-1;
+      std::vector<Real_t> cheb_node_;
+      if(cheb_deg>=0) cheb_node_=cheb_nodes<Real_t>(cheb_deg,1);
+      gridpt_cnt=cheb_node_.size()+2;
+      grid_pts.resize(gridpt_cnt);
+      for(size_t i=0;i<cheb_node_.size();i++) grid_pts[i+1]=cheb_node_[i];
+      grid_pts[0]=0.0; grid_pts[gridpt_cnt-1]=1.0;
+    }
 
-  Real_t* c=this->Coord();
-  Real_t s=pow(0.5,this->Depth());
-  size_t point_cnt=coord.size()/COORD_DIM;
-
-  int gridpt_cnt;
-  std::vector<Real_t> grid_pts;
-  {
-    int cheb_deg=lod;
-    std::vector<Real_t> cheb_node_=cheb_nodes<Real_t>(cheb_deg,1);
-    gridpt_cnt=cheb_node_.size()+2;
-    grid_pts.resize(gridpt_cnt);
-    for(size_t i=0;i<cheb_node_.size();i++) grid_pts[i+1]=cheb_node_[i];
-    grid_pts[0]=0.0; grid_pts[gridpt_cnt-1]=1.0;
-  }
-  for(int i0=0;i0<gridpt_cnt;i0++)
-  for(int i1=0;i1<gridpt_cnt;i1++)
-  for(int i2=0;i2<gridpt_cnt;i2++){
-    coord.push_back(c[0]+grid_pts[i2]*s);
-    coord.push_back(c[1]+grid_pts[i1]*s);
-    coord.push_back(c[2]+grid_pts[i0]*s);
-    for(int j=0;j<data_dof;j++) value.push_back(0.0);
-  }
-  for(int i0=0;i0<(gridpt_cnt-1);i0++)
-  for(int i1=0;i1<(gridpt_cnt-1);i1++)
-  for(int i2=0;i2<(gridpt_cnt-1);i2++){
-    for(int j0=0;j0<2;j0++)
-    for(int j1=0;j1<2;j1++)
-    for(int j2=0;j2<2;j2++)
-      connect.push_back(point_cnt + ((i0+j0)*gridpt_cnt + (i1+j1))*gridpt_cnt + (i2+j2)*1);
-    offset.push_back(connect.size());
-    types.push_back(11);
-  }
+    for(size_t nid=0;nid<nodes.size();nid++){
+      Node_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      size_t point_cnt=coord.size()/COORD_DIM;
+      Real_t* c=n->Coord();
+      Real_t s=pow(0.5,n->Depth());
+      for(int i0=0;i0<gridpt_cnt;i0++)
+      for(int i1=0;i1<gridpt_cnt;i1++)
+      for(int i2=0;i2<gridpt_cnt;i2++){
+        coord.push_back(c[0]+grid_pts[i2]*s);
+        coord.push_back(c[1]+grid_pts[i1]*s);
+        coord.push_back(c[2]+grid_pts[i0]*s);
+        for(int j=0;j<n->data_dof;j++) value.push_back(0.0);
+      }
+      for(int i0=0;i0<(gridpt_cnt-1);i0++)
+      for(int i1=0;i1<(gridpt_cnt-1);i1++)
+      for(int i2=0;i2<(gridpt_cnt-1);i2++){
+        for(int j0=0;j0<2;j0++)
+        for(int j1=0;j1<2;j1++)
+        for(int j2=0;j2<2;j2++)
+          connect.push_back(point_cnt + ((i0+j0)*gridpt_cnt + (i1+j1))*gridpt_cnt + (i2+j2)*1);
+        offset.push_back(connect.size());
+        types.push_back(11);
+      }
 
-  {// Set point values.
-    Real_t* val_ptr=&value[point_cnt*data_dof];
-    std::vector<Real_t> x(gridpt_cnt);
-    std::vector<Real_t> y(gridpt_cnt);
-    std::vector<Real_t> z(gridpt_cnt);
-    for(int i=0;i<gridpt_cnt;i++){
-      x[i]=c[0]+s*grid_pts[i];
-      y[i]=c[1]+s*grid_pts[i];
-      z[i]=c[2]+s*grid_pts[i];
+      {// Set point values.
+        Real_t* val_ptr=&value[point_cnt*n->data_dof];
+        std::vector<Real_t> x(gridpt_cnt);
+        std::vector<Real_t> y(gridpt_cnt);
+        std::vector<Real_t> z(gridpt_cnt);
+        for(int i=0;i<gridpt_cnt;i++){
+          x[i]=c[0]+s*grid_pts[i];
+          y[i]=c[1]+s*grid_pts[i];
+          z[i]=c[2]+s*grid_pts[i];
+        }
+        n->ReadVal(x, y, z, val_ptr);
+        //Rearrrange data
+        //(x1,x2,x3,...,y1,y2,...z1,...) => (x1,y1,z1,x2,y2,z2,...)
+        Matrix<Real_t> M(n->data_dof,gridpt_cnt*gridpt_cnt*gridpt_cnt,val_ptr,false);
+        M=M.Transpose();
+      }
     }
-    this->ReadVal(x, y, z, val_ptr);
-    //Rearrrange data
-    //(x1,x2,x3,...,y1,y2,...z1,...) => (x1,y1,z1,x2,y2,z2,...)
-    Matrix<Real_t> M(data_dof,gridpt_cnt*gridpt_cnt*gridpt_cnt,val_ptr,false);
-    M=M.Transpose();
   }
+  AppendVTUData(vtu_data, new_data);
 }
 
 template <class Real_t>

+ 6 - 0
include/fmm_node.hpp

@@ -152,6 +152,12 @@ class FMM_Node: public Node{
    */
   virtual void InitMultipole(PackedData data, bool own_data=true);
 
+  /**
+   * \brief Append node VTU data to vectors.
+   */
+  template <class VTUData_t, class VTUNode_t>
+  static void VTU_Data(VTUData_t& vtu_data, std::vector<VTUNode_t*>& nodes, int lod);
+
   Vector<Real_t> src_coord;  //Point sources.
   Vector<Real_t> src_value;
   Vector<size_t> src_scatter;

+ 110 - 0
include/fmm_node.txx

@@ -187,4 +187,114 @@ void FMM_Node<Node>::InitMultipole(PackedData data, bool own_data){
   }
 };
 
+template <class Node>
+template <class VTUData_t, class VTUNode_t>
+void FMM_Node<Node>::VTU_Data(VTUData_t& vtu_data, std::vector<VTUNode_t*>& nodes, int lod){
+  typedef typename VTUData_t::VTKReal_t VTKReal_t;
+  Node::VTU_Data(vtu_data, nodes, lod);
+
+  VTUData_t src_data;
+  VTUData_t surf_data;
+  VTUData_t trg_data;
+  { // Set src_data
+    VTUData_t& new_data=src_data;
+    new_data.value.resize(1);
+    new_data.name.push_back("src_value");
+    std::vector<VTKReal_t>& coord=new_data.coord;
+    std::vector<VTKReal_t>& value=new_data.value[0];
+
+    std::vector<int32_t>& connect=new_data.connect;
+    std::vector<int32_t>& offset =new_data.offset;
+    std::vector<uint8_t>& types  =new_data.types;
+
+    size_t point_cnt=0;
+    size_t connect_cnt=0;
+    for(size_t nid=0;nid<nodes.size();nid++){
+      VTUNode_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      for(size_t i=0;i<n->src_coord.Dim();i+=3){
+        coord.push_back(n->src_coord[i+0]);
+        coord.push_back(n->src_coord[i+1]);
+        coord.push_back(n->src_coord[i+2]);
+        connect_cnt++;
+        connect.push_back(point_cnt);
+        offset.push_back(connect_cnt);
+        types.push_back(1);
+        point_cnt++;
+      }
+      for(size_t i=0;i<n->src_value.Dim();i++){
+        value.push_back(n->src_value[i]);
+      }
+    }
+  }
+  { // Set surf_data
+    VTUData_t& new_data=surf_data;
+    new_data.value.resize(1);
+    new_data.name.push_back("surf_value");
+    std::vector<VTKReal_t>& coord=new_data.coord;
+    std::vector<VTKReal_t>& value=new_data.value[0];
+
+    std::vector<int32_t>& connect=new_data.connect;
+    std::vector<int32_t>& offset =new_data.offset;
+    std::vector<uint8_t>& types  =new_data.types;
+
+    size_t point_cnt=0;
+    size_t connect_cnt=0;
+    for(size_t nid=0;nid<nodes.size();nid++){
+      VTUNode_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      for(size_t i=0;i<n->surf_coord.Dim();i+=3){
+        coord.push_back(n->surf_coord[i+0]);
+        coord.push_back(n->surf_coord[i+1]);
+        coord.push_back(n->surf_coord[i+2]);
+        connect_cnt++;
+        connect.push_back(point_cnt);
+        offset.push_back(connect_cnt);
+        types.push_back(1);
+        point_cnt++;
+      }
+      for(size_t i=0;i<n->surf_value.Dim();i++){
+        value.push_back(n->surf_value[i]);
+      }
+    }
+  }
+  { // Set trg_data
+    VTUData_t& new_data=trg_data;
+    new_data.value.resize(1);
+    new_data.name.push_back("trg_value");
+    std::vector<VTKReal_t>& coord=new_data.coord;
+    std::vector<VTKReal_t>& value=new_data.value[0];
+
+    std::vector<int32_t>& connect=new_data.connect;
+    std::vector<int32_t>& offset =new_data.offset;
+    std::vector<uint8_t>& types  =new_data.types;
+
+    size_t point_cnt=0;
+    size_t connect_cnt=0;
+    for(size_t nid=0;nid<nodes.size();nid++){
+      VTUNode_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      for(size_t i=0;i<n->trg_coord.Dim();i+=3){
+        coord.push_back(n->trg_coord[i+0]);
+        coord.push_back(n->trg_coord[i+1]);
+        coord.push_back(n->trg_coord[i+2]);
+        connect_cnt++;
+        connect.push_back(point_cnt);
+        offset.push_back(connect_cnt);
+        types.push_back(1);
+        point_cnt++;
+      }
+      for(size_t i=0;i<n->trg_value.Dim();i++){
+        value.push_back(n->trg_value[i]);
+      }
+    }
+  }
+  AppendVTUData(vtu_data, src_data);
+  AppendVTUData(vtu_data, surf_data);
+  AppendVTUData(vtu_data, trg_data);
+}
+
 }//end namespace

+ 2 - 1
include/mpi_node.hpp

@@ -168,7 +168,8 @@ class MPI_Node: public TreeNode{
   /**
    * \brief Append node VTU data to vectors.
    */
-  virtual void VTU_Data(std::vector<Real_t>& coord, std::vector<Real_t>& value, std::vector<int32_t>& connect, std::vector<int32_t>& offset, std::vector<uint8_t>& types, int lod=-1);
+  template <class VTUData_t, class Node_t>
+  static void VTU_Data(VTUData_t& vtu_data, std::vector<Node_t*>& nodes, int lod);
 
   Vector<Real_t> pt_coord;   //coordinates of points
   Vector<Real_t> pt_value;   //value at points

+ 122 - 42
include/mpi_node.txx

@@ -417,50 +417,130 @@ void MPI_Node<T>::ReadVal(std::vector<Real_t> x,std::vector<Real_t> y, std::vect
   }
 }
 
-template <class T>
-void MPI_Node<T>::VTU_Data(std::vector<Real_t>& coord, std::vector<Real_t>& value, std::vector<int32_t>& connect, std::vector<int32_t>& offset, std::vector<uint8_t>& types, int lod){
-  if(!pt_coord.Dim()) return;
-  size_t point_cnt=coord.size()/COORD_DIM;
-  size_t connect_cnt=connect.size();
-  for(size_t i=0;i<pt_coord.Dim();i+=3){
-    coord.push_back(pt_coord[i+0]);
-    coord.push_back(pt_coord[i+1]);
-    coord.push_back(pt_coord[i+2]);
-    connect_cnt++;
-    connect.push_back(point_cnt);
-    offset.push_back(connect_cnt);
-    types.push_back(1);
-    point_cnt++;
-  }
-  for(size_t i=0;i<pt_value.Dim();i++) value.push_back(pt_value[i]);
-
-  Real_t* c=this->Coord();
-  Real_t s=pow(0.5,this->Depth());
-  size_t data_dof=pt_value.Dim()/(pt_coord.Dim()/dim);
-  for(int i=0;i<8;i++){
-    coord.push_back(c[0]+(i&1?1:0)*s);
-    coord.push_back(c[1]+(i&2?1:0)*s);
-    coord.push_back(c[2]+(i&4?1:0)*s);
-    for(int j=0;j<data_dof;j++) value.push_back(0.0);
-    connect.push_back(point_cnt+i);
+template <class VTUData_t>
+void AppendVTUData(VTUData_t& vtu_data, VTUData_t& new_data){
+  typedef typename VTUData_t::VTKReal_t VTKReal_t;
+
+  { // Append new_data to vtk_data
+    std::vector<VTKReal_t>&               new_coord=new_data.coord;
+    std::vector<std::string>&             new_name =new_data.name;
+    std::vector<std::vector<VTKReal_t> >& new_value=new_data.value;
+
+    std::vector<int32_t>& new_connect=new_data.connect;
+    std::vector<int32_t>& new_offset =new_data.offset;
+    std::vector<uint8_t>& new_types  =new_data.types;
+
+
+
+    std::vector<VTKReal_t>&               vtu_coord=vtu_data.coord;
+    std::vector<std::string>&             vtu_name =vtu_data.name;
+    std::vector<std::vector<VTKReal_t> >& vtu_value=vtu_data.value;
+
+    std::vector<int32_t>& vtu_connect=vtu_data.connect;
+    std::vector<int32_t>& vtu_offset =vtu_data.offset;
+    std::vector<uint8_t>& vtu_types  =vtu_data.types;
+
+
+
+    size_t old_pts=vtu_coord.size()/COORD_DIM;
+    size_t new_pts=new_coord.size()/COORD_DIM;
+
+    // New points
+    for(size_t i=0;i<new_coord.size();i++){
+      vtu_coord.push_back(new_coord[i]);
+    }
+
+    // Resize old DataArrays
+    for(size_t i=0;i<vtu_value.size();i++){
+      size_t curr_size=vtu_value[i].size();
+      size_t new_size=(curr_size*(old_pts+new_pts))/old_pts;
+      vtu_value[i].resize(new_size,0);
+    }
+
+    // Add new DataArrays
+    for(size_t i=0;i<new_name.size();i++){
+      vtu_name.push_back(new_name[i]);
+      vtu_value.push_back(std::vector<VTKReal_t>());
+      std::vector<VTKReal_t>& value=vtu_value.back();
+      if(new_pts) value.resize((new_value[i].size()*old_pts)/new_pts);
+      for(size_t j=0;j<new_value[i].size();j++){
+        value.push_back(new_value[i][j]);
+      }
+    }
+
+    size_t connect_update=old_pts;
+    size_t offset_update=vtu_connect.size();
+
+    for(size_t i=0;i<new_connect.size();i++){
+      vtu_connect.push_back(connect_update+new_connect[i]);
+    }
+    for(size_t i=0;i<new_offset.size();i++){
+      vtu_offset.push_back(offset_update+new_offset[i]);
+    }
+    for(size_t i=0;i<new_types.size();i++){
+      vtu_types.push_back(new_types[i]);
+    }
   }
-  offset.push_back(8+connect_cnt);
-  types.push_back(11);
-
-  {// Set point values.
-    Real_t* val_ptr=&value[point_cnt*data_dof];
-    std::vector<Real_t> x(2);
-    std::vector<Real_t> y(2);
-    std::vector<Real_t> z(2);
-    x[0]=c[0]; x[1]=c[0]+s;
-    y[0]=c[1]; y[1]=c[1]+s;
-    z[0]=c[2]; z[1]=c[2]+s;
-    this->ReadVal(x, y, z, val_ptr);
-    //Rearrrange data
-    //(x1,x2,x3,...,y1,y2,...z1,...) => (x1,y1,z1,x2,y2,z2,...)
-    Matrix<Real_t> M(data_dof,8,val_ptr,false);
-    M=M.Transpose();
+}
+
+template <class T>
+template <class VTUData_t, class Node_t>
+void MPI_Node<T>::VTU_Data(VTUData_t& vtu_data, std::vector<Node_t*>& nodes, int lod){
+  typedef typename VTUData_t::VTKReal_t VTKReal_t;
+
+  VTUData_t new_data;
+  { // Set new data
+    new_data.value.resize(1);
+    new_data.name.push_back("pt_value");
+    std::vector<VTKReal_t>& coord=new_data.coord;
+    std::vector<VTKReal_t>& value=new_data.value[0];
+
+    std::vector<int32_t>& connect=new_data.connect;
+    std::vector<int32_t>& offset =new_data.offset;
+    std::vector<uint8_t>& types  =new_data.types;
+
+    size_t point_cnt=0;
+    size_t connect_cnt=0;
+    for(size_t nid=0;nid<nodes.size();nid++){
+      Node_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      for(size_t i=0;i<n->pt_coord.Dim();i+=3){
+        coord.push_back(n->pt_coord[i+0]);
+        coord.push_back(n->pt_coord[i+1]);
+        coord.push_back(n->pt_coord[i+2]);
+        connect_cnt++;
+        connect.push_back(point_cnt);
+        offset.push_back(connect_cnt);
+        types.push_back(1);
+        point_cnt++;
+      }
+      for(size_t i=0;i<n->pt_value.Dim();i++){
+        value.push_back(n->pt_value[i]);
+      }
+    }
+    size_t value_dof=value.size()/point_cnt;
+    assert(value_dof*point_cnt==value.size());
+    for(size_t nid=0;nid<nodes.size();nid++){
+      Node_t* n=nodes[nid];
+      if(n->IsGhost() || !n->IsLeaf()) continue;
+
+      Real_t* c=n->Coord();
+      Real_t s=pow(0.5,n->Depth());
+      for(int i=0;i<8;i++){
+        coord.push_back(c[0]+(i&1?1:0)*s);
+        coord.push_back(c[1]+(i&2?1:0)*s);
+        coord.push_back(c[2]+(i&4?1:0)*s);
+        for(int j=0;j<value_dof;j++) value.push_back(0.0);
+        connect.push_back(point_cnt+i);
+        connect_cnt++;
+      }
+      offset.push_back(connect_cnt);
+      types.push_back(11);
+      point_cnt+=8;
+    }
   }
+  AppendVTUData(vtu_data, new_data);
 }
 
 }//end namespace

+ 16 - 0
include/mpi_tree.hpp

@@ -8,6 +8,7 @@
 
 #include <mpi.h>
 #include <vector>
+#include <string>
 
 #include <pvfmm_common.hpp>
 #include <mortonid.hpp>
@@ -105,6 +106,21 @@ class MPI_Tree: public Tree<TreeNode>{
    */
   void ConstructLET(BoundaryType bndry=FreeSpace);
 
+  template <class VTKReal>
+  struct VTUData_t{
+    typedef VTKReal VTKReal_t;
+
+    // Point data
+    std::vector<VTKReal_t> coord;
+    std::vector<std::vector<VTKReal_t> > value;
+    std::vector<std::string> name;
+
+    // Cell data
+    std::vector<int32_t> connect;
+    std::vector<int32_t> offset;
+    std::vector<uint8_t> types;
+  };
+
   /**
    * \brief Write to a <fname>.vtu file.
    */

+ 36 - 34
include/mpi_tree.txx

@@ -2098,37 +2098,29 @@ inline bool isLittleEndian(){
 
 template <class TreeNode>
 void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
-  typedef double VTKData_t;
+  typedef double VTKReal_t;
+
   int myrank, np;
   MPI_Comm_size(*Comm(),&np);
   MPI_Comm_rank(*Comm(),&myrank);
 
-  std::vector<Real_t> coord_;  //Coordinates of octant corners.
-  std::vector<Real_t> value_;  //Data value at points.
-  std::vector<VTKData_t> coord;  //Coordinates of octant corners.
-  std::vector<VTKData_t> value;  //Data value at points.
-  std::vector<int32_t> mpi_rank;  //MPI_Rank at points.
-  std::vector<int32_t> connect;  //Cell connectivity.
-  std::vector<int32_t> offset ;  //Cell offset.
-  std::vector<uint8_t> types  ;  //Cell types.
+  VTUData_t<VTKReal_t> vtu_data;
+  TreeNode::VTU_Data(vtu_data, this->GetNodeList(), lod);
+
+  std::vector<VTKReal_t>&               coord=vtu_data.coord;
+  std::vector<std::string>&             name =vtu_data.name;
+  std::vector<std::vector<VTKReal_t> >& value=vtu_data.value;
+
+  std::vector<int32_t>& connect=vtu_data.connect;
+  std::vector<int32_t>& offset =vtu_data.offset;
+  std::vector<uint8_t>& types  =vtu_data.types;
 
-  //Build list of octant corner points.
-  Node_t* n=this->PreorderFirst();
-  while(n!=NULL){
-    if(!n->IsGhost() && n->IsLeaf())
-      n->VTU_Data(coord_, value_, connect, offset, types, lod);
-    n=this->PreorderNxt(n);
-  }
-  for(size_t i=0;i<coord_.size();i++) coord.push_back(coord_[i]);
-  for(size_t i=0;i<value_.size();i++) value.push_back(value_[i]);
   int pt_cnt=coord.size()/COORD_DIM;
-  int dof=(pt_cnt?value.size()/pt_cnt:0);
-  assert(value.size()==(size_t)pt_cnt*dof);
   int cell_cnt=types.size();
 
-  mpi_rank.resize(pt_cnt);
+  std::vector<int32_t> mpi_rank;  //MPI_Rank at points.
   int new_myrank=myrank;//rand();
-  for(int i=0;i<pt_cnt;i++) mpi_rank[i]=new_myrank;
+  mpi_rank.resize(pt_cnt,new_myrank);
 
   //Open file for writing.
   std::stringstream vtufname;
@@ -2148,15 +2140,19 @@ void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
 
   //---------------------------------------------------------------------------
   vtufile<<"      <Points>\n";
-  vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKData_t)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
-  data_size+=sizeof(uint32_t)+coord.size()*sizeof(VTKData_t);
+  vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKReal_t)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
+  data_size+=sizeof(uint32_t)+coord.size()*sizeof(VTKReal_t);
   vtufile<<"      </Points>\n";
   //---------------------------------------------------------------------------
   vtufile<<"      <PointData>\n";
-  vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKData_t)*8<<"\" NumberOfComponents=\""<<dof<<"\" Name=\"value\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
-  data_size+=sizeof(uint32_t)+value.size()*sizeof(VTKData_t);
-  vtufile<<"        <DataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
-  data_size+=sizeof(uint32_t)+mpi_rank.size()*sizeof(int32_t);
+  for(size_t i=0;i<name.size();i++) if(value[i].size()){ // value
+    vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKReal_t)*8<<"\" NumberOfComponents=\""<<value[i].size()/pt_cnt<<"\" Name=\""<<name[i]<<"\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
+    data_size+=sizeof(uint32_t)+value[i].size()*sizeof(VTKReal_t);
+  }
+  { // mpi_rank
+    vtufile<<"        <DataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
+    data_size+=sizeof(uint32_t)+mpi_rank.size()*sizeof(int32_t);
+  }
   vtufile<<"      </PointData>\n";
   //---------------------------------------------------------------------------
   //---------------------------------------------------------------------------
@@ -2170,7 +2166,7 @@ void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
   vtufile<<"      </Cells>\n";
   //---------------------------------------------------------------------------
   //vtufile<<"      <CellData>\n";
-  //vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKData_t)*8<<"\" Name=\"Velocity\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
+  //vtufile<<"        <DataArray type=\"Float"<<sizeof(VTKReal_t)*8<<"\" Name=\"Velocity\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
   //vtufile<<"      </CellData>\n";
   //---------------------------------------------------------------------------
 
@@ -2181,8 +2177,10 @@ void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
   vtufile<<"    _";
 
   int32_t block_size;
-  block_size=coord   .size()*sizeof(VTKData_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&coord   [0], coord   .size()*sizeof(VTKData_t));
-  block_size=value   .size()*sizeof(VTKData_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&value   [0], value   .size()*sizeof(VTKData_t));
+  block_size=coord   .size()*sizeof(VTKReal_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&coord   [0], coord   .size()*sizeof(VTKReal_t));
+  for(size_t i=0;i<name.size();i++) if(value[i].size()){ // value
+    block_size=value[i].size()*sizeof(VTKReal_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&value[i][0], value[i].size()*sizeof(VTKReal_t));
+  }
   block_size=mpi_rank.size()*sizeof(  int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&mpi_rank[0], mpi_rank.size()*sizeof(  int32_t));
 
   block_size=connect.size()*sizeof(int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&connect[0], connect.size()*sizeof(int32_t));
@@ -2206,11 +2204,15 @@ void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
   pvtufile<<"<VTKFile type=\"PUnstructuredGrid\">\n";
   pvtufile<<"  <PUnstructuredGrid GhostLevel=\"0\">\n";
   pvtufile<<"      <PPoints>\n";
-  pvtufile<<"        <PDataArray type=\"Float"<<sizeof(VTKData_t)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\"/>\n";
+  pvtufile<<"        <PDataArray type=\"Float"<<sizeof(VTKReal_t)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\"/>\n";
   pvtufile<<"      </PPoints>\n";
   pvtufile<<"      <PPointData>\n";
-  pvtufile<<"        <PDataArray type=\"Float"<<sizeof(VTKData_t)*8<<"\" NumberOfComponents=\""<<dof<<"\" Name=\"value\"/>\n";
-  pvtufile<<"        <PDataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\"/>\n";
+  for(size_t i=0;i<name.size();i++) if(value[i].size()){ // value
+    pvtufile<<"        <PDataArray type=\"Float"<<sizeof(VTKReal_t)*8<<"\" NumberOfComponents=\""<<value[i].size()/pt_cnt<<"\" Name=\""<<name[i]<<"\"/>\n";
+  }
+  { // mpi_rank
+    pvtufile<<"        <PDataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\"/>\n";
+  }
   pvtufile<<"      </PPointData>\n";
   {
     // Extract filename from path.

+ 1 - 1
include/pvfmm.hpp

@@ -29,7 +29,7 @@ typedef FMM_Node<Cheb_Node<double> > ChebFMM_Node;
 typedef FMM_Cheb<ChebFMM_Node>       ChebFMM;
 typedef FMM_Tree<ChebFMM>            ChebFMM_Tree;
 typedef ChebFMM_Node::NodeData       ChebFMM_Data;
-typedef void (*ChebFn)(const double* , int , double*);
+typedef ChebFMM_Node::Function_t<double> ChebFn;
 
 ChebFMM_Tree* ChebFMM_CreateTree(int cheb_deg, int data_dim, ChebFn fn_ptr, std::vector<double>& trg_coord, MPI_Comm& comm,
                                  double tol=1e-6, int max_pts=100, BoundaryType bndry=FreeSpace, int init_depth=0){