fe_evaluation_cell.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 by the hyper.deal authors
//
// This file is part of the hyper.deal library.
//
// The hyper.deal 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 3.0 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 hyper.deal.
//
// ---------------------------------------------------------------------
 
#ifndef HYPERDEAL_NDIM_FEEVALUATION_CELL
#define HYPERDEAL_NDIM_FEEVALUATION_CELL
 
#include <hyper.deal/base/config.h>
 
#include <deal.II/matrix_free/fe_evaluation.h>
 
#include <hyper.deal/matrix_free/fe_evaluation_base.h>
#include <hyper.deal/matrix_free/read_write_operation.h>
#include <hyper.deal/matrix_free/vector_access_internal.h>
 
namespace hyperdeal
{
  template <int dim_x,
            int dim_v,
            int degree,
            int n_points,
            typename Number,
            typename VectorizedArrayType>
  class FEEvaluation : public FEEvaluationBase<dim_x,
                                               dim_v,
                                               degree,
                                               n_points,
                                               Number,
                                               VectorizedArrayType>
  {
  public:
    // clang-format off
    static const unsigned int dim                    = dim_x + dim_v;
    static const unsigned int static_dofs_per_cell_x = dealii::Utilities::pow(degree + 1, dim_x);
    static const unsigned int static_dofs_per_cell_v = dealii::Utilities::pow(degree + 1, dim_v);
    static const unsigned int static_dofs_per_cell   = static_dofs_per_cell_x * static_dofs_per_cell_v;
    static const unsigned int n_q_points_x           = dealii::Utilities::pow(n_points, dim_x);
    static const unsigned int n_q_points_v           = dealii::Utilities::pow(n_points, dim_v);
    static const unsigned int n_q_points             = n_q_points_x * n_q_points_v;
    // clang-format on
 
    using PARENT = FEEvaluationBase<dim_x,
                                    dim_v,
                                    degree,
                                    n_points,
                                    Number,
                                    VectorizedArrayType>;
 
    static const unsigned int n_vectors   = PARENT::n_vectors;
    static const unsigned int n_vectors_v = PARENT::n_vectors_v;
 
    FEEvaluation(
      const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &matrix_free,
      const unsigned int                                           dof_no_x,
      const unsigned int                                           dof_no_v,
      const unsigned int                                           quad_no_x,
      const unsigned int                                           quad_no_v)
      : PARENT(matrix_free)
      , phi_x(this->matrix_free_x, dof_no_x, quad_no_x)
      , phi_v(this->matrix_free_v, dof_no_v, quad_no_v)
    {
      this->shape_values = &phi_x.get_shape_info().data[0].shape_values_eo;
      this->shape_gradients =
        &phi_x.get_shape_info().data[0].shape_gradients_collocation_eo;
    }
 
 
    void
    reinit(typename MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::ID
             cell_index) override
    {
      PARENT::reinit(cell_index);
      phi_x.reinit(this->macro_cell_x);
      phi_v.reinit(this->macro_cell_v);
    }
 
 
    void
    read_dof_values(
      const dealii::LinearAlgebra::distributed::Vector<Number> &src,
      VectorizedArrayType *                                     data) const
    {
      internal::MatrixFreeFunctions::ReadWriteOperation<Number>(
        this->matrix_free.get_dof_info(),
        this->matrix_free.get_face_info(),
        this->matrix_free.get_shape_info())
        .template process_cell<dim, degree>(
          internal::MatrixFreeFunctions::VectorReader<Number,
                                                      VectorizedArrayType>(),
          src.shared_vector_data(),
          data,
          this->macro);
    }
 
 
 
    void
    read_dof_values(
      const dealii::LinearAlgebra::distributed::Vector<Number> &src)
    {
      read_dof_values(src, &this->data[0]);
    }
 
 
 
    void
    set_dof_values(
      dealii::LinearAlgebra::distributed::Vector<Number> &dst) const
    {
      internal::MatrixFreeFunctions::ReadWriteOperation<Number>(
        this->matrix_free.get_dof_info(),
        this->matrix_free.get_face_info(),
        this->matrix_free.get_shape_info())
        .template process_cell<dim, degree>(
          internal::MatrixFreeFunctions::VectorSetter<Number,
                                                      VectorizedArrayType>(),
          dst.shared_vector_data(),
          &this->data[0],
          this->macro);
    }
 
 
 
    inline unsigned int
    n_vectorization_lanes_filled() const
    {
      return this->matrix_free_x.n_active_entries_per_cell_batch(
        this->macro_cell_x);
    }
 
 
    inline VectorizedArrayType
    submit_inplace(const VectorizedArrayType data, const unsigned int q) const
    {
      return submit_inplace(data, q, q % n_q_points_x, q / n_q_points_x);
    }
 
 
 
    inline VectorizedArrayType
    submit_inplace(const VectorizedArrayType data,
                   unsigned int /*q*/,
                   unsigned int qx,
                   unsigned int qv) const
    {
      return data * phi_x.JxW(qx) *
             phi_v.JxW(qv)[n_vectors_v == 1 ? 0 : this->lane_y];
    }
 
 
    inline DEAL_II_ALWAYS_INLINE //
      VectorizedArrayType
      JxW(unsigned int qx, unsigned int qv) const
    {
      return phi_x.JxW(qx) * phi_v.JxW(qv)[n_vectors_v == 1 ? 0 : this->lane_y];
    }
 
 
 
    inline DEAL_II_ALWAYS_INLINE //
      dealii::Tensor<1, dim_x, VectorizedArrayType>
      get_gradient_x(const VectorizedArrayType *__restrict grad_in,
                     const unsigned int q,
                     const unsigned int qx,
                     const unsigned int qv) const
    {
      (void)qv;
 
      dealii::Tensor<1, dim_x, VectorizedArrayType> result;
 
      const auto jacobian = phi_x.inverse_jacobian(qx);
 
      for (auto d = 0u; d < dim_x; d++)
        {
          result[d] = jacobian[d][0] * grad_in[q];
          for (auto e = 1u; e < dim_x; ++e)
            result[d] +=
              (jacobian[d][e] * grad_in[q + e * n_q_points_x * n_q_points_v]);
        }
 
      return result;
    }
 
 
 
    inline DEAL_II_ALWAYS_INLINE //
      dealii::Tensor<1, dim_v, VectorizedArrayType>
      get_gradient_v(const VectorizedArrayType *__restrict grad_in,
                     const unsigned int q,
                     const unsigned int qx,
                     const unsigned int qv) const
    {
      (void)qx;
 
      dealii::Tensor<1, dim_v, VectorizedArrayType> result;
 
      const auto jacobian = phi_v.inverse_jacobian(qv);
 
      for (auto d = 0u; d < dim_v; d++)
        {
          result[d] =
            jacobian[d][0][n_vectors_v == 1 ? 0 : this->lane_y] * grad_in[q];
          for (auto e = 1u; e < dim_v; ++e)
            result[d] += (jacobian[d][e][n_vectors_v == 1 ? 0 : this->lane_y] *
                          grad_in[q + e * n_q_points_x * n_q_points_v]);
        }
 
      return result;
    }
 
 
 
    template <bool do_add>
    inline DEAL_II_ALWAYS_INLINE //
      void
      submit_value(VectorizedArrayType *__restrict data_ptr_out,
                   const VectorizedArrayType values_in,
                   const unsigned int        q,
                   const unsigned int        qx,
                   const unsigned int        qv) const
    {
      const auto jxw =
        phi_x.JxW(qx) * phi_v.JxW(qv)[n_vectors_v == 1 ? 0 : this->lane_y];
 
      if (do_add)
        data_ptr_out[q] += values_in * jxw;
      else
        data_ptr_out[q] = values_in * jxw;
    }
 
 
 
    inline DEAL_II_ALWAYS_INLINE //
      void
      submit_gradient_x(VectorizedArrayType *__restrict data_ptr_out,
                        const VectorizedArrayType *__restrict grad_in,
                        const unsigned int q,
                        const unsigned int qx,
                        const unsigned int qv) const
    {
      const auto jxw =
        phi_x.JxW(qx) * phi_v.JxW(qv)[n_vectors_v == 1 ? 0 : this->lane_y];
      const auto jacobian = phi_x.inverse_jacobian(qx);
 
      for (auto d = 0u; d < dim_x; d++)
        {
          auto new_val = jacobian[0][d] * grad_in[0];
          for (auto e = 1u; e < dim_x; ++e)
            new_val += (jacobian[e][d] * grad_in[e]);
          data_ptr_out[q + d * n_q_points_x * n_q_points_v] = new_val * jxw;
        }
    }
 
 
 
    inline DEAL_II_ALWAYS_INLINE //
      void
      submit_gradient_v(VectorizedArrayType *__restrict data_ptr_out,
                        const VectorizedArrayType *__restrict grad_in,
                        const unsigned int q,
                        const unsigned int qx,
                        const unsigned int qv) const
    {
      const auto jxw =
        phi_x.JxW(qx) * phi_v.JxW(qv)[n_vectors_v == 1 ? 0 : this->lane_y];
      const auto jacobian = phi_v.inverse_jacobian(qv);
 
      for (auto d = 0u; d < dim_v; d++)
        {
          auto new_val =
            jacobian[0][d][n_vectors_v == 1 ? 0 : this->lane_y] * grad_in[0];
          for (auto e = 1u; e < dim_v; ++e)
            new_val += (jacobian[e][d][n_vectors_v == 1 ? 0 : this->lane_y] *
                        grad_in[e]);
          data_ptr_out[q + d * n_q_points_x * n_q_points_v] = new_val * jxw;
        }
    }
 
 
 
    inline dealii::Point<dim, VectorizedArrayType>
    get_quadrature_point(const unsigned int q) const
    {
      dealii::Point<dim, VectorizedArrayType> temp;
 
      auto t1 = this->phi_x.quadrature_point(q % n_q_points_x);
      for (int i = 0; i < dim_x; i++)
        temp[i] = t1[i];
      auto t2 = this->phi_v.quadrature_point(q / n_q_points_x);
      for (int i = 0; i < dim_v; i++)
        temp[i + dim_x] = t2[i][n_vectors_v == 1 ? 0 : this->lane_y];
 
      return temp;
    }
 
    inline dealii::Point<dim, VectorizedArrayType>
    get_quadrature_point(const unsigned int qx, const unsigned int qv) const
    {
      dealii::Point<dim, VectorizedArrayType> temp;
 
      auto t1 = this->phi_x.quadrature_point(qx);
      for (int i = 0; i < dim_x; i++)
        temp[i] = t1[i];
      auto t2 = this->phi_v.quadrature_point(qv);
      for (int i = 0; i < dim_v; i++)
        temp[i + dim_x] = t2[i][n_vectors_v == 1 ? 0 : this->lane_y];
 
      return temp;
    }
 
    inline dealii::Point<dim_v, VectorizedArrayType>
    get_quadrature_point_v(const unsigned int qv) const
    {
      dealii::Point<dim_v, VectorizedArrayType> temp;
 
      auto t2 = this->phi_v.quadrature_point(qv);
      for (int i = 0; i < dim_v; i++)
        temp[i] = t2[i][n_vectors_v == 1 ? 0 : this->lane_y];
 
      return temp;
    }
 
  protected:
    // clang-format off
    dealii::FEEvaluation<dim_x, degree, n_points, 1, Number, typename PARENT::VectorizedArrayTypeX> phi_x;
    dealii::FEEvaluation<dim_v, degree, n_points, 1, Number, typename PARENT::VectorizedArrayTypeV> phi_v;
    // clang-format on
  };
 
} // namespace hyperdeal
 
#endif