graphics/vtkGraphics.hh

Go to the documentation of this file.

/*
 * VtkGraphics.hh
 *
 *  Created on: 21.09.2011
 *      Author: Robert Gruhlke und Christian Heier
 */

#ifndef VtkGraphics_HH_
#define VtkGraphics_HH_
#include "basics.hh"
#include "matrixBasis.hh"
#include "basis.hh"

namespace graphics {

  class VtkGraphics: public OutputBase {
  public:
    VtkGraphics();

    template<class F, class G>
    VtkGraphics(const concepts::Space<G>& spc, const std::string filename,
        const concepts::Vector<F>& sol,
        const concepts::ElementFunction<F, G>* fun = 0, const uint dim = 2);

    template<class F, class G>
    VtkGraphics(const concepts::Space<G>& spc, const std::string filename,
        const concepts::ElementFormula<F, G>& frm, const uint dim = 2);

    template<class G>
    VtkGraphics(const concepts::Space<G>& spc, const std::string filename,
        const uint dim = 2);

    VtkGraphics(concepts::Mesh& msh, const std::string filename,
        const uint dim = 2, const uint points = 5);

    virtual ~VtkGraphics();

    template<class G>
    inline void clearMapping(
        std::map<std::string, concepts::DenseMatrix<G>*>& dense_ptr_) {

      delete (dense_ptr_["x"]);

      if (dim_ > 1)
        delete (dense_ptr_["y"]);
      if (dim_ > 2)
        delete (dense_ptr_["z"]);

      if (dim_ > 1) {
        delete (dense_ptr_["edgemsh"]);
        delete (dense_ptr_["vtxmsh"]);
      }
      delete (dense_ptr_["Attr"]);
      delete (dense_ptr_["edgAttr"]);

      //stores the gradient
      //if (gradAndValue_)
      //  delete (dense_ptr_["grad"]);

      //store informations to plot solution
      if (data_) {
        delete (dense_ptr_["p"]);
        delete (dense_ptr_["u"]);
        // Cell attribute on each quadrature points
        delete (dense_ptr_["pattr"]);
        delete (dense_ptr_["w"]);
        delete (dense_ptr_["msh"]);
      }
      dense_ptr_.clear();
    }
    ;

  private:

    // Create the matrices
    template<class G>
    void createMatrices_(bool data, uint count[], std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_);

    //Store data as .mat file
    template<class G, class F>
    void storeData_(bool data, std::string filename, std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_,
        concepts::DenseMatrix<F>* functionValue = 0);

    void writeHead_(std::string&, std::ostringstream& stream);
    template<class G, class F>
    void writePoints_(std::ostringstream& stream, std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_,
        concepts::DenseMatrix<F>* functionValue = 0);
    template<class G>
    void writeMesh_(std::ostringstream& stream, std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_);
    template<class G>
    void writeEdgeMesh_(std::ostringstream& stream, std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_);
    template<class F>
    void writeValues_(std::ostringstream& stream,
        concepts::DenseMatrix<F>* functionValue = 0);

    template<class G>
    void writeWeightsAndAttr_(std::ostringstream& stream, std::map<std::string,
        concepts::DenseMatrix<G>*>& dense_ptr_);

    std::string vtkEnding_(const std::string& filename) {
      uint size = filename.size();
      if (size < 4 || std::string(filename, size - 4, 4) != ".vtk")
        return (filename + ".vtk");
      return filename;
    }

    // for destructor
    uint data_;

    // Spatial dimension of the space
    uint dim_;

    bool isComplex_;

    // dimension of the output: 1 if output is solution | dim_ if output is gradient
    uint outputDimension_;

  };

  template<class F, class G>
  VtkGraphics::VtkGraphics(const concepts::Space<G>& spc,
      const std::string filename, const concepts::Vector<F>& sol,
      const concepts::ElementFunction<F, G>* fun, const uint dim) :
    OutputBase(vtkEnding_(filename)) {

    dim_ = dim;

    // Post processor (global)
    concepts::GlobalPostprocess<G> postProcess(spc);

    //output dimension for ex. n = 1 if fun = Value or n = 2 if fun = Grad
    if (fun)
      outputDimension_ = fun->n();
    else
      outputDimension_ = 1;

    //is functionValue Complex or Real?
    if (typeid(F) == typeid(concepts::Cmplx))
      isComplex_ = true;
    if (typeid(F) == typeid(concepts::Real))
      isComplex_ = false;

    //declare mapping
    std::map<std::string, concepts::DenseMatrix<G>*> dense_ptr_;

    //declare FunctionValues u
    concepts::DenseMatrix<F>* functionValue;

    // Post processor on cell
    MatrixBaseDataCell<F, G> bdc(dense_ptr_, functionValue, sol, fun);

    //calculate size of the matrices / size of data
    postProcess(bdc);

    //Now the we know the size of the matrices, so create them.
    createMatrices_(true, bdc.dimensionsOfMatrices, dense_ptr_);

    //Now the we know the number of evaluations of the function, so create the valuevector.
    functionValue = new concepts::DenseMatrix<F>(outputDimension_,
        bdc.dimensionsOfMatrices[0]);

    //the MatrixBaseDataCell should know that everything is counted i.e. it can read data in the matrices
    bdc.counted = true;

    //run GlobalPostprocess with bdc to fill the concepts::DenseMatrix<> with mesh/solution informations.
    postProcess(bdc);

    //store the matrices
    storeData_(true, filename, dense_ptr_, functionValue);

    //deletes dense matrix pointer
    delete (functionValue);

    //Deletes the Dense matrices in the mapping
    clearMapping(dense_ptr_);
  }

  template<class F, class G>
  VtkGraphics::VtkGraphics(const concepts::Space<G>& spc,
      const std::string filename, const concepts::ElementFormula<F, G>& frm,
      const uint dim) :
    OutputBase(vtkEnding_(filename)) {

    //save dimension
    dim_ = dim;

    if (typeid(F) == typeid(concepts::Cmplx))
      isComplex_ = true;
    else
      isComplex_ = false;

    //output dimension for ex. n = 1 if fun = Value or n = 2 if fun = Grad
    outputDimension_ = Size<F>::n();

    // Post processor (global)
    concepts::GlobalPostprocess<G> postProcess(spc);

    //declare mapping
    std::map<std::string, concepts::DenseMatrix<G>*> dense_ptr_;

    //declare FunctionValues u
    concepts::DenseMatrix<F>* functionValue;

    // Post processor on cell
    MatrixBaseElementFormulaCell<F, G>
        bdc(dense_ptr_, functionValue, dim_, frm);

    //calculate size of the matrices / size of data
    postProcess(bdc);

    //Now we know the size of the matrices, so lets inatialize them.
    createMatrices_(true, bdc.dimensionsOfMatrices, dense_ptr_);

    //Now the we know the number of evaluations of the function, so create the valuevector.
    functionValue = new concepts::DenseMatrix<F>(outputDimension_,
        bdc.dimensionsOfMatrices[0]);

    //the MatrixBaseMeshCell should know that everything is counted i.e. it can read data in the matrices
    bdc.counted = true;

    //run GlobalPostprocess with bdc to fill the concepts::DenseMatrix<> with mesh/solution informations.
    postProcess(bdc);

    //store the matrices
    storeData_(true, filename, dense_ptr_, functionValue);

    //Delete the Dense matrices in the mapping and function value
    delete (functionValue);
    clearMapping(dense_ptr_);

  }

  template<class G>
  VtkGraphics::VtkGraphics(const concepts::Space<G>& spc,
      const std::string filename, const uint dim) :
    OutputBase(vtkEnding_(filename)) {

    dim_ = dim;
    isComplex_ = false;
    outputDimension_ = 0;

    // Post processor (global)
    concepts::GlobalPostprocess<G> postProcess(spc);

    //declare mapping
    std::map<std::string, concepts::DenseMatrix<G>*> dense_ptr_;

    // Post processor on cell
    MatrixBaseMeshCell bdc(dense_ptr_, dim_);

    //calculate size of the matrices / size of data
    postProcess(bdc);

    //Now the we know the size of the matrices, so create them.
    createMatrices_(false, bdc.dimensionsOfMatrices, dense_ptr_);

    //the MatrixBaseDataCell should know that everything is counted i.e. it can read data in the matrices
    bdc.counted = true;

    //run GlobalPostprocess with bdc to fill the concepts::DenseMatrix<> with mesh/solution informations.
    postProcess(bdc);

    //store the matrices
    storeData_<concepts::Real, concepts::Real> (false, filename, dense_ptr_);

    //Deletes the Dense matrices in the mapping
    clearMapping(dense_ptr_);
  }

  template<class G>
  void VtkGraphics::createMatrices_(bool data, uint count[], std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_) {

    data_ = data;
    // coordinates

    dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("x",
        new concepts::DenseMatrix<G>(1, count[0])));
    if (dim_ > 1)
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("y",
          new concepts::DenseMatrix<G>(1, count[0])));
    if (dim_ > 2)
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("z",
          new concepts::DenseMatrix<G>(1, count[0])));

    //TODO add FaceMesh

    dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>(
        "edgAttr", new concepts::DenseMatrix<G>(1, count[2])));

    if (dim_ > 1) {
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>(
          "edgemsh", new concepts::DenseMatrix<G>(2, count[2])));
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>(
          "vtxmsh", new concepts::DenseMatrix<G>(5, count[3])));
    }
    dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("Attr",
        new concepts::DenseMatrix<G>(1, count[3])));

    //stores the gradient
    //if (gradAndValue_)
    //  dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("grad",
    //      new concepts::DenseMatrix<G> (dim_, count[0])));

    //store informations to plot solution
    if (data) {
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("p",
          new concepts::DenseMatrix<G>(2, count[3])));
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("u",
          new concepts::DenseMatrix<G>(outputDimension_, count[0])));
      // Cell attribute on each quadrature points
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>(
          "pattr", new concepts::DenseMatrix<G>(1, count[0])));
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>("w",
          new concepts::DenseMatrix<G>(1, count[0])));
      dense_ptr_.insert(std::pair<std::string, concepts::DenseMatrix<G>*>(
          "msh", new concepts::DenseMatrix<G>(5, count[1])));
    }

  }//end Create matrices


  //Save the matrices that hold informations
  template<class G, class F>
  void VtkGraphics::storeData_(bool data, std::string filename, std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_,
      concepts::DenseMatrix<F>* functionValue) {

    std::ostringstream stream;

    //write head of vtk file
    writeHead_(filename, stream);

    //write the figure-points in the vtk file
    writePoints_(stream, dense_ptr_, functionValue);

    //writes the mesh in the vtk file
    //(i.e. the quads = [point1, point2, point3, point4] that builds the mesh)
    if (data)
      writeMesh_(stream, dense_ptr_);
    //writes data for the edge
    writeEdgeMesh_(stream, dense_ptr_);

    //writes data in the vtk file. The value/gradient/laplace of the solution,
    //the integration weights of the integration rule and attributes
    if (data) {
      writeValues_(stream, functionValue);
      writeWeightsAndAttr_(stream, dense_ptr_);
    }

    //store the written data in ofs
    *ofs_ << stream.str();
  }

  template<class G>
  void VtkGraphics::writeMesh_(std::ostringstream& stream, std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_) {
    stream << "POLYGONS" << " " << dense_ptr_["msh"]->dimY() << " "
        << dense_ptr_["msh"]->dimY() * 5 << "\n";

    for (uint i = 0; i < dense_ptr_["msh"]->dimY(); ++i)
      stream << "4 " << dense_ptr_["msh"]->operator()(0, i) - 1 << " "
          << dense_ptr_["msh"]->operator()(1, i) - 1 << " "
          << dense_ptr_["msh"]->operator()(2, i) - 1 << " "
          << dense_ptr_["msh"]->operator()(3, i) - 1 << '\n';
  }

  template<class G>
  void VtkGraphics::writeEdgeMesh_(std::ostringstream& stream, std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_) {
    stream << "LINES" << " " << dense_ptr_["edgemsh"]->dimY() << " "
        << dense_ptr_["edgemsh"]->dimY() * 3 << "\n";

    for (uint i = 0; i < dense_ptr_["edgemsh"]->dimY(); ++i)
      stream << "2 " << dense_ptr_["edgemsh"]->operator()(0, i) - 1 << " "
          << dense_ptr_["edgemsh"]->operator()(1, i) - 1 << "\n";
  }

  template<class F>
  void VtkGraphics::writeValues_(std::ostringstream& stream,
      concepts::DenseMatrix<F>* functionValue) {
    stream << "POINT_DATA " << functionValue->dimY() << "\n";
    if (isComplex_) {
      stream << "VECTORS solution FLOAT\n";
    } else {
      stream << "SCALARS solution FLOAT\n";
      stream << "LOOKUP_TABLE default\n";
    }
    //stream << "LOOKUP_TABLE default\n";

    if (isComplex_) {
      concepts::DenseMatrix<concepts::Cmplx>* functionValueCmplx =
          dynamic_cast<concepts::DenseMatrix<concepts::Cmplx>*> (functionValue);

      for (uint i = 0; i < functionValue->dimY(); ++i)
        stream << (functionValueCmplx->operator()(0, i)).real() << " "
            << (functionValueCmplx->operator()(0, i)).imag() << " 0\n";
    } else {
      for (uint i = 0; i < functionValue->dimY(); ++i)
        stream << functionValue->operator()(0, i) << "\n";
    }

  }

  template<class G>
  void VtkGraphics::writeWeightsAndAttr_(std::ostringstream& stream, std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_) {
    stream << "SCALARS weights FLOAT\n";
    stream << "LOOKUP_TABLE default\n";
    for (uint i = 0; i < dense_ptr_["w"]->dimY(); ++i)
      stream << dense_ptr_["w"]->operator()(0, i) << "\n";

    stream << "SCALARS pointAttributes FLOAT\n";
    stream << "LOOKUP_TABLE default\n";
    for (uint i = 0; i < dense_ptr_["pattr"]->dimY(); ++i)
      stream << dense_ptr_["pattr"]->operator()(0, i) << "\n";

  }

  template<class G, class F>
  void VtkGraphics::writePoints_(std::ostringstream& stream, std::map<
      std::string, concepts::DenseMatrix<G>*>& dense_ptr_,
      concepts::DenseMatrix<F>* functionValue) {

    stream << "POINTS " << dense_ptr_["x"]->dimY() << " float\n";
    if (isComplex_ && functionValue != 0) {
      //std::cout << "isComplex" << std::endl;
      //concepts::DenseMatrix<concepts::Cmplx> * functionValueCmplx =
      //    dynamic_cast<concepts::DenseMatrix<concepts::Cmplx>*> (functionValue);
      for (uint i = 0; i < dense_ptr_["x"]->dimY(); ++i) {
        stream << dense_ptr_["x"]->operator()(0, i) << " "
            << dense_ptr_["y"]->operator()(0, i);
        // if (outputDimension_ == 1) {
        //   stream << " " << (functionValueCmplx->operator()(0, i)).real()
        //       << "\n";
        // } else
        stream << " " << 0 << "\n";
      }
    } else {
      for (uint i = 0; i < dense_ptr_["x"]->dimY(); ++i) {
        stream << dense_ptr_["x"]->operator()(0, i) << " "
            << dense_ptr_["y"]->operator()(0, i);
        if (outputDimension_ == 1 && functionValue != 0)
          stream << " " << functionValue->operator()(0, i) << "\n";
        else
          stream << " " << 0 << "\n";
      }

    }

  }

} // namespace graphics


#endif /* VtkGraphics_HH_ */