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

#ifndef dealii_data_out_dof_data_h
#define dealii_data_out_dof_data_h



#include <deal.II/base/config.h>

#include <deal.II/base/data_out_base.h>
#include <deal.II/base/smartpointer.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/mapping.h>

#include <deal.II/grid/tria.h>

#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/q_collection.h>

#include <deal.II/numerics/data_component_interpretation.h>
#include <deal.II/numerics/data_postprocessor.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.");

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

  } // namespace DataOutImplementation
} // namespace Exceptions


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;
    };
  } // namespace DataOutImplementation
} // namespace internal


// 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:
  using cell_iterator =
    typename Triangulation<DoFHandlerType::dimension,
                           DoFHandlerType::space_dimension>::cell_iterator;
  using active_cell_iterator = typename Triangulation<
    DoFHandlerType::dimension,
    DoFHandlerType::space_dimension>::active_cell_iterator;

public:
  enum DataVectorType
  {
    type_dof_data,

    type_cell_data,

    type_automatic
  };

  DataOut_DoFData();

  virtual ~DataOut_DoFData() override;

  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:
  using Patch = ::DataOutBase::Patch<patch_dim, patch_space_dim>;

  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 override;

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

  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,
               DataComponentInterpretation::DataComponentInterpretation>>
  get_nonscalar_data_ranges() const override;

  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_nonscalar_data_ranges().size() ==
           source.get_nonscalar_data_ranges().size(),
         ExcMessage("Both sources need to declare the same components "
                    "as vectors."));
  for (unsigned int i = 0; i < get_nonscalar_data_ranges().size(); ++i)
    {
      Assert(std::get<0>(get_nonscalar_data_ranges()[i]) ==
               std::get<0>(source.get_nonscalar_data_ranges()[i]),
             ExcMessage("Both sources need to declare the same components "
                        "as vectors."));
      Assert(std::get<1>(get_nonscalar_data_ranges()[i]) ==
               std::get<1>(source.get_nonscalar_data_ranges()[i]),
             ExcMessage("Both sources need to declare the same components "
                        "as vectors."));
      Assert(std::get<2>(get_nonscalar_data_ranges()[i]) ==
               std::get<2>(source.get_nonscalar_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