data_out_dof_data.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 1999 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#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/mg_level_object.h>
#include <deal.II/base/observer_pointer.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. This is "
                     "generally done using the 'attach_dof_handler()' or "
                     "'attach_triangulation()' member functions.");
 
    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. This is "
                     "generally done using the 'attach_dof_handler()' "
                     "member function.");
 
    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 (when"
                   << " the data is of type type_dof_data), or equal to the"
                   << " number of active cells (when the data is of type "
                   << " type_cell_data).");
    DeclException2(
      ExcInvalidCharacter,
      std::string,
      std::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 <int dim, int spacedim>
    class DataEntryBase
    {
    public:
      DataEntryBase(const DoFHandler<dim, spacedim> *dofs,
                    const std::vector<std::string>  &names,
                    const std::vector<
                      DataComponentInterpretation::DataComponentInterpretation>
                      &data_component_interpretation);
 
      DataEntryBase(const DoFHandler<dim, spacedim>   *dofs,
                    const DataPostprocessor<spacedim> *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<dim, spacedim> &fe_patch_values,
                          const ComponentExtractor           extract_component,
                          std::vector<double> &patch_values) const = 0;
 
      virtual void
      get_function_values(
        const FEValuesBase<dim, spacedim>   &fe_patch_values,
        const ComponentExtractor             extract_component,
        std::vector<::Vector<double>> &patch_values_system) const = 0;
 
      virtual void
      get_function_gradients(
        const FEValuesBase<dim, spacedim> &fe_patch_values,
        const ComponentExtractor           extract_component,
        std::vector<Tensor<1, spacedim>>  &patch_gradients) const = 0;
 
      virtual void
      get_function_gradients(const FEValuesBase<dim, spacedim> &fe_patch_values,
                             const ComponentExtractor extract_component,
                             std::vector<std::vector<Tensor<1, spacedim>>>
                               &patch_gradients_system) const = 0;
 
      virtual void
      get_function_hessians(
        const FEValuesBase<dim, spacedim> &fe_patch_values,
        const ComponentExtractor           extract_component,
        std::vector<Tensor<2, spacedim>>  &patch_hessians) const = 0;
 
      virtual void
      get_function_hessians(const FEValuesBase<dim, spacedim> &fe_patch_values,
                            const ComponentExtractor extract_component,
                            std::vector<std::vector<Tensor<2, spacedim>>>
                              &patch_hessians_system) const = 0;
 
      virtual bool
      is_complex_valued() const = 0;
 
      virtual void
      clear() = 0;
 
      virtual std::size_t
      memory_consumption() const = 0;
 
      ObserverPointer<const DoFHandler<dim, spacedim>> dof_handler;
 
      const std::vector<std::string> names;
 
      const std::vector<
        DataComponentInterpretation::DataComponentInterpretation>
        data_component_interpretation;
 
      ObserverPointer<const ::DataPostprocessor<spacedim>> 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 unsigned int               n_datasets,
        const unsigned int               n_subdivisions,
        const std::vector<unsigned int> &n_postprocessor_outputs,
        const ::hp::MappingCollection<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);
 
      void
      reinit_all_fe_values(
        std::vector<std::shared_ptr<DataEntryBase<dim, spacedim>>> &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 <int dim,
          int patch_dim,
          int spacedim       = dim,
          int patch_spacedim = patch_dim>
class DataOut_DoFData : public DataOutInterface<patch_dim, patch_spacedim>
{
public:
  using cell_iterator = typename Triangulation<dim, spacedim>::cell_iterator;
 
public:
  enum DataVectorType
  {
    type_dof_data,
 
    type_cell_data,
 
    type_automatic
  };
 
  DataOut_DoFData();
 
  virtual ~DataOut_DoFData() override;
 
  void
  attach_dof_handler(const DoFHandler<dim, spacedim> &);
 
  void
  attach_triangulation(const Triangulation<dim, spacedim> &);
 
  template <typename 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 = {});
 
  template <typename 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 = {});
 
  template <typename VectorType>
  void
  add_data_vector(
    const DoFHandler<dim, spacedim> &dof_handler,
    const VectorType                &data,
    const std::vector<std::string>  &names,
    const std::vector<DataComponentInterpretation::DataComponentInterpretation>
      &data_component_interpretation = {});
 
 
  template <typename VectorType>
  void
  add_data_vector(
    const DoFHandler<dim, spacedim> &dof_handler,
    const VectorType                &data,
    const std::string               &name,
    const std::vector<DataComponentInterpretation::DataComponentInterpretation>
      &data_component_interpretation = {});
 
  template <typename VectorType>
  void
  add_data_vector(const VectorType                  &data,
                  const DataPostprocessor<spacedim> &data_postprocessor);
 
  template <typename VectorType>
  void
  add_data_vector(const DoFHandler<dim, spacedim>   &dof_handler,
                  const VectorType                  &data,
                  const DataPostprocessor<spacedim> &data_postprocessor);
 
  template <typename VectorType>
  void
  add_mg_data_vector(
    const DoFHandler<dim, spacedim> &dof_handler,
    const MGLevelObject<VectorType> &data,
    const std::vector<std::string>  &names,
    const std::vector<DataComponentInterpretation::DataComponentInterpretation>
      &data_component_interpretation = std::vector<
        DataComponentInterpretation::DataComponentInterpretation>());
 
  template <typename VectorType>
  void
  add_mg_data_vector(const DoFHandler<dim, spacedim> &dof_handler,
                     const MGLevelObject<VectorType> &data,
                     const std::string               &name);
 
  void
  clear_data_vectors();
 
  void
  clear_input_data_references();
 
  template <int dim2, int spacedim2>
  void
  merge_patches(
    const DataOut_DoFData<dim2, patch_dim, spacedim2, patch_spacedim> &source,
    const Point<patch_spacedim> &shift = Point<patch_spacedim>());
 
  virtual void
  clear();
 
  std::size_t
  memory_consumption() const;
 
protected:
  using Patch = ::DataOutBase::Patch<patch_dim, patch_spacedim>;
 
  ObserverPointer<const Triangulation<dim, spacedim>> triangulation;
 
  ObserverPointer<const DoFHandler<dim, spacedim>> dofs;
 
  std::vector<std::shared_ptr<
    internal::DataOutImplementation::DataEntryBase<dim, spacedim>>>
    dof_data;
 
  std::vector<std::shared_ptr<
    internal::DataOutImplementation::DataEntryBase<dim, spacedim>>>
    cell_data;
 
  std::vector<Patch> patches;
 
public:
  virtual const std::vector<Patch> &
  get_patches() const override;
 
protected:
  virtual std::vector<std::string>
  get_dataset_names() const override;
 
  virtual std::vector<
    std::tuple<unsigned int,
               unsigned int,
               std::string,
               DataComponentInterpretation::DataComponentInterpretation>>
  get_nonscalar_data_ranges() const override;
 
protected:
  void
  remove_ignored_components();
 
  std::vector<std::shared_ptr<::hp::FECollection<dim, spacedim>>>
  get_fes() const;
 
  // Make all template siblings friends. Needed for the merge_patches()
  // function.
  template <int, int, int, int>
  friend class DataOut_DoFData;
 
  template <int, class>
  friend class MGDataOut;
 
private:
  template <typename VectorType>
  void
  add_data_vector_internal(
    const DoFHandler<dim, spacedim> *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 <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <typename VectorType>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::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 <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <typename VectorType>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::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 <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <typename VectorType>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::add_data_vector(
  const DoFHandler<dim, spacedim> &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 <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <typename VectorType>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::add_data_vector(
  const DoFHandler<dim, spacedim> &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 <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <typename VectorType>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::add_data_vector(
  const VectorType                  &vec,
  const DataPostprocessor<spacedim> &data_postprocessor)
{
  Assert(dofs != nullptr,
         Exceptions::DataOutImplementation::ExcNoDoFHandlerSelected());
  add_data_vector(*dofs, vec, data_postprocessor);
}
 
 
 
template <int dim, int patch_dim, int spacedim, int patch_spacedim>
template <int dim2, int spacedim2>
void
DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::merge_patches(
  const DataOut_DoFData<dim2, patch_dim, spacedim2, patch_spacedim> &source,
  const Point<patch_spacedim>                                       &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. Ideally, we would check that all input
  // patches from both collections are all compatible, but we'll be content
  // with checking that just the first ones from both sources are.
  //
  // We check compatibility by testing that both sets of patches result
  // from the same number of subdivisions, and that they have the same
  // number of source vectors (they really should, since we already checked
  // that there are the same number of source components above, but you
  // never know). This implies that the data should have the same number of
  // columns. They should really have the same number of rows as well,
  // but depending on whether a patch has points included or not, the
  // number of rows may or may not include coordinates for the points,
  // and the comparison has to account for that because in each source
  // stream, the patches may include some that have points included.
  Assert(patches[0].n_subdivisions == source_patches[0].n_subdivisions,
         Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
  Assert(patches[0].data.n_cols() == source_patches[0].data.n_cols(),
         Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
  Assert((patches[0].data.n_rows() +
          (patches[0].points_are_available ? 0 : patch_spacedim)) ==
           (source_patches[0].data.n_rows() +
            (source_patches[0].points_are_available ? 0 : patch_spacedim)),
         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_spacedim>())
    for (unsigned int i = old_n_patches; i < patches.size(); ++i)
      for (const unsigned int v : GeometryInfo<patch_dim>::vertex_indices())
        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 (const unsigned int n : GeometryInfo<patch_dim>::face_indices())
      if (patches[i].neighbors[n] != Patch::no_neighbor)
        patches[i].neighbors[n] += old_n_patches;
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif