fe_interface_values.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2018 - 2020 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_interface_values_h
#define dealii_fe_interface_values_h

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

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

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim = dim>
class FEInterfaceValues
{
public:
  const unsigned int n_quadrature_points;


  FEInterfaceValues(const Mapping<dim, spacedim> &      mapping,
                    const FiniteElement<dim, spacedim> &fe,
                    const Quadrature<dim - 1> &         quadrature,
                    const UpdateFlags                   update_flags);

  FEInterfaceValues(const FiniteElement<dim, spacedim> &fe,
                    const Quadrature<dim - 1> &         quadrature,
                    const UpdateFlags                   update_flags);

  template <class CellIteratorType>
  void
  reinit(const CellIteratorType &                         cell,
         const unsigned int                               face_no,
         const unsigned int                               sub_face_no,
         const typename identity<CellIteratorType>::type &cell_neighbor,
         const unsigned int                               face_no_neighbor,
         const unsigned int                               sub_face_no_neighbor);

  template <class CellIteratorType>
  void
  reinit(const CellIteratorType &cell, const unsigned int face_no);

  const FEFaceValuesBase<dim, spacedim> &
  get_fe_face_values(const unsigned int cell_index) const;

  const Quadrature<dim - 1> &
  get_quadrature() const;

  UpdateFlags
  get_update_flags() const;

  bool
  at_boundary() const;

  double
  JxW(const unsigned int quadrature_point) const;

  const std::vector<double> &
  get_JxW_values() const;

  const std::vector<Tensor<1, spacedim>> &
  get_normal_vectors() const;

  const std::vector<Point<spacedim>> &
  get_quadrature_points() const;


  unsigned
  n_current_interface_dofs() const;

  std::vector<types::global_dof_index>
  get_interface_dof_indices() const;

  std::array<unsigned int, 2>
  interface_dof_to_dof_indices(const unsigned int interface_dof_index) const;

  Tensor<1, spacedim>
  normal(const unsigned int q_point_index) const;

  double
  shape_value(const bool         here_or_there,
              const unsigned int interface_dof_index,
              const unsigned int q_point,
              const unsigned int component = 0) const;

  double
  jump(const unsigned int interface_dof_index,
       const unsigned int q_point,
       const unsigned int component = 0) const;

  double
  average(const unsigned int interface_dof_index,
          const unsigned int q_point,
          const unsigned int component = 0) const;

  Tensor<1, dim>
  average_gradient(const unsigned int interface_dof_index,
                   const unsigned int q_point,
                   const unsigned int component = 0) const;

  Tensor<2, dim>
  average_hessian(const unsigned int interface_dof_index,
                  const unsigned int q_point,
                  const unsigned int component = 0) const;

  Tensor<1, dim>
  jump_gradient(const unsigned int interface_dof_index,
                const unsigned int q_point,
                const unsigned int component = 0) const;

  Tensor<2, dim>
  jump_hessian(const unsigned int interface_dof_index,
               const unsigned int q_point,
               const unsigned int component = 0) const;

  Tensor<3, dim>
  jump_3rd_derivative(const unsigned int interface_dof_index,
                      const unsigned int q_point,
                      const unsigned int component = 0) const;

private:
  std::vector<types::global_dof_index> interface_dof_indices;

  std::vector<std::array<unsigned int, 2>> dofmap;

  FEFaceValues<dim> internal_fe_face_values;

  FESubfaceValues<dim> internal_fe_subface_values;

  FEFaceValues<dim> internal_fe_face_values_neighbor;

  FESubfaceValues<dim> internal_fe_subface_values_neighbor;

  FEFaceValuesBase<dim> *fe_face_values;

  FEFaceValuesBase<dim> *fe_face_values_neighbor;
};



#ifndef DOXYGEN

/*---------------------- Inline functions ---------------------*/

template <int dim, int spacedim>
FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
  const Mapping<dim, spacedim> &      mapping,
  const FiniteElement<dim, spacedim> &fe,
  const Quadrature<dim - 1> &         quadrature,
  const UpdateFlags                   update_flags)
  : n_quadrature_points(quadrature.size())
  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
  , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
  , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
  , fe_face_values(nullptr)
  , fe_face_values_neighbor(nullptr)
{}



template <int dim, int spacedim>
FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
  const FiniteElement<dim, spacedim> &fe,
  const Quadrature<dim - 1> &         quadrature,
  const UpdateFlags                   update_flags)
  : n_quadrature_points(quadrature.size())
  , internal_fe_face_values(StaticMappingQ1<dim, spacedim>::mapping,
                            fe,
                            quadrature,
                            update_flags)
  , internal_fe_subface_values(StaticMappingQ1<dim, spacedim>::mapping,
                               fe,
                               quadrature,
                               update_flags)
  , internal_fe_face_values_neighbor(StaticMappingQ1<dim, spacedim>::mapping,
                                     fe,
                                     quadrature,
                                     update_flags)
  , internal_fe_subface_values_neighbor(StaticMappingQ1<dim, spacedim>::mapping,
                                        fe,
                                        quadrature,
                                        update_flags)
  , fe_face_values(nullptr)
  , fe_face_values_neighbor(nullptr)
{}



template <int dim, int spacedim>
template <class CellIteratorType>
void
FEInterfaceValues<dim, spacedim>::reinit(
  const CellIteratorType &                         cell,
  const unsigned int                               face_no,
  const unsigned int                               sub_face_no,
  const typename identity<CellIteratorType>::type &cell_neighbor,
  const unsigned int                               face_no_neighbor,
  const unsigned int                               sub_face_no_neighbor)
{
  if (sub_face_no == numbers::invalid_unsigned_int)
    {
      internal_fe_face_values.reinit(cell, face_no);
      fe_face_values = &internal_fe_face_values;
    }
  else
    {
      internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
      fe_face_values = &internal_fe_subface_values;
    }
  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
    {
      internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
      fe_face_values_neighbor = &internal_fe_face_values_neighbor;
    }
  else
    {
      internal_fe_subface_values_neighbor.reinit(cell_neighbor,
                                                 face_no_neighbor,
                                                 sub_face_no_neighbor);
      fe_face_values_neighbor = &internal_fe_subface_values_neighbor;
    }

  // Set up dof mapping and remove duplicates (for continuous elements).
  {
    // Get dof indices first:
    std::vector<types::global_dof_index> v(
      fe_face_values->get_fe().n_dofs_per_cell());
    cell->get_active_or_mg_dof_indices(v);
    std::vector<types::global_dof_index> v2(
      fe_face_values_neighbor->get_fe().n_dofs_per_cell());
    cell_neighbor->get_active_or_mg_dof_indices(v2);



    // Fill a map from the global dof index to the left and right
    // local index.
    std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
                                          tempmap;
    std::pair<unsigned int, unsigned int> invalid_entry(
      numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);

    for (unsigned int i = 0; i < v.size(); ++i)
      {
        // If not already existing, add an invalid entry:
        auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
        result.first->second.first = i;
      }

    for (unsigned int i = 0; i < v2.size(); ++i)
      {
        // If not already existing, add an invalid entry:
        auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
        result.first->second.second = i;
      }

    // Transfer from the map to the sorted std::vectors.
    dofmap.resize(tempmap.size());
    interface_dof_indices.resize(tempmap.size());
    unsigned int idx = 0;
    for (auto &x : tempmap)
      {
        interface_dof_indices[idx] = x.first;
        dofmap[idx]                = {{x.second.first, x.second.second}};
        ++idx;
      }
  }
}



template <int dim, int spacedim>
template <class CellIteratorType>
void
FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
                                         const unsigned int      face_no)
{
  internal_fe_face_values.reinit(cell, face_no);
  fe_face_values          = &internal_fe_face_values;
  fe_face_values_neighbor = nullptr;

  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
  cell->get_active_or_mg_dof_indices(interface_dof_indices);


  dofmap.resize(interface_dof_indices.size());

  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
    {
      dofmap[i] = {{i, numbers::invalid_unsigned_int}};
    }
}



template <int dim, int spacedim>
inline double
FEInterfaceValues<dim, spacedim>::JxW(const unsigned int q) const
{
  Assert(fe_face_values != nullptr,
         ExcMessage("This call requires a call to reinit() first."));
  return fe_face_values->JxW(q);
}



template <int dim, int spacedim>
const std::vector<double> &
FEInterfaceValues<dim, spacedim>::get_JxW_values() const
{
  Assert(fe_face_values != nullptr,
         ExcMessage("This call requires a call to reinit() first."));
  return fe_face_values->get_JxW_values();
}



template <int dim, int spacedim>
const std::vector<Tensor<1, spacedim>> &
FEInterfaceValues<dim, spacedim>::get_normal_vectors() const
{
  Assert(fe_face_values != nullptr,
         ExcMessage("This call requires a call to reinit() first."));
  return fe_face_values->get_normal_vectors();
}



template <int dim, int spacedim>
const Quadrature<dim - 1> &
FEInterfaceValues<dim, spacedim>::get_quadrature() const
{
  return internal_fe_face_values.get_quadrature();
}



template <int dim, int spacedim>
const std::vector<Point<spacedim>> &
FEInterfaceValues<dim, spacedim>::get_quadrature_points() const
{
  Assert(fe_face_values != nullptr,
         ExcMessage("This call requires a call to reinit() first."));
  return fe_face_values->get_quadrature_points();
}



template <int dim, int spacedim>
UpdateFlags
FEInterfaceValues<dim, spacedim>::get_update_flags() const
{
  return internal_fe_face_values.get_update_flags();
}



template <int dim, int spacedim>
unsigned
FEInterfaceValues<dim, spacedim>::n_current_interface_dofs() const
{
  Assert(
    interface_dof_indices.size() > 0,
    ExcMessage(
      "n_current_interface_dofs() is only available after a call to reinit()."));
  return interface_dof_indices.size();
}



template <int dim, int spacedim>
bool
FEInterfaceValues<dim, spacedim>::at_boundary() const
{
  return fe_face_values_neighbor == nullptr;
}



template <int dim, int spacedim>
std::vector<types::global_dof_index>
FEInterfaceValues<dim, spacedim>::get_interface_dof_indices() const
{
  return interface_dof_indices;
}



template <int dim, int spacedim>
std::array<unsigned int, 2>
FEInterfaceValues<dim, spacedim>::interface_dof_to_dof_indices(
  const unsigned int interface_dof_index) const
{
  AssertIndexRange(interface_dof_index, n_current_interface_dofs());
  return dofmap[interface_dof_index];
}



template <int dim, int spacedim>
const FEFaceValuesBase<dim, spacedim> &
FEInterfaceValues<dim, spacedim>::get_fe_face_values(
  const unsigned int cell_index) const
{
  AssertIndexRange(cell_index, 2);
  Assert(
    cell_index == 0 || !at_boundary(),
    ExcMessage(
      "You are on a boundary, so you can only ask for the first FEFaceValues object."));

  return (cell_index == 0) ? *fe_face_values : *fe_face_values_neighbor;
}



template <int dim, int spacedim>
Tensor<1, spacedim>
FEInterfaceValues<dim, spacedim>::normal(const unsigned int q_point_index) const
{
  return fe_face_values->normal_vector(q_point_index);
}



template <int dim, int spacedim>
double
FEInterfaceValues<dim, spacedim>::shape_value(
  const bool         here_or_there,
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (here_or_there && dof_pair[0] != numbers::invalid_unsigned_int)
    return get_fe_face_values(0).shape_value_component(dof_pair[0],
                                                       q_point,
                                                       component);
  if (!here_or_there && dof_pair[1] != numbers::invalid_unsigned_int)
    return get_fe_face_values(1).shape_value_component(dof_pair[1],
                                                       q_point,
                                                       component);

  return 0.0;
}



template <int dim, int spacedim>
double
FEInterfaceValues<dim, spacedim>::jump(const unsigned int interface_dof_index,
                                       const unsigned int q_point,
                                       const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  double value = 0.0;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += get_fe_face_values(0).shape_value_component(dof_pair[0],
                                                         q_point,
                                                         component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value -= get_fe_face_values(1).shape_value_component(dof_pair[1],
                                                         q_point,
                                                         component);
  return value;
}



template <int dim, int spacedim>
double
FEInterfaceValues<dim, spacedim>::average(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return 1.0 * get_fe_face_values(0).shape_value_component(dof_pair[0],
                                                             q_point,
                                                             component);

  double value = 0.0;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(0).shape_value_component(dof_pair[0],
                                                               q_point,
                                                               component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(1).shape_value_component(dof_pair[1],
                                                               q_point,
                                                               component);

  return value;
}



template <int dim, int spacedim>
Tensor<1, dim>
FEInterfaceValues<dim, spacedim>::average_gradient(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return get_fe_face_values(0).shape_grad_component(dof_pair[0],
                                                      q_point,
                                                      component);

  Tensor<1, dim> value;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
                                                              q_point,
                                                              component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
                                                              q_point,
                                                              component);

  return value;
}



template <int dim, int spacedim>
Tensor<2, dim>
FEInterfaceValues<dim, spacedim>::average_hessian(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
                                                         q_point,
                                                         component);

  Tensor<2, dim> value;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
                                                                 q_point,
                                                                 component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value += 0.5 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
                                                                 q_point,
                                                                 component);

  return value;
}



template <int dim, int spacedim>
Tensor<1, dim>
FEInterfaceValues<dim, spacedim>::jump_gradient(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return get_fe_face_values(0).shape_grad_component(dof_pair[0],
                                                      q_point,
                                                      component);

  Tensor<1, dim> value;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
                                                              q_point,
                                                              component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value += -1.0 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
                                                               q_point,
                                                               component);

  return value;
}



template <int dim, int spacedim>
Tensor<2, dim>
FEInterfaceValues<dim, spacedim>::jump_hessian(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
                                                         q_point,
                                                         component);

  Tensor<2, dim> value;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
                                                                 q_point,
                                                                 component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
                                                                  q_point,
                                                                  component);

  return value;
}


template <int dim, int spacedim>
Tensor<3, dim>
FEInterfaceValues<dim, spacedim>::jump_3rd_derivative(
  const unsigned int interface_dof_index,
  const unsigned int q_point,
  const unsigned int component) const
{
  const auto dof_pair = dofmap[interface_dof_index];

  if (at_boundary())
    return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
                                                                q_point,
                                                                component);

  Tensor<3, dim> value;

  if (dof_pair[0] != numbers::invalid_unsigned_int)
    value +=
      1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
                                                                 q_point,
                                                                 component);
  if (dof_pair[1] != numbers::invalid_unsigned_int)
    value +=
      -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1],
                                                                  q_point,
                                                                  component);

  return value;
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif