data_out_dof_data.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_data_out_dof_data_h
#define dealii_data_out_dof_data_h



#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/data_out_base.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/tria.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/numerics/data_postprocessor.h>
#include <deal.II/numerics/data_component_interpretation.h>

#include <memory>

DEAL_II_NAMESPACE_OPEN

namespace Exceptions
{
  namespace DataOutImplementation
  {
    DeclException1 (ExcInvalidNumberOfSubdivisions,
                    int,
                    << "The number of subdivisions per patch, " << arg1
                    << ", is not valid. It needs to be greater or equal to "
                    "one, or zero if you want it to be determined "
                    "automatically.");

    DeclExceptionMsg (ExcNoTriangulationSelected,
                      "For the operation you are attempting, you first need to "
                      "tell the DataOut or related object which DoFHandler or "
                      "triangulation you would like to work on.");

    DeclExceptionMsg (ExcNoDoFHandlerSelected,
                      "For the operation you are attempting, you first need to "
                      "tell the DataOut or related object which DoFHandler "
                      "you would like to work on.");

    DeclException3 (ExcInvalidVectorSize,
                    int, int, int,
                    << "The vector has size " << arg1
                    << " but the DoFHandler object says that there are " << arg2
                    << " degrees of freedom and there are " << arg3
                    << " active cells. The size of your vector needs to be"
                    << " either equal to the number of degrees of freedom, or"
                    << " equal to the number of active cells.");
    DeclException2 (ExcInvalidCharacter,
                    std::string, size_t,
                    << "Please use only the characters [a-zA-Z0-9_<>()] for" << std::endl
                    << "description strings since some graphics formats will only accept these."
                    << std::endl
                    << "The string you gave was <" << arg1
                    << ">, within which the invalid character is <" << arg1[arg2]
                    << ">." << std::endl);
    DeclExceptionMsg (ExcOldDataStillPresent,
                      "When attaching a triangulation or DoFHandler object, it is "
                      "not allowed if old data vectors are still referenced. If "
                      "you want to reuse an object of the current type, you first "
                      "need to call the 'clear_data_vector()' function.");
    DeclException2 (ExcInvalidNumberOfNames,
                    int, int,
                    << "You have to give one name per component in your "
                    << "data vector. The number you gave was " << arg1
                    << ", but the number of components is " << arg2
                    << ".");
    DeclExceptionMsg (ExcIncompatibleDatasetNames,
                      "While merging sets of patches, the two sets to be merged "
                      "need to refer to data that agrees on the names of the "
                      "various variables represented. In other words, you "
                      "cannot merge sets of patches that originate from "
                      "entirely unrelated simulations.");
    DeclExceptionMsg (ExcIncompatiblePatchLists,
                      "While merging sets of patches, the two sets to be merged "
                      "need to refer to data that agrees on the number of "
                      "subdivisions and other properties. In other words, you "
                      "cannot merge sets of patches that originate from "
                      "entirely unrelated simulations.");

    DeclException2 (ExcInvalidVectorDeclaration,
                    int, std::string,
                    << "When declaring that a number of components in a data "
                    << "set to be output logically form a vector instead of "
                    << "simply a set of scalar fields, you need to specify "
                    << "this for all relevant components. Furthermore, "
                    << "vectors must always consist of exactly <dim> "
                    << "components. However, the vector component at "
                    << "position " << arg1 << " with name <" << arg2
                    << "> does not satisfy these conditions.");
  }
}


namespace internal
{
  namespace DataOutImplementation
  {
    enum class ComponentExtractor
    {
      real_part,
      imaginary_part
    };


    template <typename DoFHandlerType>
    class DataEntryBase
    {
    public:
      DataEntryBase (const DoFHandlerType           *dofs,
                     const std::vector<std::string> &names,
                     const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation);

      DataEntryBase (const DoFHandlerType *dofs,
                     const DataPostprocessor<DoFHandlerType::space_dimension> *data_postprocessor);

      virtual ~DataEntryBase () = default;

      virtual
      double
      get_cell_data_value (const unsigned int cell_number,
                           const ComponentExtractor extract_component) const = 0;

      virtual
      void
      get_function_values
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<double> &patch_values) const = 0;

      virtual
      void
      get_function_values
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<::Vector<double> > &patch_values_system) const = 0;

      virtual
      void
      get_function_gradients
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<Tensor<1,DoFHandlerType::space_dimension> > &patch_gradients) const = 0;

      virtual
      void
      get_function_gradients
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<std::vector<Tensor<1,DoFHandlerType::space_dimension> > > &patch_gradients_system) const = 0;

      virtual
      void
      get_function_hessians
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<Tensor<2,DoFHandlerType::space_dimension> > &patch_hessians) const = 0;

      virtual
      void
      get_function_hessians
      (const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &fe_patch_values,
       const ComponentExtractor extract_component,
       std::vector<std::vector< Tensor<2,DoFHandlerType::space_dimension> > > &patch_hessians_system) const = 0;

      virtual bool is_complex_valued () const = 0;

      virtual void clear () = 0;

      virtual std::size_t memory_consumption () const = 0;

      SmartPointer<const DoFHandlerType> dof_handler;

      const std::vector<std::string> names;

      const std::vector<DataComponentInterpretation::DataComponentInterpretation>
      data_component_interpretation;

      SmartPointer<const ::DataPostprocessor<DoFHandlerType::space_dimension> > postprocessor;

      unsigned int n_output_variables;
    };


    template <int dim, int spacedim>
    struct ParallelDataBase
    {
      ParallelDataBase (const unsigned int n_datasets,
                        const unsigned int n_subdivisions,
                        const std::vector<unsigned int> &n_postprocessor_outputs,
                        const Mapping<dim,spacedim> &mapping,
                        const std::vector<std::shared_ptr<::hp::FECollection<dim,spacedim> > > &finite_elements,
                        const UpdateFlags update_flags,
                        const bool        use_face_values);

      ParallelDataBase (const ParallelDataBase &data);

      template <typename DoFHandlerType>
      void reinit_all_fe_values(std::vector<std::shared_ptr<DataEntryBase<DoFHandlerType> > > &dof_data,
                                const typename ::Triangulation<dim,spacedim>::cell_iterator &cell,
                                const unsigned int face = numbers::invalid_unsigned_int);

      const FEValuesBase<dim,spacedim> &
      get_present_fe_values(const unsigned int dataset) const;

      void resize_system_vectors(const unsigned int n_components);

      const unsigned int n_datasets;
      const unsigned int n_subdivisions;

      DataPostprocessorInputs::Scalar<spacedim>          patch_values_scalar;
      DataPostprocessorInputs::Vector<spacedim>          patch_values_system;
      std::vector<std::vector<::Vector<double> > > postprocessed_values;

      const ::hp::MappingCollection<dim,spacedim> mapping_collection;
      const std::vector<std::shared_ptr<::hp::FECollection<dim,spacedim> > > finite_elements;
      const UpdateFlags update_flags;

      std::vector<std::shared_ptr<::hp::FEValues<dim,spacedim> > > x_fe_values;
      std::vector<std::shared_ptr<::hp::FEFaceValues<dim,spacedim> > > x_fe_face_values;
    };
  }
}


//TODO: Most of the documentation of DataOut_DoFData applies to DataOut.

template <typename DoFHandlerType, int patch_dim, int patch_space_dim=patch_dim>
class DataOut_DoFData : public DataOutInterface<patch_dim,patch_space_dim>
{
public:

  typedef typename Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension>::cell_iterator cell_iterator;
  typedef typename Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension>::active_cell_iterator active_cell_iterator;

public:

  enum DataVectorType
  {
    type_dof_data,

    type_cell_data,

    type_automatic
  };

  DataOut_DoFData ();

  virtual ~DataOut_DoFData ();

  void attach_dof_handler (const DoFHandlerType &);

  void attach_triangulation (const Triangulation<DoFHandlerType::dimension,
                             DoFHandlerType::space_dimension> &);

  template <class VectorType>
  void add_data_vector (const VectorType               &data,
                        const std::vector<std::string> &names,
                        const DataVectorType           type = type_automatic,
                        const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation
                        = std::vector<DataComponentInterpretation::DataComponentInterpretation>());

  template <class VectorType>
  void add_data_vector (const VectorType     &data,
                        const std::string    &name,
                        const DataVectorType type = type_automatic,
                        const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation
                        = std::vector<DataComponentInterpretation::DataComponentInterpretation>());

  template <class VectorType>
  void add_data_vector (const DoFHandlerType           &dof_handler,
                        const VectorType               &data,
                        const std::vector<std::string> &names,
                        const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation
                        = std::vector<DataComponentInterpretation::DataComponentInterpretation>());


  template <class VectorType>
  void add_data_vector (const DoFHandlerType          &dof_handler,
                        const VectorType  &data,
                        const std::string &name,
                        const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation
                        = std::vector<DataComponentInterpretation::DataComponentInterpretation>());

  template <class VectorType>
  void add_data_vector (const VectorType                             &data,
                        const DataPostprocessor<DoFHandlerType::space_dimension> &data_postprocessor);

  template <class VectorType>
  void add_data_vector (const DoFHandlerType                         &dof_handler,
                        const VectorType                             &data,
                        const DataPostprocessor<DoFHandlerType::space_dimension> &data_postprocessor);

  void clear_data_vectors ();

  void clear_input_data_references ();

  template <typename DoFHandlerType2>
  void merge_patches (const DataOut_DoFData<DoFHandlerType2,patch_dim,patch_space_dim> &source,
                      const Point<patch_space_dim> &shift = Point<patch_space_dim>());

  virtual void clear ();

  std::size_t memory_consumption () const;

protected:
  typedef ::DataOutBase::Patch<patch_dim,patch_space_dim> Patch;

  SmartPointer<const Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> > triangulation;

  SmartPointer<const DoFHandlerType> dofs;

  std::vector<std::shared_ptr<internal::DataOutImplementation::DataEntryBase<DoFHandlerType> > >  dof_data;

  std::vector<std::shared_ptr<internal::DataOutImplementation::DataEntryBase<DoFHandlerType> > >  cell_data;

  std::vector<Patch> patches;

  virtual
  const std::vector<Patch> &get_patches () const;

  virtual
  std::vector<std::string> get_dataset_names () const;

  std::vector<std::shared_ptr<::hp::FECollection<DoFHandlerType::dimension,DoFHandlerType::space_dimension> > >
  get_fes() const;

  virtual
  std::vector<std::tuple<unsigned int, unsigned int, std::string> >
  get_vector_data_ranges () const;

  template <class, int, int>
  friend class DataOut_DoFData;
private:
  template <class VectorType>
  void
  add_data_vector_internal
  (const DoFHandlerType           *dof_handler,
   const VectorType               &data,
   const std::vector<std::string> &names,
   const DataVectorType            type,
   const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation,
   const bool                      deduce_output_names);
};



// -------------------- template and inline functions ------------------------
template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename VectorType>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
add_data_vector
(const VectorType                         &vec,
 const std::string                        &name,
 const DataVectorType                      type,
 const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation)
{
  Assert (triangulation != nullptr, Exceptions::DataOutImplementation::ExcNoTriangulationSelected ());
  std::vector<std::string> names(1, name);
  add_data_vector_internal (dofs, vec, names, type, data_component_interpretation, true);
}



template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename VectorType>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
add_data_vector
(const VectorType                         &vec,
 const std::vector<std::string>           &names,
 const DataVectorType                      type,
 const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation)
{
  Assert (triangulation != nullptr, Exceptions::DataOutImplementation::ExcNoTriangulationSelected ());
  add_data_vector_internal(dofs, vec, names, type, data_component_interpretation, false);
}



template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename VectorType>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
add_data_vector
(const DoFHandlerType           &dof_handler,
 const VectorType               &data,
 const std::string              &name,
 const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation)
{
  std::vector<std::string> names(1, name);
  add_data_vector_internal (&dof_handler, data, names, type_dof_data, data_component_interpretation, true);
}



template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename VectorType>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
add_data_vector
(const DoFHandlerType           &dof_handler,
 const VectorType               &data,
 const std::vector<std::string> &names,
 const std::vector<DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretation)
{
  add_data_vector_internal(&dof_handler, data, names, type_dof_data, data_component_interpretation, false);
}



template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename VectorType>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
add_data_vector (const VectorType                       &vec,
                 const DataPostprocessor<DoFHandlerType::space_dimension> &data_postprocessor)
{
  Assert (dofs != nullptr, Exceptions::DataOutImplementation::ExcNoDoFHandlerSelected ());
  add_data_vector(*dofs, vec, data_postprocessor);
}



template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
template <typename DoFHandlerType2>
void
DataOut_DoFData<DoFHandlerType,patch_dim,patch_space_dim>::
merge_patches (const DataOut_DoFData<DoFHandlerType2,patch_dim,patch_space_dim> &source,
               const Point<patch_space_dim> &shift)
{
  const std::vector<Patch> &source_patches = source.get_patches ();
  Assert ((patches.size () != 0) &&
          (source_patches.size () != 0),
          ExcMessage ("When calling this function, both the current "
                      "object and the one being merged need to have a "
                      "nonzero number of patches associated with it. "
                      "Either you called this function on objects that "
                      "are empty, or you may have forgotten to call "
                      "the 'build_patches()' function."));
  // check equality of component
  // names
  Assert (get_dataset_names() == source.get_dataset_names(),
          Exceptions::DataOutImplementation::ExcIncompatibleDatasetNames());
  // make sure patches are compatible. we'll
  // assume that if the first respective
  // patches are ok that all the other ones
  // are ok as well
  Assert (patches[0].n_subdivisions == source_patches[0].n_subdivisions,
          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
  Assert (patches[0].data.n_rows() == source_patches[0].data.n_rows(),
          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
  Assert (patches[0].data.n_cols() == source_patches[0].data.n_cols(),
          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());

  // check equality of the vector data
  // specifications
  Assert (get_vector_data_ranges().size() ==
          source.get_vector_data_ranges().size(),
          ExcMessage ("Both sources need to declare the same components "
                      "as vectors."));
  for (unsigned int i=0; i<get_vector_data_ranges().size(); ++i)
    {
      Assert (std::get<0>(get_vector_data_ranges()[i]) ==
              std::get<0>(source.get_vector_data_ranges()[i]),
              ExcMessage ("Both sources need to declare the same components "
                          "as vectors."));
      Assert (std::get<1>(get_vector_data_ranges()[i]) ==
              std::get<1>(source.get_vector_data_ranges()[i]),
              ExcMessage ("Both sources need to declare the same components "
                          "as vectors."));
      Assert (std::get<2>(get_vector_data_ranges()[i]) ==
              std::get<2>(source.get_vector_data_ranges()[i]),
              ExcMessage ("Both sources need to declare the same components "
                          "as vectors."));
    }

  // merge patches. store old number
  // of elements, since we need to
  // adjust patch numbers, etc
  // afterwards
  const unsigned int old_n_patches = patches.size();
  patches.insert (patches.end(),
                  source_patches.begin(),
                  source_patches.end());

  // perform shift, if so desired
  if (shift != Point<patch_space_dim>())
    for (unsigned int i=old_n_patches; i<patches.size(); ++i)
      for (unsigned int v=0; v<GeometryInfo<patch_dim>::vertices_per_cell; ++v)
        patches[i].vertices[v] += shift;

  // adjust patch numbers
  for (unsigned int i=old_n_patches; i<patches.size(); ++i)
    patches[i].patch_index += old_n_patches;

  // adjust patch neighbors
  for (unsigned int i=old_n_patches; i<patches.size(); ++i)
    for (unsigned int n=0; n<GeometryInfo<patch_dim>::faces_per_cell; ++n)
      if (patches[i].neighbors[n] != Patch::no_neighbor)
        patches[i].neighbors[n] += old_n_patches;
}


DEAL_II_NAMESPACE_CLOSE

#endif