fe_dgp_monomial.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 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/std_cxx14/memory.h>

#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_tools.h>

#include <sstream>


DEAL_II_NAMESPACE_OPEN

namespace internal
{
  namespace FE_DGPMonomial
  {
    namespace
    {
      // storage of hand-chosen support
      // points
      //
      // For dim=2, dofs_per_cell of
      // FE_DGPMonomial(k) is given by
      // 0.5(k+1)(k+2), i.e.
      //
      // k    0  1  2  3  4  5  6  7
      // dofs 1  3  6 10 15 21 28 36
      //
      // indirect access of unit points:
      // the points for degree k are
      // located at
      //
      // points[start_index[k]..start_index[k+1]-1]
      const unsigned int start_index2d[6] = {0, 1, 4, 10, 20, 35};
      const double       points2d[35][2]  = {
        {0, 0},       {0, 0},       {1, 0},       {0, 1},       {0, 0},
        {1, 0},       {0, 1},       {1, 1},       {0.5, 0},     {0, 0.5},
        {0, 0},       {1, 0},       {0, 1},       {1, 1},       {1. / 3., 0},
        {2. / 3., 0}, {0, 1. / 3.}, {0, 2. / 3.}, {0.5, 1},     {1, 0.5},
        {0, 0},       {1, 0},       {0, 1},       {1, 1},       {0.25, 0},
        {0.5, 0},     {0.75, 0},    {0, 0.25},    {0, 0.5},     {0, 0.75},
        {1. / 3., 1}, {2. / 3., 1}, {1, 1. / 3.}, {1, 2. / 3.}, {0.5, 0.5}};

      // For dim=3, dofs_per_cell of
      // FE_DGPMonomial(k) is given by
      // 1./6.(k+1)(k+2)(k+3), i.e.
      //
      // k    0  1  2  3  4  5  6   7
      // dofs 1  4 10 20 35 56 84 120
      const unsigned int start_index3d[6] = {0, 1, 5, 15 /*,35*/};
      const double       points3d[35][3]  = {{0, 0, 0},
                                      {0, 0, 0},
                                      {1, 0, 0},
                                      {0, 1, 0},
                                      {0, 0, 1},
                                      {0, 0, 0},
                                      {1, 0, 0},
                                      {0, 1, 0},
                                      {0, 0, 1},
                                      {0.5, 0, 0},
                                      {0, 0.5, 0},
                                      {0, 0, 0.5},
                                      {1, 1, 0},
                                      {1, 0, 1},
                                      {0, 1, 1}};


      template <int dim>
      void
      generate_unit_points(const unsigned int, std::vector<Point<dim>> &);

      template <>
      void
      generate_unit_points(const unsigned int k, std::vector<Point<1>> &p)
      {
        Assert(p.size() == k + 1, ExcDimensionMismatch(p.size(), k + 1));
        const double h = 1. / k;
        for (unsigned int i = 0; i < p.size(); ++i)
          p[i](0) = i * h;
      }

      template <>
      void
      generate_unit_points(const unsigned int k, std::vector<Point<2>> &p)
      {
        Assert(k <= 4, ExcNotImplemented());
        Assert(p.size() == start_index2d[k + 1] - start_index2d[k],
               ExcInternalError());
        for (unsigned int i = 0; i < p.size(); ++i)
          {
            p[i](0) = points2d[start_index2d[k] + i][0];
            p[i](1) = points2d[start_index2d[k] + i][1];
          }
      }

      template <>
      void
      generate_unit_points(const unsigned int k, std::vector<Point<3>> &p)
      {
        Assert(k <= 2, ExcNotImplemented());
        Assert(p.size() == start_index3d[k + 1] - start_index3d[k],
               ExcInternalError());
        for (unsigned int i = 0; i < p.size(); ++i)
          {
            p[i](0) = points3d[start_index3d[k] + i][0];
            p[i](1) = points3d[start_index3d[k] + i][1];
            p[i](2) = points3d[start_index3d[k] + i][2];
          }
      }
    } // namespace
  }   // namespace FE_DGPMonomial
} // namespace internal



template <int dim>
FE_DGPMonomial<dim>::FE_DGPMonomial(const unsigned int degree)
  : FE_Poly<PolynomialsP<dim>, dim>(
      PolynomialsP<dim>(degree),
      FiniteElementData<dim>(get_dpo_vector(degree),
                             1,
                             degree,
                             FiniteElementData<dim>::L2),
      std::vector<bool>(
        FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
        true),
      std::vector<ComponentMask>(
        FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
        std::vector<bool>(1, true)))
{
  Assert(this->poly_space.n() == this->dofs_per_cell, ExcInternalError());
  Assert(this->poly_space.degree() == this->degree, ExcInternalError());

  // DG doesn't have constraints, so
  // leave them empty

  // Reinit the vectors of
  // restriction and prolongation
  // matrices to the right sizes
  this->reinit_restriction_and_prolongation_matrices();
  // Fill prolongation matrices with embedding operators
  FETools::compute_embedding_matrices(*this, this->prolongation);
  // Fill restriction matrices with L2-projection
  FETools::compute_projection_matrices(*this, this->restriction);
}



template <int dim>
std::string
FE_DGPMonomial<dim>::get_name() const
{
  // note that the
  // FETools::get_fe_by_name
  // function depends on the
  // particular format of the string
  // this function returns, so they
  // have to be kept in synch

  std::ostringstream namebuf;
  namebuf << "FE_DGPMonomial<" << dim << ">(" << this->degree << ")";

  return namebuf.str();
}



template <int dim>
std::unique_ptr<FiniteElement<dim, dim>>
FE_DGPMonomial<dim>::clone() const
{
  return std_cxx14::make_unique<FE_DGPMonomial<dim>>(*this);
}



// TODO: Remove this function and use the one in FETools, if needed
template <int dim>
void
FE_DGPMonomial<dim>::get_interpolation_matrix(
  const FiniteElement<dim> &source_fe,
  FullMatrix<double> &      interpolation_matrix) const
{
  const FE_DGPMonomial<dim> *source_dgp_monomial =
    dynamic_cast<const FE_DGPMonomial<dim> *>(&source_fe);

  if (source_dgp_monomial)
    {
      // ok, source_fe is a DGP_Monomial
      // element. Then, the interpolation
      // matrix is simple
      const unsigned int m = interpolation_matrix.m();
      const unsigned int n = interpolation_matrix.n();
      (void)m;
      (void)n;
      Assert(m == this->dofs_per_cell,
             ExcDimensionMismatch(m, this->dofs_per_cell));
      Assert(n == source_dgp_monomial->dofs_per_cell,
             ExcDimensionMismatch(n, source_dgp_monomial->dofs_per_cell));

      const unsigned int min_mn =
        interpolation_matrix.m() < interpolation_matrix.n() ?
          interpolation_matrix.m() :
          interpolation_matrix.n();

      for (unsigned int i = 0; i < min_mn; ++i)
        interpolation_matrix(i, i) = 1.;
    }
  else
    {
      std::vector<Point<dim>> unit_points(this->dofs_per_cell);
      internal::FE_DGPMonomial::generate_unit_points(this->degree, unit_points);

      FullMatrix<double> source_fe_matrix(unit_points.size(),
                                          source_fe.dofs_per_cell);
      for (unsigned int j = 0; j < source_fe.dofs_per_cell; ++j)
        for (unsigned int k = 0; k < unit_points.size(); ++k)
          source_fe_matrix(k, j) = source_fe.shape_value(j, unit_points[k]);

      FullMatrix<double> this_matrix(this->dofs_per_cell, this->dofs_per_cell);
      for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
        for (unsigned int k = 0; k < unit_points.size(); ++k)
          this_matrix(k, j) = this->poly_space.compute_value(j, unit_points[k]);

      this_matrix.gauss_jordan();

      this_matrix.mmult(interpolation_matrix, source_fe_matrix);
    }
}



template <int dim>
void
FE_DGPMonomial<dim>::initialize_restriction()
{
  Assert(false, ExcNotImplemented());
}


//---------------------------------------------------------------------------
// Auxiliary functions
//---------------------------------------------------------------------------


template <int dim>
std::vector<unsigned int>
FE_DGPMonomial<dim>::get_dpo_vector(const unsigned int deg)
{
  std::vector<unsigned int> dpo(dim + 1, 0U);
  dpo[dim] = deg + 1;
  for (unsigned int i = 1; i < dim; ++i)
    {
      dpo[dim] *= deg + 1 + i;
      dpo[dim] /= i + 1;
    }
  return dpo;
}


template <int dim>
void
FE_DGPMonomial<dim>::get_face_interpolation_matrix(
  const FiniteElement<dim> &x_source_fe,
  FullMatrix<double> &      interpolation_matrix) const
{
  // this is only implemented, if the source
  // FE is also a DGPMonomial element. in that case,
  // both elements have no dofs on their
  // faces and the face interpolation matrix
  // is necessarily empty -- i.e. there isn't
  // much we need to do here.
  (void)interpolation_matrix;
  AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
                (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
                 nullptr),
              typename FiniteElement<dim>::ExcInterpolationNotImplemented());

  Assert(interpolation_matrix.m() == 0,
         ExcDimensionMismatch(interpolation_matrix.m(), 0));
  Assert(interpolation_matrix.n() == 0,
         ExcDimensionMismatch(interpolation_matrix.n(), 0));
}



template <int dim>
void
FE_DGPMonomial<dim>::get_subface_interpolation_matrix(
  const FiniteElement<dim> &x_source_fe,
  const unsigned int,
  FullMatrix<double> &interpolation_matrix) const
{
  // this is only implemented, if the source
  // FE is also a DGPMonomial element. in that case,
  // both elements have no dofs on their
  // faces and the face interpolation matrix
  // is necessarily empty -- i.e. there isn't
  // much we need to do here.
  (void)interpolation_matrix;
  AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
                (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
                 nullptr),
              typename FiniteElement<dim>::ExcInterpolationNotImplemented());

  Assert(interpolation_matrix.m() == 0,
         ExcDimensionMismatch(interpolation_matrix.m(), 0));
  Assert(interpolation_matrix.n() == 0,
         ExcDimensionMismatch(interpolation_matrix.n(), 0));
}



template <int dim>
bool
FE_DGPMonomial<dim>::hp_constraints_are_implemented() const
{
  return true;
}



template <int dim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGPMonomial<dim>::hp_vertex_dof_identities(
  const FiniteElement<dim> &fe_other) const
{
  // there are no such constraints for DGPMonomial
  // elements at all
  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
    return std::vector<std::pair<unsigned int, unsigned int>>();
  else
    {
      Assert(false, ExcNotImplemented());
      return std::vector<std::pair<unsigned int, unsigned int>>();
    }
}



template <int dim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGPMonomial<dim>::hp_line_dof_identities(
  const FiniteElement<dim> &fe_other) const
{
  // there are no such constraints for DGPMonomial
  // elements at all
  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
    return std::vector<std::pair<unsigned int, unsigned int>>();
  else
    {
      Assert(false, ExcNotImplemented());
      return std::vector<std::pair<unsigned int, unsigned int>>();
    }
}



template <int dim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGPMonomial<dim>::hp_quad_dof_identities(
  const FiniteElement<dim> &fe_other) const
{
  // there are no such constraints for DGPMonomial
  // elements at all
  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
    return std::vector<std::pair<unsigned int, unsigned int>>();
  else
    {
      Assert(false, ExcNotImplemented());
      return std::vector<std::pair<unsigned int, unsigned int>>();
    }
}



template <int dim>
FiniteElementDomination::Domination
FE_DGPMonomial<dim>::compare_for_domination(const FiniteElement<dim> &fe_other,
                                            const unsigned int codim) const
{
  Assert(codim <= dim, ExcImpossibleInDim(dim));

  // vertex/line/face domination
  // ---------------------------
  if (codim > 0)
    // this is a discontinuous element, so by definition there will
    // be no constraints wherever this element comes together with
    // any other kind of element
    return FiniteElementDomination::no_requirements;

  // cell domination
  // ---------------
  if (const FE_DGPMonomial<dim> *fe_monomial_other =
        dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other))
    {
      if (this->degree < fe_monomial_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_monomial_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_Nothing<dim> *fe_nothing =
             dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
    {
      if (fe_nothing->is_dominating())
        return FiniteElementDomination::other_element_dominates;
      else
        // the FE_Nothing has no degrees of freedom and it is typically used
        // in a context where we don't require any continuity along the
        // interface
        return FiniteElementDomination::no_requirements;
    }

  Assert(false, ExcNotImplemented());
  return FiniteElementDomination::neither_element_dominates;
}



template <>
bool
FE_DGPMonomial<1>::has_support_on_face(const unsigned int,
                                       const unsigned int face_index) const
{
  return face_index == 1 || (face_index == 0 && this->degree == 0);
}



template <>
bool
FE_DGPMonomial<2>::has_support_on_face(const unsigned int shape_index,
                                       const unsigned int face_index) const
{
  bool support_on_face = false;
  if (face_index == 1 || face_index == 2)
    support_on_face = true;
  else
    {
      const std::array<unsigned int, 2> degrees =
        this->poly_space.directional_degrees(shape_index);

      if ((face_index == 0 && degrees[1] == 0) ||
          (face_index == 3 && degrees[0] == 0))
        support_on_face = true;
    }
  return support_on_face;
}



template <>
bool
FE_DGPMonomial<3>::has_support_on_face(const unsigned int shape_index,
                                       const unsigned int face_index) const
{
  bool support_on_face = false;
  if (face_index == 1 || face_index == 3 || face_index == 4)
    support_on_face = true;
  else
    {
      const std::array<unsigned int, 3> degrees =
        this->poly_space.directional_degrees(shape_index);

      if ((face_index == 0 && degrees[1] == 0) ||
          (face_index == 2 && degrees[2] == 0) ||
          (face_index == 5 && degrees[0] == 0))
        support_on_face = true;
    }
  return support_on_face;
}



template <int dim>
std::size_t
FE_DGPMonomial<dim>::memory_consumption() const
{
  Assert(false, ExcNotImplemented());
  return 0;
}



// explicit instantiations
#include "fe_dgp_monomial.inst"


DEAL_II_NAMESPACE_CLOSE