data_out_rotation.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/numerics/data_out_rotation.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


//TODO: Update documentation
//TODO: Unify code for dimensions


//TODO: build_some_patches isn't going to work if first_cell/next_cell
//don't iterate over all cells and if cell data is requested. in that
//case, we need to calculate cell_number as in the DataOut class

// Not implemented for 3D


namespace internal
{
  namespace DataOutRotationImplementation
  {
    template <int dim, int spacedim>
    ParallelData<dim,spacedim>::
    ParallelData (const unsigned int n_datasets,
                  const unsigned int n_subdivisions,
                  const unsigned int n_patches_per_circle,
                  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)
      :
      internal::DataOutImplementation::
      ParallelDataBase<dim,spacedim> (n_datasets,
                                      n_subdivisions,
                                      n_postprocessor_outputs,
                                      mapping,
                                      finite_elements,
                                      update_flags,
                                      false),
      n_patches_per_circle (n_patches_per_circle)
    {}



    template <int dim, int spacedim>
    void
    append_patch_to_list (const std::vector<DataOutBase::Patch<dim+1,spacedim+1> > &new_patches,
                          std::vector<DataOutBase::Patch<dim+1,spacedim+1> > &patches)
    {
      for (unsigned int i=0; i<new_patches.size(); ++i)
        {
          patches.push_back (new_patches[i]);
          patches.back().patch_index = patches.size()-1;
        }
    }
  }
}



template <int dim, typename DoFHandlerType>
void
DataOutRotation<dim,DoFHandlerType>::
build_one_patch (const cell_iterator                                                 *cell,
                 internal::DataOutRotationImplementation::ParallelData<dimension, space_dimension> &data,
                 std::vector<DataOutBase::Patch<dimension+1,space_dimension+1> >     &my_patches)
{
  if (dim == 3)
    {
      // would this function make any sense after all? who would want to
      // output/compute in four space dimensions?
      Assert (false, ExcNotImplemented());
      return;
    }

  Assert ((*cell)->is_locally_owned(),
          ExcNotImplemented());

  const unsigned int n_patches_per_circle = data.n_patches_per_circle;

  // another abbreviation denoting the number of q_points in each direction
  const unsigned int n_points = data.n_subdivisions+1;

  // set up an array that holds the directions in the plane of rotation in
  // which we will put points in the whole domain (not the rotationally
  // reduced one in which the computation took place. for simplicity add the
  // initial direction at the end again
  std::vector<Point<dimension+1> > angle_directions (n_patches_per_circle+1);
  for (unsigned int i=0; i<=n_patches_per_circle; ++i)
    {
      angle_directions[i][dimension-1] = std::cos(2*numbers::PI *
                                                  i/n_patches_per_circle);
      angle_directions[i][dimension] = std::sin(2*numbers::PI *
                                                i/n_patches_per_circle);
    }

  for (unsigned int angle=0; angle<n_patches_per_circle; ++angle)
    {
      // first compute the vertices of the patch. note that they will have to
      // be computed from the vertices of the cell, which has one dimension
      // less, however.
      switch (dimension)
        {
        case 1:
        {
          const double r1 = (*cell)->vertex(0)(0),
                       r2 = (*cell)->vertex(1)(0);
          Assert (r1 >= 0, ExcRadialVariableHasNegativeValues(r1));
          Assert (r2 >= 0, ExcRadialVariableHasNegativeValues(r2));

          my_patches[angle].vertices[0] = r1*angle_directions[angle];
          my_patches[angle].vertices[1] = r2*angle_directions[angle];
          my_patches[angle].vertices[2] = r1*angle_directions[angle+1];
          my_patches[angle].vertices[3] = r2*angle_directions[angle+1];

          break;
        };

        case 2:
        {
          for (unsigned int vertex=0;
               vertex<GeometryInfo<dimension>::vertices_per_cell;
               ++vertex)
            {
              const Point<dimension> v = (*cell)->vertex(vertex);

              // make sure that the radial variable is nonnegative
              Assert (v(0) >= 0, ExcRadialVariableHasNegativeValues(v(0)));

              // now set the vertices of the patch
              my_patches[angle].vertices[vertex] = v(0) * angle_directions[angle];
              my_patches[angle].vertices[vertex][0] = v(1);

              my_patches[angle].vertices[vertex+GeometryInfo<dimension>::vertices_per_cell]
                = v(0) * angle_directions[angle+1];
              my_patches[angle].vertices[vertex+GeometryInfo<dimension>::vertices_per_cell][0]
                = v(1);
            };

          break;
        };

        default:
          Assert (false, ExcNotImplemented());
        };

      // then fill in data
      if (data.n_datasets > 0)
        {
          unsigned int offset=0;

          data.reinit_all_fe_values(this->dof_data, *cell);
          // first fill dof_data
          for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
            {
              const FEValuesBase<dimension> &fe_patch_values
                = data.get_present_fe_values(dataset);
              const unsigned int n_components
                = fe_patch_values.get_fe().n_components();
              const DataPostprocessor<dim> *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 (fe_patch_values,
                                                                      internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                      data.patch_values_scalar.solution_values);
                      if (update_flags & update_gradients)
                        this->dof_data[dataset]->get_function_gradients (fe_patch_values,
                                                                         internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                         data.patch_values_scalar.solution_gradients);
                      if (update_flags & update_hessians)
                        this->dof_data[dataset]->get_function_hessians (fe_patch_values,
                                                                        internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                        data.patch_values_scalar.solution_hessians);

                      if (update_flags & update_quadrature_points)
                        data.patch_values_scalar.evaluation_points = fe_patch_values.get_quadrature_points();

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

                      postprocessor->
                      evaluate_scalar_field(data.patch_values_scalar,
                                            data.postprocessed_values[dataset]);
                    }
                  else
                    {
                      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 (fe_patch_values,
                                                                      internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                      data.patch_values_system.solution_values);
                      if (update_flags & update_gradients)
                        this->dof_data[dataset]->get_function_gradients (fe_patch_values,
                                                                         internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                         data.patch_values_system.solution_gradients);
                      if (update_flags & update_hessians)
                        this->dof_data[dataset]->get_function_hessians (fe_patch_values,
                                                                        internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                        data.patch_values_system.solution_hessians);

                      if (update_flags & update_quadrature_points)
                        data.patch_values_system.evaluation_points = fe_patch_values.get_quadrature_points();

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

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

                  for (unsigned int component=0;
                       component<this->dof_data[dataset]->n_output_variables;
                       ++component)
                    {
                      switch (dimension)
                        {
                        case 1:
                          for (unsigned int x=0; x<n_points; ++x)
                            for (unsigned int y=0; y<n_points; ++y)
                              my_patches[angle].data(offset+component,
                                                     x*n_points + y)
                                = data.postprocessed_values[dataset][x](component);
                          break;

                        case 2:
                          for (unsigned int x=0; x<n_points; ++x)
                            for (unsigned int y=0; y<n_points; ++y)
                              for (unsigned int z=0; z<n_points; ++z)
                                my_patches[angle].data(offset+component,
                                                       x*n_points*n_points +
                                                       y*n_points +
                                                       z)
                                  = data.postprocessed_values[dataset][x*n_points+z](component);
                          break;

                        default:
                          Assert (false, ExcNotImplemented());
                        }
                    }
                }
              else if (n_components == 1)
                {
                  this->dof_data[dataset]->get_function_values (fe_patch_values,
                                                                internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                data.patch_values_scalar.solution_values);

                  switch (dimension)
                    {
                    case 1:
                      for (unsigned int x=0; x<n_points; ++x)
                        for (unsigned int y=0; y<n_points; ++y)
                          my_patches[angle].data(offset,
                                                 x*n_points + y)
                            = data.patch_values_scalar.solution_values[x];
                      break;

                    case 2:
                      for (unsigned int x=0; x<n_points; ++x)
                        for (unsigned int y=0; y<n_points; ++y)
                          for (unsigned int z=0; z<n_points; ++z)
                            my_patches[angle].data(offset,
                                                   x*n_points*n_points +
                                                   y +
                                                   z*n_points)
                              = data.patch_values_scalar.solution_values[x*n_points+z];
                      break;

                    default:
                      Assert (false, ExcNotImplemented());
                    }
                }
              else
                // system of components
                {
                  data.resize_system_vectors(n_components);
                  this->dof_data[dataset]->get_function_values (fe_patch_values,
                                                                internal::DataOutImplementation::ComponentExtractor::real_part,
                                                                data.patch_values_system.solution_values);

                  for (unsigned int component=0; component<n_components;
                       ++component)
                    {
                      switch (dimension)
                        {
                        case 1:
                          for (unsigned int x=0; x<n_points; ++x)
                            for (unsigned int y=0; y<n_points; ++y)
                              my_patches[angle].data(offset+component,
                                                     x*n_points + y)
                                = data.patch_values_system.solution_values[x](component);
                          break;

                        case 2:
                          for (unsigned int x=0; x<n_points; ++x)
                            for (unsigned int y=0; y<n_points; ++y)
                              for (unsigned int z=0; z<n_points; ++z)
                                my_patches[angle].data(offset+component,
                                                       x*n_points*n_points +
                                                       y*n_points +
                                                       z)
                                  = data.patch_values_system.solution_values[x*n_points+z](component);
                          break;

                        default:
                          Assert (false, ExcNotImplemented());
                        }
                    }
                }
              offset+=this->dof_data[dataset]->n_output_variables;
            }

          // then do the cell data
          for (unsigned int dataset=0; dataset<this->cell_data.size(); ++dataset)
            {
              // we need to get at the number of the cell to which this face
              // belongs in order to access the cell data. this is not readily
              // available, so choose the following rather inefficient way:
              Assert ((*cell)->active(),
                      ExcMessage("Cell must be active for cell data"));
              const unsigned int cell_number
                = std::distance (this->triangulation->begin_active(),
                                 typename Triangulation<dimension,space_dimension>::active_cell_iterator(*cell));
              const double value
                = this->cell_data[dataset]->get_cell_data_value (cell_number,
                                                                 internal::DataOutImplementation::ComponentExtractor::real_part);
              switch (dimension)
                {
                case 1:
                  for (unsigned int x=0; x<n_points; ++x)
                    for (unsigned int y=0; y<n_points; ++y)
                      my_patches[angle].data(dataset+offset,
                                             x*n_points +
                                             y)
                        = value;
                  break;

                case 2:
                  for (unsigned int x=0; x<n_points; ++x)
                    for (unsigned int y=0; y<n_points; ++y)
                      for (unsigned int z=0; z<n_points; ++z)
                        my_patches[angle].data(dataset+offset,
                                               x*n_points*n_points +
                                               y*n_points +
                                               z)
                          = value;
                  break;

                default:
                  Assert (false, ExcNotImplemented());
                }
            }
        }
    }
}



template <int dim, typename DoFHandlerType>
void DataOutRotation<dim,DoFHandlerType>::build_patches (const unsigned int n_patches_per_circle,
                                                         const unsigned int nnnn_subdivisions)
{
  // Check consistency of redundant
  // template parameter
  Assert (dim==dimension, ExcDimensionMismatch(dim, dimension));
  Assert (this->triangulation != nullptr,
          Exceptions::DataOutImplementation::ExcNoTriangulationSelected());

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

  this->validate_dataset_names();

  unsigned int n_datasets=this->cell_data.size();
  for (unsigned int i=0; i<this->dof_data.size(); ++i)
    n_datasets+= this->dof_data[i]->n_output_variables;

  UpdateFlags update_flags=update_values | 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."));

  // first count the cells we want to
  // create patches of and make sure
  // there is enough memory for that
  std::vector<cell_iterator> all_cells;
  for (cell_iterator cell=first_cell(); cell != this->triangulation->end();
       cell = next_cell(cell))
    all_cells.push_back (cell);

  // then also take into account that
  // we want more than one patch to
  // come out of every cell, as they
  // are repeated around the axis of
  // rotation
  this->patches.clear();
  this->patches.reserve (all_cells.size() * n_patches_per_circle);


  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;

  internal::DataOutRotationImplementation::ParallelData<dimension, space_dimension>
  thread_data (n_datasets,
               n_subdivisions, n_patches_per_circle,
               n_postprocessor_outputs,
               StaticMappingQ1<dimension,space_dimension>::mapping,
               this->get_fes(),
               update_flags);
  std::vector<DataOutBase::Patch<dimension+1,space_dimension+1> >
  new_patches (n_patches_per_circle);
  for (unsigned int i=0; i<new_patches.size(); ++i)
    {
      new_patches[i].n_subdivisions = n_subdivisions;
      new_patches[i].data.reinit (n_datasets,
                                  Utilities::fixed_power<dimension+1>(n_subdivisions+1));
    }

  // now build the patches in parallel
  WorkStream::run (&all_cells[0],
                   &all_cells[0]+all_cells.size(),
                   std::bind(&DataOutRotation<dim,DoFHandlerType>::build_one_patch,
                             this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3),
                   std::bind(&internal::DataOutRotationImplementation
                             ::append_patch_to_list<dim,space_dimension>,
                             std::placeholders::_1, std::ref(this->patches)),
                   thread_data,
                   new_patches);
}



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


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



// explicit instantiations
#include "data_out_rotation.inst"


DEAL_II_NAMESPACE_CLOSE