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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_output_h
#define dealii_mesh_worker_output_h

#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/utilities.h>

#include <deal.II/lac/block_vector.h>

#include <deal.II/meshworker/dof_info.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);
    }
  } // namespace Assembler
} // namespace MeshWorker

DEAL_II_NAMESPACE_CLOSE

#endif