output.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_output_h
#define dealii_mesh_worker_output_h

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/base/mg_level_object.h>


DEAL_II_NAMESPACE_OPEN

namespace MeshWorker
{
  namespace Assembler
  {

    class GnuplotPatch
    {
    public:
      GnuplotPatch();

      void initialize (const unsigned int n_points,
                       const unsigned int n_vectors);

      void initialize_stream (std::ostream &stream);

      template <int dim>
      void initialize_info(DoFInfo<dim> &info, bool face);

      template <int dim>
      void assemble(const DoFInfo<dim> &info);

      template <int dim>
      void assemble(const DoFInfo<dim> &info1,
                    const DoFInfo<dim> &info2);

    private:
      template <typename T>
      void write(const T &t) const;

      void write_endl () const;

      unsigned int n_vectors;
      unsigned int n_points;

      std::ostream *os;
    };

//----------------------------------------------------------------------//

    template <typename T>
    inline void
    GnuplotPatch::write(const T &d) const
    {
      if (os == nullptr)
        deallog << d;
      else
        (*os) << d;
    }


    inline void
    GnuplotPatch::write_endl() const
    {
      if (os == nullptr)
        deallog << std::endl;
      else
        (*os) << std::endl;
    }


    inline
    GnuplotPatch::GnuplotPatch()
      :
      n_vectors(numbers::invalid_unsigned_int),
      n_points(numbers::invalid_unsigned_int),
      os(nullptr)
    {}


    inline void
    GnuplotPatch::initialize (const unsigned int np,
                              const unsigned int nv)
    {
      n_vectors = nv;
      n_points = np;
    }


    inline void
    GnuplotPatch::initialize_stream (std::ostream &stream)
    {
      os = &stream;
    }


    template <int dim>
    inline void
    GnuplotPatch::initialize_info(DoFInfo<dim> &info, bool face)
    {
      if (face)
        info.initialize_quadrature(Utilities::fixed_power<dim-1>(n_points), n_vectors+dim);
      else
        info.initialize_quadrature(Utilities::fixed_power<dim>(n_points), n_vectors+dim);
    }


    template <int dim>
    inline void
    GnuplotPatch::assemble(const DoFInfo<dim> &info)
    {
      const unsigned int np = info.n_quadrature_points();
      const unsigned int nv = info.n_quadrature_values();
      const unsigned int patch_dim = (info.face_number == numbers::invalid_unsigned_int)
                                     ? dim : (dim-1);
      const unsigned int row_length = n_points;
      // If patches are 1D, end the
      // patch after a row, else end
      // it after a square
      const unsigned int row_length2 = (patch_dim==1) ? row_length : (row_length*row_length);

//      AssertDimension(np, Utilities::fixed_power<dim>(n_points));
      AssertDimension(nv, n_vectors+dim);


      for (unsigned int k=0; k<np; ++k)
        {
          if (k % row_length == 0)
            write_endl();
          if (k % row_length2 == 0)
            write_endl();

          for (unsigned int i=0; i<nv; ++i)
            {
              write(info.quadrature_value(k,i));
              write('\t');
            }
          write_endl();
        }
    }


    template <int dim>
    inline void
    GnuplotPatch::assemble(const DoFInfo<dim> &info1, const DoFInfo<dim> &info2)
    {
      assemble(info1);
      assemble(info2);
    }
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif