fe_data.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2023 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_data_h
#define dealii_fe_data_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
 
#include <deal.II/grid/reference_cell.h>
 
#include <deal.II/lac/block_indices.h>
 
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations:
#ifndef DOXYGEN
template <int dim>
class FiniteElementData;
#endif
 
namespace FiniteElementDomination
{
  enum Domination
  {
    this_element_dominates,
    other_element_dominates,
    neither_element_dominates,
    either_element_can_dominate,
    no_requirements
  };
 
 
  inline Domination
  operator&(const Domination d1, const Domination d2);
} // namespace FiniteElementDomination
 
namespace internal
{
  struct GenericDoFsPerObject
  {
    std::vector<std::vector<unsigned int>> dofs_per_object_exclusive;
 
    std::vector<std::vector<unsigned int>> dofs_per_object_inclusive;
 
    std::vector<std::vector<unsigned int>> object_index;
 
    std::vector<std::vector<unsigned int>> first_object_index_on_face;
 
    template <int dim>
    static GenericDoFsPerObject
    generate(const FiniteElementData<dim> &fe);
  };
} // namespace internal
 
template <int dim>
class FiniteElementData
{
public:
  enum Conformity
  {
    unknown = 0x00,
 
    L2 = 0x01,
 
    Hcurl = 0x02,
 
    Hdiv = 0x04,
 
    H1 = Hcurl | Hdiv,
 
    H2 = 0x0e
  };
 
  static constexpr unsigned int dimension = dim;
 
private:
  const ReferenceCell reference_cell_kind;
 
  const unsigned int number_of_unique_2d_subobjects;
 
  const unsigned int number_unique_faces;
 
public:
  const unsigned int dofs_per_vertex;
 
  const unsigned int dofs_per_line;
 
private:
  const std::vector<unsigned int> n_dofs_on_quad;
 
public:
  const unsigned int dofs_per_quad;
 
private:
  const unsigned int dofs_per_quad_max;
 
public:
  const unsigned int dofs_per_hex;
 
  const unsigned int first_line_index;
 
private:
  const std::vector<unsigned int> first_index_of_quads;
 
public:
  const unsigned int first_quad_index;
 
  const unsigned int first_hex_index;
 
private:
  const std::vector<unsigned int> first_line_index_of_faces;
 
public:
  const unsigned int first_face_line_index;
 
private:
  const std::vector<unsigned int> first_quad_index_of_faces;
 
public:
  const unsigned int first_face_quad_index;
 
private:
  const std::vector<unsigned int> n_dofs_on_face;
 
public:
  const unsigned int dofs_per_face;
 
private:
  const unsigned int dofs_per_face_max;
 
public:
  const unsigned int dofs_per_cell;
 
  const unsigned int components;
 
  const unsigned int degree;
 
  const Conformity conforming_space;
 
  const BlockIndices block_indices_data;
 
  FiniteElementData(const std::vector<unsigned int> &dofs_per_object,
                    const unsigned int               n_components,
                    const unsigned int               degree,
                    const Conformity                 conformity = unknown,
                    const BlockIndices &block_indices = BlockIndices());
 
  FiniteElementData(const std::vector<unsigned int> &dofs_per_object,
                    const ReferenceCell              reference_cell,
                    const unsigned int               n_components,
                    const unsigned int               degree,
                    const Conformity                 conformity = unknown,
                    const BlockIndices &block_indices = BlockIndices());
 
  FiniteElementData(const internal::GenericDoFsPerObject &data,
                    const ReferenceCell                   reference_cell,
                    const unsigned int                    n_components,
                    const unsigned int                    degree,
                    const Conformity                      conformity = unknown,
                    const BlockIndices &block_indices = BlockIndices());
 
  ReferenceCell
  reference_cell() const;
 
  unsigned int
  n_unique_2d_subobjects() const;
 
  unsigned int
  n_unique_faces() const;
 
  unsigned int
  n_dofs_per_vertex() const;
 
  unsigned int
  n_dofs_per_line() const;
 
  unsigned int
  n_dofs_per_quad(unsigned int face_no = 0) const;
 
  unsigned int
  max_dofs_per_quad() const;
 
  unsigned int
  n_dofs_per_hex() const;
 
  unsigned int
  n_dofs_per_face(unsigned int face_no = 0, unsigned int child = 0) const;
 
  unsigned int
  max_dofs_per_face() const;
 
  unsigned int
  n_dofs_per_cell() const;
 
  template <int structdim>
  unsigned int
  n_dofs_per_object(const unsigned int i = 0) const;
 
  unsigned int
  n_components() const;
 
  unsigned int
  n_blocks() const;
 
  const BlockIndices &
  block_indices() const;
 
  unsigned int
  tensor_degree() const;
 
  bool
  conforms(const Conformity) const;
 
  bool
  operator==(const FiniteElementData &) const;
 
  unsigned int
  get_first_line_index() const;
 
  unsigned int
  get_first_quad_index(const unsigned int quad_no = 0) const;
 
  unsigned int
  get_first_hex_index() const;
 
  unsigned int
  get_first_face_line_index(const unsigned int face_no = 0) const;
 
  unsigned int
  get_first_face_quad_index(const unsigned int face_no = 0) const;
};
 
namespace internal
{
  internal::GenericDoFsPerObject
  expand(const unsigned int               dim,
         const std::vector<unsigned int> &dofs_per_object,
         const ReferenceCell              reference_cell);
} // namespace internal
 
 
 
// --------- inline and template functions ---------------
 
 
#ifndef DOXYGEN
 
namespace FiniteElementDomination
{
  inline Domination
  operator&(const Domination d1, const Domination d2)
  {
    // go through the entire list of possibilities. note that if we were into
    // speed, obfuscation and cared enough, we could implement this operator
    // by doing a bitwise & (and) if we gave these values to the enum values:
    // neither_element_dominates=0, this_element_dominates=1,
    // other_element_dominates=2, either_element_can_dominate=3
    // =this_element_dominates|other_element_dominates
    switch (d1)
      {
        case this_element_dominates:
          if ((d2 == this_element_dominates) ||
              (d2 == either_element_can_dominate) || (d2 == no_requirements))
            return this_element_dominates;
          else
            return neither_element_dominates;
 
        case other_element_dominates:
          if ((d2 == other_element_dominates) ||
              (d2 == either_element_can_dominate) || (d2 == no_requirements))
            return other_element_dominates;
          else
            return neither_element_dominates;
 
        case neither_element_dominates:
          return neither_element_dominates;
 
        case either_element_can_dominate:
          if (d2 == no_requirements)
            return either_element_can_dominate;
          else
            return d2;
 
        case no_requirements:
          return d2;
 
        default:
          // shouldn't get here
          Assert(false, ExcInternalError());
      }
 
    return neither_element_dominates;
  }
} // namespace FiniteElementDomination
 
 
template <int dim>
inline ReferenceCell
FiniteElementData<dim>::reference_cell() const
{
  return reference_cell_kind;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_unique_2d_subobjects() const
{
  return number_of_unique_2d_subobjects;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_unique_faces() const
{
  return number_unique_faces;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_vertex() const
{
  return dofs_per_vertex;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_line() const
{
  return dofs_per_line;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_quad(unsigned int face_no) const
{
  return n_dofs_on_quad[n_dofs_on_quad.size() == 1 ? 0 : face_no];
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::max_dofs_per_quad() const
{
  return dofs_per_quad_max;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_hex() const
{
  return dofs_per_hex;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_face(unsigned int face_no,
                                        unsigned int child_no) const
{
  (void)child_no;
 
  return n_dofs_on_face[n_dofs_on_face.size() == 1 ? 0 : face_no];
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::max_dofs_per_face() const
{
  return dofs_per_face_max;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_cell() const
{
  return dofs_per_cell;
}
 
 
 
template <int dim>
template <int structdim>
inline unsigned int
FiniteElementData<dim>::n_dofs_per_object(const unsigned int i) const
{
  switch (structdim)
    {
      case 0:
        return n_dofs_per_vertex();
      case 1:
        return n_dofs_per_line();
      case 2:
        return n_dofs_per_quad((structdim == 2 && dim == 3) ? i : 0);
      case 3:
        return n_dofs_per_hex();
      default:
        Assert(false, ExcInternalError());
    }
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_components() const
{
  return components;
}
 
 
 
template <int dim>
inline const BlockIndices &
FiniteElementData<dim>::block_indices() const
{
  return block_indices_data;
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::n_blocks() const
{
  return block_indices_data.size();
}
 
 
 
template <int dim>
inline unsigned int
FiniteElementData<dim>::tensor_degree() const
{
  return degree;
}
 
 
template <int dim>
inline bool
FiniteElementData<dim>::conforms(const Conformity space) const
{
  return ((space & conforming_space) == space);
}
 
 
 
template <int dim>
unsigned int
FiniteElementData<dim>::get_first_line_index() const
{
  return first_line_index;
}
 
template <int dim>
unsigned int
FiniteElementData<dim>::get_first_quad_index(const unsigned int quad_no) const
{
  if (first_index_of_quads.size() == 1)
    return first_index_of_quads[0] + quad_no * n_dofs_per_quad(0);
  else
    return first_index_of_quads[quad_no];
}
 
template <int dim>
unsigned int
FiniteElementData<dim>::get_first_hex_index() const
{
  return first_hex_index;
}
 
template <int dim>
unsigned int
FiniteElementData<dim>::get_first_face_line_index(
  const unsigned int face_no) const
{
  return first_line_index_of_faces[first_line_index_of_faces.size() == 1 ?
                                     0 :
                                     face_no];
}
 
template <int dim>
unsigned int
FiniteElementData<dim>::get_first_face_quad_index(
  const unsigned int face_no) const
{
  return first_quad_index_of_faces[first_quad_index_of_faces.size() == 1 ?
                                     0 :
                                     face_no];
}
 
template <int dim>
internal::GenericDoFsPerObject
internal::GenericDoFsPerObject::generate(const FiniteElementData<dim> &fe)
{
  const auto reference_cell = fe.reference_cell();
 
  internal::GenericDoFsPerObject result;
 
  result.dofs_per_object_exclusive.resize(4);
  result.dofs_per_object_inclusive.resize(4);
  result.object_index.resize(4);
 
  unsigned int counter = 0;
 
  for (const unsigned int v : reference_cell.vertex_indices())
    {
      const auto c = fe.template n_dofs_per_object<0>(v);
 
      result.dofs_per_object_exclusive[0].emplace_back(c);
      result.dofs_per_object_inclusive[0].emplace_back(c);
      result.object_index[0].emplace_back(counter);
 
      counter += c;
    }
 
  if (dim >= 2)
    for (const unsigned int l : reference_cell.line_indices())
      {
        const auto c = fe.template n_dofs_per_object<1>(l);
 
        result.dofs_per_object_exclusive[1].emplace_back(c);
        result.dofs_per_object_inclusive[1].emplace_back(
          c + 2 * fe.template n_dofs_per_object<0>());
        result.object_index[1].emplace_back(counter);
 
        counter += c;
      }
 
  if (dim == 3)
    for (const unsigned int f : reference_cell.face_indices())
      {
        const auto c = fe.template n_dofs_per_object<2>(f);
 
        result.dofs_per_object_exclusive[2].emplace_back(c);
        result.dofs_per_object_inclusive[2].emplace_back(fe.n_dofs_per_face(f));
        result.object_index[2].emplace_back(counter);
 
        counter += c;
      }
 
  {
    const auto c = fe.template n_dofs_per_object<dim>();
 
    result.dofs_per_object_exclusive[dim].emplace_back(c);
    result.dofs_per_object_inclusive[dim].emplace_back(fe.n_dofs_per_cell());
    result.object_index[dim].emplace_back(counter);
 
    counter += c;
  }
 
  for (unsigned int d = dim + 1; d <= 3; ++d)
    {
      result.dofs_per_object_exclusive[d].emplace_back(0);
      result.dofs_per_object_inclusive[d].emplace_back(0);
      result.object_index[d].emplace_back(counter);
    }
 
  result.first_object_index_on_face.resize(3);
  for (const unsigned int face_no : reference_cell.face_indices())
    {
      result.first_object_index_on_face[0].emplace_back(0);
 
      result.first_object_index_on_face[1].emplace_back(
        fe.get_first_face_line_index(face_no));
 
      result.first_object_index_on_face[2].emplace_back(
        fe.get_first_face_quad_index(face_no));
    }
 
  return result;
}
 
 
#endif // DOXYGEN
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif