fe_series.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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_fe_series_H
#define dealii_fe_series_H



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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.h>
#include <deal.II/base/table_indices.h>
#include <deal.II/base/tensor.h>

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

#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>

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

#include <memory>
#include <string>
#include <vector>


DEAL_II_NAMESPACE_OPEN




namespace FESeries
{
  template <int dim, int spacedim = dim>
  class Fourier : public Subscriptor
  {
  public:
    using CoefficientType = typename std::complex<double>;

    Fourier(const unsigned int                     size_in_each_direction,
            const hp::FECollection<dim, spacedim> &fe_collection,
            const hp::QCollection<dim> &           q_collection);

    template <typename Number>
    void
    calculate(const ::Vector<Number> &local_dof_values,
              const unsigned int            cell_active_fe_index,
              Table<dim, CoefficientType> & fourier_coefficients);

  private:
    SmartPointer<const hp::FECollection<dim, spacedim>> fe_collection;

    SmartPointer<const hp::QCollection<dim>> q_collection;

    Table<dim, Tensor<1, dim>> k_vectors;

    std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;

    std::vector<CoefficientType> unrolled_coefficients;
  };

  template <int dim, int spacedim = dim>
  class Legendre : public Subscriptor
  {
  public:
    using CoefficientType = double;

    Legendre(const unsigned int                     size_in_each_direction,
             const hp::FECollection<dim, spacedim> &fe_collection,
             const hp::QCollection<dim> &           q_collection);

    template <typename Number>
    void
    calculate(const ::Vector<Number> &local_dof_values,
              const unsigned int            cell_active_fe_index,
              Table<dim, CoefficientType> & legendre_coefficients);

  private:
    const unsigned int N;

    SmartPointer<const hp::FECollection<dim, spacedim>> fe_collection;

    SmartPointer<const hp::QCollection<dim>> q_collection;

    std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;

    std::vector<CoefficientType> unrolled_coefficients;
  };


  template <int dim, typename CoefficientType>
  std::pair<std::vector<unsigned int>, std::vector<double>>
  process_coefficients(const Table<dim, CoefficientType> &coefficients,
                       const std::function<std::pair<bool, unsigned int>(
                         const TableIndices<dim> &)> &    predicate,
                       const VectorTools::NormType        norm);



  std::pair<double, double>
  linear_regression(const std::vector<double> &x, const std::vector<double> &y);

} // namespace FESeries

#ifndef DOXYGEN

// -------------------  inline and template functions ----------------

namespace internal
{
  namespace FESeriesImplementation
  {
    template <int dim, typename CoefficientType>
    void
    fill_map_index(
      const Table<dim, CoefficientType> &coefficients,
      const TableIndices<dim> &          ind,
      const std::function<
        std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
      std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
    {
      const std::pair<bool, unsigned int> pred_pair = predicate(ind);
      // don't add a value if predicate is false
      if (pred_pair.first == false)
        return;

      const unsigned int     pred_value  = pred_pair.second;
      const CoefficientType &coeff_value = coefficients(ind);
      // If pred_value is not in the pred_to_values map, the element will be
      // created. Otherwise a reference to the existing element is returned.
      pred_to_values[pred_value].push_back(coeff_value);
    }

    template <typename CoefficientType>
    void
    fill_map(
      const Table<1, CoefficientType> &coefficients,
      const std::function<
        std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
      std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
    {
      for (unsigned int i = 0; i < coefficients.size(0); i++)
        {
          const TableIndices<1> ind(i);
          fill_map_index(coefficients, ind, predicate, pred_to_values);
        }
    }

    template <typename CoefficientType>
    void
    fill_map(
      const Table<2, CoefficientType> &coefficients,
      const std::function<
        std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
      std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
    {
      for (unsigned int i = 0; i < coefficients.size(0); i++)
        for (unsigned int j = 0; j < coefficients.size(1); j++)
          {
            const TableIndices<2> ind(i, j);
            fill_map_index(coefficients, ind, predicate, pred_to_values);
          }
    }

    template <typename CoefficientType>
    void
    fill_map(
      const Table<3, CoefficientType> &coefficients,
      const std::function<
        std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
      std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
    {
      for (unsigned int i = 0; i < coefficients.size(0); i++)
        for (unsigned int j = 0; j < coefficients.size(1); j++)
          for (unsigned int k = 0; k < coefficients.size(2); k++)
            {
              const TableIndices<3> ind(i, j, k);
              fill_map_index(coefficients, ind, predicate, pred_to_values);
            }
    }


    template <typename Number>
    double
    complex_mean_value(const Number &value)
    {
      return value;
    }

    template <typename Number>
    double
    complex_mean_value(const std::complex<Number> &value)
    {
      AssertThrow(false,
                  ExcMessage(
                    "FESeries::process_coefficients() can not be used with"
                    "complex-valued coefficients and VectorTools::mean norm."));
      return std::abs(value);
    }
  } // namespace FESeriesImplementation
} // namespace internal


template <int dim, typename CoefficientType>
std::pair<std::vector<unsigned int>, std::vector<double>>
FESeries::process_coefficients(
  const Table<dim, CoefficientType> &coefficients,
  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
    &                         predicate,
  const VectorTools::NormType norm)
{
  std::vector<unsigned int> predicate_values;
  std::vector<double>       norm_values;

  // first, parse all table elements into a map of predicate values and
  // coefficients. We could have stored (predicate values ->TableIndicies) map,
  // but its processing would have been much harder later on.
  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
  internal::FESeriesImplementation::fill_map(coefficients,
                                             predicate,
                                             pred_to_values);

  // now go through the map and populate the @p norm_values based on @p norm:
  for (const auto &pred_to_value : pred_to_values)
    {
      predicate_values.push_back(pred_to_value.first);
      Vector<CoefficientType> values(pred_to_value.second.cbegin(),
                                     pred_to_value.second.cend());

      switch (norm)
        {
          case VectorTools::L2_norm:
            {
              norm_values.push_back(values.l2_norm());
              break;
            }
          case VectorTools::L1_norm:
            {
              norm_values.push_back(values.l1_norm());
              break;
            }
          case VectorTools::Linfty_norm:
            {
              norm_values.push_back(values.linfty_norm());
              break;
            }
          case VectorTools::mean:
            {
              norm_values.push_back(
                internal::FESeriesImplementation::complex_mean_value(
                  values.mean_value()));
              break;
            }
          default:
            AssertThrow(false, ExcNotImplemented());
            break;
        }
    }

  return std::make_pair(predicate_values, norm_values);
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif // dealii_fe_series_H