data_out.cc

// ---------------------------------------------------------------------
//
// 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.
//
// ---------------------------------------------------------------------

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

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

#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>

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

#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_values.h>

#include <deal.II/numerics/data_out.h>

#include <sstream>

DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace DataOutImplementation
  {
    template <int dim, int spacedim>
    ParallelData<dim, spacedim>::ParallelData(
      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 std::vector<std::vector<unsigned int>> &cell_to_patch_index_map)
      : ParallelDataBase<dim, spacedim>(n_datasets,
                                        n_subdivisions,
                                        n_postprocessor_outputs,
                                        mapping,
                                        finite_elements,
                                        update_flags,
                                        false)
      , cell_to_patch_index_map(&cell_to_patch_index_map)
    {}
  } // namespace DataOutImplementation
} // namespace internal



template <int dim, typename DoFHandlerType>
void
DataOut<dim, DoFHandlerType>::build_one_patch(
  const std::pair<cell_iterator, unsigned int> *cell_and_index,
  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
                                                DoFHandlerType::space_dimension>
    &                    scratch_data,
  const unsigned int     n_subdivisions,
  const CurvedCellRegion curved_cell_region)
{
  // first create the output object that we will write into
  ::DataOutBase::Patch<DoFHandlerType::dimension,
                               DoFHandlerType::space_dimension>
    patch;
  patch.n_subdivisions = n_subdivisions;

  // set the vertices of the patch. if the mapping does not preserve locations
  // (e.g. MappingQEulerian), we need to compute the offset of the vertex for
  // the graphical output. Otherwise, we can just use the vertex info.
  for (unsigned int vertex = 0;
       vertex < GeometryInfo<DoFHandlerType::dimension>::vertices_per_cell;
       ++vertex)
    if (scratch_data.mapping_collection[0].preserves_vertex_locations())
      patch.vertices[vertex] = cell_and_index->first->vertex(vertex);
    else
      patch.vertices[vertex] =
        scratch_data.mapping_collection[0].transform_unit_to_real_cell(
          cell_and_index->first,
          GeometryInfo<DoFHandlerType::dimension>::unit_cell_vertex(vertex));

  // create DoFHandlerType::active_cell_iterator and initialize FEValues
  scratch_data.reinit_all_fe_values(this->dof_data, cell_and_index->first);

  const FEValuesBase<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
    &fe_patch_values = scratch_data.get_present_fe_values(0);

  const unsigned int n_q_points = fe_patch_values.n_quadrature_points;

  // depending on the requested output of curved cells, if necessary
  // append the quadrature points to the last rows of the patch.data
  // member. This is the case if we want to produce curved cells at the
  // boundary and this cell actually is at the boundary, or else if we
  // want to produce curved cells everywhere
  //
  // note: a cell is *always* at the boundary if dim<spacedim
  if (curved_cell_region == curved_inner_cells ||
      (curved_cell_region == curved_boundary &&
       (cell_and_index->first->at_boundary() ||
        (DoFHandlerType::dimension != DoFHandlerType::space_dimension))))
    {
      Assert(patch.space_dim == DoFHandlerType::space_dimension,
             ExcInternalError());

      // set the flag indicating that for this cell the points are
      // explicitly given
      patch.points_are_available = true;

      // then resize the patch.data member in order to have enough memory for
      // the quadrature points as well, and copy the quadrature points there
      const std::vector<Point<DoFHandlerType::space_dimension>> &q_points =
        fe_patch_values.get_quadrature_points();

      patch.data.reinit(scratch_data.n_datasets +
                          DoFHandlerType::space_dimension,
                        n_q_points);
      for (unsigned int i = 0; i < DoFHandlerType::space_dimension; ++i)
        for (unsigned int q = 0; q < n_q_points; ++q)
          patch.data(patch.data.size(0) - DoFHandlerType::space_dimension + i,
                     q) = q_points[q][i];
    }
  else
    {
      patch.data.reinit(scratch_data.n_datasets, n_q_points);
      patch.points_are_available = false;
    }


  if (scratch_data.n_datasets > 0)
    {
      // counter for data records
      unsigned int offset = 0;

      // first fill dof_data
      for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
        {
          const FEValuesBase<DoFHandlerType::dimension,
                             DoFHandlerType::space_dimension>
            &this_fe_patch_values = scratch_data.get_present_fe_values(dataset);
          const unsigned int n_components =
            this_fe_patch_values.get_fe().n_components();

          const DataPostprocessor<DoFHandlerType::space_dimension>
            *postprocessor = this->dof_data[dataset]->postprocessor;

          if (postprocessor != nullptr)
            {
              // we have to postprocess the data, so determine, which fields
              // have to be updated
              const UpdateFlags update_flags =
                postprocessor->get_needed_update_flags();
              if (n_components == 1)
                {
                  // at each point there is only one component of value,
                  // gradient etc.
                  if (update_flags & update_values)
                    this->dof_data[dataset]->get_function_values(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_scalar.solution_values);
                  if (update_flags & update_gradients)
                    this->dof_data[dataset]->get_function_gradients(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_scalar.solution_gradients);
                  if (update_flags & update_hessians)
                    this->dof_data[dataset]->get_function_hessians(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_scalar.solution_hessians);

                  if (update_flags & update_quadrature_points)
                    scratch_data.patch_values_scalar.evaluation_points =
                      this_fe_patch_values.get_quadrature_points();

                  const typename DoFHandlerType::active_cell_iterator dh_cell(
                    &cell_and_index->first->get_triangulation(),
                    cell_and_index->first->level(),
                    cell_and_index->first->index(),
                    this->dof_data[dataset]->dof_handler);
                  scratch_data.patch_values_scalar
                    .template set_cell<DoFHandlerType>(dh_cell);

                  postprocessor->evaluate_scalar_field(
                    scratch_data.patch_values_scalar,
                    scratch_data.postprocessed_values[dataset]);
                }
              else
                {
                  scratch_data.resize_system_vectors(n_components);

                  // at each point there is a vector valued function and its
                  // derivative...
                  if (update_flags & update_values)
                    this->dof_data[dataset]->get_function_values(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_system.solution_values);
                  if (update_flags & update_gradients)
                    this->dof_data[dataset]->get_function_gradients(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_system.solution_gradients);
                  if (update_flags & update_hessians)
                    this->dof_data[dataset]->get_function_hessians(
                      this_fe_patch_values,
                      internal::DataOutImplementation::ComponentExtractor::
                        real_part,
                      scratch_data.patch_values_system.solution_hessians);

                  if (update_flags & update_quadrature_points)
                    scratch_data.patch_values_system.evaluation_points =
                      this_fe_patch_values.get_quadrature_points();

                  const typename DoFHandlerType::active_cell_iterator dh_cell(
                    &cell_and_index->first->get_triangulation(),
                    cell_and_index->first->level(),
                    cell_and_index->first->index(),
                    this->dof_data[dataset]->dof_handler);
                  scratch_data.patch_values_system
                    .template set_cell<DoFHandlerType>(dh_cell);

                  postprocessor->evaluate_vector_field(
                    scratch_data.patch_values_system,
                    scratch_data.postprocessed_values[dataset]);
                }

              for (unsigned int q = 0; q < n_q_points; ++q)
                for (unsigned int component = 0;
                     component < this->dof_data[dataset]->n_output_variables;
                     ++component)
                  patch.data(offset + component, q) =
                    scratch_data.postprocessed_values[dataset][q](component);
            }
          else
            // use the given data vector directly, without a postprocessor.
            // again, we treat single component functions separately for
            // efficiency reasons.
            if (n_components == 1)
            {
              // first output the real part of the solution vector
              this->dof_data[dataset]->get_function_values(
                this_fe_patch_values,
                internal::DataOutImplementation::ComponentExtractor::real_part,
                scratch_data.patch_values_scalar.solution_values);
              for (unsigned int q = 0; q < n_q_points; ++q)
                patch.data(offset, q) =
                  scratch_data.patch_values_scalar.solution_values[q];

              // and if there is one, also output the imaginary part
              if (this->dof_data[dataset]->is_complex_valued() == true)
                {
                  this->dof_data[dataset]->get_function_values(
                    this_fe_patch_values,
                    internal::DataOutImplementation::ComponentExtractor::
                      imaginary_part,
                    scratch_data.patch_values_scalar.solution_values);
                  for (unsigned int q = 0; q < n_q_points; ++q)
                    patch.data(offset + 1, q) =
                      scratch_data.patch_values_scalar.solution_values[q];
                }
            }
          else
            {
              scratch_data.resize_system_vectors(n_components);

              // same as above: first the real part
              const unsigned int stride =
                (this->dof_data[dataset]->is_complex_valued() ? 2 : 1);
              this->dof_data[dataset]->get_function_values(
                this_fe_patch_values,
                internal::DataOutImplementation::ComponentExtractor::real_part,
                scratch_data.patch_values_system.solution_values);
              for (unsigned int component = 0; component < n_components;
                   ++component)
                for (unsigned int q = 0; q < n_q_points; ++q)
                  patch.data(offset + component * stride, q) =
                    scratch_data.patch_values_system.solution_values[q](
                      component);

              // and if there is one, also output the imaginary part
              if (this->dof_data[dataset]->is_complex_valued() == true)
                {
                  this->dof_data[dataset]->get_function_values(
                    this_fe_patch_values,
                    internal::DataOutImplementation::ComponentExtractor::
                      imaginary_part,
                    scratch_data.patch_values_system.solution_values);
                  for (unsigned int component = 0; component < n_components;
                       ++component)
                    for (unsigned int q = 0; q < n_q_points; ++q)
                      patch.data(offset + component * stride + 1, q) =
                        scratch_data.patch_values_system.solution_values[q](
                          component);
                }
            }

          // increment the counter for the actual data record
          offset += this->dof_data[dataset]->n_output_variables *
                    (this->dof_data[dataset]->is_complex_valued() ? 2 : 1);
        }

      // then do the cell data. only compute the number of a cell if needed;
      // also make sure that we only access cell data if the
      // first_cell/next_cell functions only return active cells
      if (this->cell_data.size() != 0)
        {
          Assert(!cell_and_index->first->has_children(), ExcNotImplemented());

          for (unsigned int dataset = 0; dataset < this->cell_data.size();
               ++dataset)
            {
              // as above, first output the real part
              {
                const double value =
                  this->cell_data[dataset]->get_cell_data_value(
                    cell_and_index->second,
                    internal::DataOutImplementation::ComponentExtractor::
                      real_part);
                for (unsigned int q = 0; q < n_q_points; ++q)
                  patch.data(offset, q) = value;
              }

              // and if there is one, also output the imaginary part
              if (this->cell_data[dataset]->is_complex_valued() == true)
                {
                  const double value =
                    this->cell_data[dataset]->get_cell_data_value(
                      cell_and_index->second,
                      internal::DataOutImplementation::ComponentExtractor::
                        imaginary_part);
                  for (unsigned int q = 0; q < n_q_points; ++q)
                    patch.data(offset + 1, q) = value;
                }

              offset += (this->cell_data[dataset]->is_complex_valued() ? 2 : 1);
            }
        }
    }


  for (unsigned int f = 0;
       f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
       ++f)
    {
      // let's look up whether the neighbor behind that face is noted in the
      // table of cells which we treat. this can only happen if the neighbor
      // exists, and is on the same level as this cell, but it may also happen
      // that the neighbor is not a member of the range of cells over which we
      // loop, in which case the respective entry in the
      // cell_to_patch_index_map will have the value no_neighbor. (note that
      // since we allocated only as much space in this array as the maximum
      // index of the cells we loop over, not every neighbor may have its
      // space in it, so we have to assume that it is extended by values
      // no_neighbor)
      if (cell_and_index->first->at_boundary(f) ||
          (cell_and_index->first->neighbor(f)->level() !=
           cell_and_index->first->level()))
        {
          patch.neighbors[f] = numbers::invalid_unsigned_int;
          continue;
        }

      const cell_iterator neighbor = cell_and_index->first->neighbor(f);
      Assert(static_cast<unsigned int>(neighbor->level()) <
               scratch_data.cell_to_patch_index_map->size(),
             ExcInternalError());
      if ((static_cast<unsigned int>(neighbor->index()) >=
           (*scratch_data.cell_to_patch_index_map)[neighbor->level()].size()) ||
          ((*scratch_data.cell_to_patch_index_map)[neighbor->level()]
                                                  [neighbor->index()] ==
           ::DataOutBase::Patch<DoFHandlerType::dimension>::no_neighbor))
        {
          patch.neighbors[f] = numbers::invalid_unsigned_int;
          continue;
        }

      // now, there is a neighbor, so get its patch number and set it for the
      // neighbor index
      patch.neighbors[f] =
        (*scratch_data
            .cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
    }

  const unsigned int patch_idx =
    (*scratch_data.cell_to_patch_index_map)[cell_and_index->first->level()]
                                           [cell_and_index->first->index()];
  // did we mess up the indices?
  Assert(patch_idx < this->patches.size(), ExcInternalError());
  patch.patch_index = patch_idx;

  // Put the patch into the patches vector. instead of copying the data,
  // simply swap the contents to avoid the penalty of writing into another
  // processor's memory
  this->patches[patch_idx].swap(patch);
}



template <int dim, typename DoFHandlerType>
void
DataOut<dim, DoFHandlerType>::build_patches(const unsigned int n_subdivisions)
{
  build_patches(StaticMappingQ1<DoFHandlerType::dimension,
                                DoFHandlerType::space_dimension>::mapping,
                n_subdivisions,
                no_curved_cells);
}



template <int dim, typename DoFHandlerType>
void
DataOut<dim, DoFHandlerType>::build_patches(
  const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
    &                    mapping,
  const unsigned int     n_subdivisions_,
  const CurvedCellRegion curved_region)
{
  // Check consistency of redundant template parameter
  Assert(dim == DoFHandlerType::dimension,
         ExcDimensionMismatch(dim, DoFHandlerType::dimension));

  Assert(this->triangulation != nullptr,
         Exceptions::DataOutImplementation::ExcNoTriangulationSelected());

  const unsigned int n_subdivisions =
    (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
  Assert(n_subdivisions >= 1,
         Exceptions::DataOutImplementation::ExcInvalidNumberOfSubdivisions(
           n_subdivisions));

  this->validate_dataset_names();

  // First count the cells we want to create patches of. Also fill the object
  // that maps the cell indices to the patch numbers, as this will be needed
  // for generation of neighborship information.
  // Note, there is a confusing mess of different indices here at play:
  // patch_index - the index of a patch in all_cells
  // cell->index - only unique on each level, used in cell_to_patch_index_map
  // active_index - index for a cell when counting from begin_active() using
  // ++cell cell_index - unique index of a cell counted using
  // next_locally_owned_cell()
  //              starting from first_locally_owned_cell()
  //
  // It turns out that we create one patch for each selected cell, so
  // patch_index==cell_index.
  //
  // will be cell_to_patch_index_map[cell->level][cell->index] = patch_index
  std::vector<std::vector<unsigned int>> cell_to_patch_index_map;
  cell_to_patch_index_map.resize(this->triangulation->n_levels());
  for (unsigned int l = 0; l < this->triangulation->n_levels(); ++l)
    {
      // max_index is the largest cell->index on level l
      unsigned int max_index = 0;
      for (cell_iterator cell = first_locally_owned_cell();
           cell != this->triangulation->end();
           cell = next_locally_owned_cell(cell))
        if (static_cast<unsigned int>(cell->level()) == l)
          max_index =
            std::max(max_index, static_cast<unsigned int>(cell->index()));

      cell_to_patch_index_map[l].resize(
        max_index + 1,
        ::DataOutBase::Patch<
          DoFHandlerType::dimension,
          DoFHandlerType::space_dimension>::no_neighbor);
    }

  // will be all_cells[patch_index] = pair(cell, active_index)
  std::vector<std::pair<cell_iterator, unsigned int>> all_cells;
  {
    // important: we need to compute the active_index of the cell in the range
    // 0..n_active_cells() because this is where we need to look up cell
    // data from (cell data vectors do not have the length distance computed by
    // first_locally_owned_cell/next_locally_owned_cell because this might skip
    // some values (FilteredIterator).
    active_cell_iterator active_cell  = this->triangulation->begin_active();
    unsigned int         active_index = 0;
    cell_iterator        cell         = first_locally_owned_cell();
    for (; cell != this->triangulation->end();
         cell = next_locally_owned_cell(cell))
      {
        // move forward until active_cell points at the cell (cell) we are
        // looking at to compute the current active_index
        while (active_cell != this->triangulation->end() && cell->active() &&
               active_cell_iterator(cell) != active_cell)
          {
            ++active_cell;
            ++active_index;
          }

        Assert(static_cast<unsigned int>(cell->level()) <
                 cell_to_patch_index_map.size(),
               ExcInternalError());
        Assert(static_cast<unsigned int>(cell->index()) <
                 cell_to_patch_index_map[cell->level()].size(),
               ExcInternalError());
        Assert(active_index < this->triangulation->n_active_cells(),
               ExcInternalError());
        cell_to_patch_index_map[cell->level()][cell->index()] =
          all_cells.size();

        all_cells.emplace_back(cell, active_index);
      }
  }

  this->patches.clear();
  this->patches.resize(all_cells.size());

  // now create a default object for the WorkStream object to work with
  unsigned int n_datasets = 0;
  for (unsigned int i = 0; i < this->cell_data.size(); ++i)
    n_datasets += (this->cell_data[i]->is_complex_valued() ? 2 : 1);
  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
    n_datasets += (this->dof_data[i]->n_output_variables *
                   (this->dof_data[i]->is_complex_valued() ? 2 : 1));

  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
    if (this->dof_data[dataset]->postprocessor)
      n_postprocessor_outputs[dataset] =
        this->dof_data[dataset]->n_output_variables;
    else
      n_postprocessor_outputs[dataset] = 0;

  const CurvedCellRegion curved_cell_region =
    (n_subdivisions < 2 ? no_curved_cells : curved_region);

  UpdateFlags update_flags = update_values;
  if (curved_cell_region != no_curved_cells)
    update_flags |= update_quadrature_points;

  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
    if (this->dof_data[i]->postprocessor)
      update_flags |=
        this->dof_data[i]->postprocessor->get_needed_update_flags();
  // perhaps update_normal_vectors is present, which would only be useful on
  // faces, but we may not use it here.
  Assert(
    !(update_flags & update_normal_vectors),
    ExcMessage(
      "The update of normal vectors may not be requested for evaluation of "
      "data on cells via DataPostprocessor."));

  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
                                                DoFHandlerType::space_dimension>
    thread_data(n_datasets,
                n_subdivisions,
                n_postprocessor_outputs,
                mapping,
                this->get_fes(),
                update_flags,
                cell_to_patch_index_map);

  // now build the patches in parallel
  if (all_cells.size() > 0)
    WorkStream::run(
      &all_cells[0],
      &all_cells[0] + all_cells.size(),
      std::bind(&DataOut<dim, DoFHandlerType>::build_one_patch,
                this,
                std::placeholders::_1,
                std::placeholders::_2,
                /* no std::placeholders::_3, since this function doesn't
                   actually need a copy data object -- it just writes everything
                   right into the output array */
                n_subdivisions,
                curved_cell_region),
      // no copy-local-to-global function needed here
      std::function<void(const int &)>(),
      thread_data,
      /* dummy CopyData object = */ 0,
      // experimenting shows that we can make things run a bit
      // faster if we increase the number of cells we work on
      // per item (i.e., WorkStream's chunk_size argument,
      // about 10% improvement) and the items in flight at any
      // given time (another 5% on the testcase discussed in
      // @ref workstream_paper, on 32 cores) and if
      8 * MultithreadInfo::n_threads(),
      64);
}



template <int dim, typename DoFHandlerType>
typename DataOut<dim, DoFHandlerType>::cell_iterator
DataOut<dim, DoFHandlerType>::first_cell()
{
  return this->triangulation->begin_active();
}



template <int dim, typename DoFHandlerType>
typename DataOut<dim, DoFHandlerType>::cell_iterator
DataOut<dim, DoFHandlerType>::next_cell(
  const typename DataOut<dim, DoFHandlerType>::cell_iterator &cell)
{
  // convert the iterator to an active_iterator and advance this to the next
  // active cell
  typename Triangulation<DoFHandlerType::dimension,
                         DoFHandlerType::space_dimension>::active_cell_iterator
    active_cell = cell;
  ++active_cell;
  return active_cell;
}



template <int dim, typename DoFHandlerType>
typename DataOut<dim, DoFHandlerType>::cell_iterator
DataOut<dim, DoFHandlerType>::first_locally_owned_cell()
{
  typename DataOut<dim, DoFHandlerType>::cell_iterator cell = first_cell();

  // skip cells if the current one has no children (is active) and is a ghost
  // or artificial cell
  while ((cell != this->triangulation->end()) &&
         (cell->has_children() == false) && !cell->is_locally_owned())
    cell = next_cell(cell);

  return cell;
}



template <int dim, typename DoFHandlerType>
typename DataOut<dim, DoFHandlerType>::cell_iterator
DataOut<dim, DoFHandlerType>::next_locally_owned_cell(
  const typename DataOut<dim, DoFHandlerType>::cell_iterator &old_cell)
{
  typename DataOut<dim, DoFHandlerType>::cell_iterator cell =
    next_cell(old_cell);
  while ((cell != this->triangulation->end()) &&
         (cell->has_children() == false) && !cell->is_locally_owned())
    cell = next_cell(cell);
  return cell;
}


// explicit instantiations
#include "data_out.inst"

DEAL_II_NAMESPACE_CLOSE