fe_poly_tensor.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2005 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_fe_poly_tensor_h
#define dealii_fe_poly_tensor_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/mutex.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/tensor_polynomials_base.h>
 
#include <deal.II/fe/fe.h>
 
#include <deal.II/lac/full_matrix.h>
 
#include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
template <int dim, int spacedim = dim>
class FE_PolyTensor : public FiniteElement<dim, spacedim>
{
public:
  FE_PolyTensor(const TensorPolynomialsBase<dim> &polynomials,
                const FiniteElementData<dim>     &fe_data,
                const std::vector<bool>          &restriction_is_additive_flags,
                const std::vector<ComponentMask> &nonzero_components);
 
 
  FE_PolyTensor(const FE_PolyTensor &fe);
 
  // for documentation, see the FiniteElement base class
  virtual UpdateFlags
  requires_update_flags(const UpdateFlags update_flags) const override;
 
  virtual double
  shape_value(const unsigned int i, const Point<dim> &p) const override;
 
  // documentation inherited from the base class
  virtual double
  shape_value_component(const unsigned int i,
                        const Point<dim>  &p,
                        const unsigned int component) const override;
 
  virtual Tensor<1, dim>
  shape_grad(const unsigned int i, const Point<dim> &p) const override;
 
  // documentation inherited from the base class
  virtual Tensor<1, dim>
  shape_grad_component(const unsigned int i,
                       const Point<dim>  &p,
                       const unsigned int component) const override;
 
  virtual Tensor<2, dim>
  shape_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
  // documentation inherited from the base class
  virtual Tensor<2, dim>
  shape_grad_grad_component(const unsigned int i,
                            const Point<dim>  &p,
                            const unsigned int component) const override;
 
protected:
#ifndef DOXYGEN
  class InternalData;
#endif
 
  void
  compute_fill(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const typename QProjector<dim>::DataSetDescriptor          &offset,
    const unsigned int                                          n_q_points,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                              &mapping_data,
    const typename FE_PolyTensor<dim, spacedim>::InternalData &fe_internal,
    internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
      &output_data) const;
 
  std::vector<MappingKind> mapping_kind;
 
  bool
  single_mapping_kind() const;
 
  bool
  adjust_quad_dof_sign_for_face_orientation(
    const unsigned int                 index,
    const unsigned int                 face_no,
    const types::geometric_orientation combined_orientation) const;
 
  std::vector<Table<2, bool>> adjust_quad_dof_sign_for_face_orientation_table;
 
  MappingKind
  get_mapping_kind(const unsigned int i) const;
 
  /* NOTE: The following function has its definition inlined into the class
     declaration because we otherwise run into a compiler error with MS Visual
     Studio. */
  virtual std::unique_ptr<
    typename FiniteElement<dim, spacedim>::InternalDataBase>
  get_data(
    const UpdateFlags             update_flags,
    const Mapping<dim, spacedim> &mapping,
    const Quadrature<dim>        &quadrature,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override
  {
    (void)mapping;
    (void)output_data;
    // generate a new data object and
    // initialize some fields
    std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
          data_ptr   = std::make_unique<InternalData>();
    auto &data       = dynamic_cast<InternalData &>(*data_ptr);
    data.update_each = requires_update_flags(update_flags);
 
    const unsigned int n_q_points = quadrature.size();
 
    // some scratch arrays
    std::vector<Tensor<1, dim>> values(0);
    std::vector<Tensor<2, dim>> grads(0);
    std::vector<Tensor<3, dim>> grad_grads(0);
    std::vector<Tensor<4, dim>> third_derivatives(0);
    std::vector<Tensor<5, dim>> fourth_derivatives(0);
 
    if (update_flags & (update_values | update_gradients | update_hessians))
      data.dof_sign_change.resize(this->dofs_per_cell);
 
    // initialize fields only if really
    // necessary. otherwise, don't
    // allocate memory
 
    const bool update_transformed_shape_values =
      std::any_of(this->mapping_kind.begin(),
                  this->mapping_kind.end(),
                  [](const MappingKind t) { return t != mapping_none; });
 
    const bool update_transformed_shape_grads =
      std::any_of(this->mapping_kind.begin(),
                  this->mapping_kind.end(),
                  [](const MappingKind t) {
                    return (t == mapping_raviart_thomas || t == mapping_piola ||
                            t == mapping_nedelec || t == mapping_contravariant);
                  });
 
    const bool update_transformed_shape_hessian_tensors =
      update_transformed_shape_values;
 
    if (update_flags & update_values)
      {
        values.resize(this->n_dofs_per_cell());
        data.shape_values.reinit(this->n_dofs_per_cell(), n_q_points);
        if (update_transformed_shape_values)
          data.transformed_shape_values.resize(n_q_points);
      }
 
    if (update_flags & update_gradients)
      {
        grads.resize(this->n_dofs_per_cell());
        data.shape_grads.reinit(this->n_dofs_per_cell(), n_q_points);
        data.transformed_shape_grads.resize(n_q_points);
 
        if (update_transformed_shape_grads)
          data.untransformed_shape_grads.resize(n_q_points);
      }
 
    if (update_flags & update_hessians)
      {
        grad_grads.resize(this->n_dofs_per_cell());
        data.shape_grad_grads.reinit(this->n_dofs_per_cell(), n_q_points);
        data.transformed_shape_hessians.resize(n_q_points);
        if (update_transformed_shape_hessian_tensors)
          data.untransformed_shape_hessian_tensors.resize(n_q_points);
      }
 
    // Compute shape function values
    // and derivatives and hessians on
    // the reference cell.
    // Make sure, that for the
    // node values N_i holds
    // N_i(v_j)=\delta_ij for all basis
    // functions v_j
    if (update_flags & (update_values | update_gradients))
      for (unsigned int k = 0; k < n_q_points; ++k)
        {
          poly_space->evaluate(quadrature.point(k),
                               values,
                               grads,
                               grad_grads,
                               third_derivatives,
                               fourth_derivatives);
 
          if (update_flags & update_values)
            {
              if (inverse_node_matrix.n_cols() == 0)
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  data.shape_values[i][k] = values[i];
              else
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  {
                    Tensor<1, dim> add_values;
                    for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
                      add_values += inverse_node_matrix(j, i) * values[j];
                    data.shape_values[i][k] = add_values;
                  }
            }
 
          if (update_flags & update_gradients)
            {
              if (inverse_node_matrix.n_cols() == 0)
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  {
                    if constexpr (dim == spacedim)
                      data.shape_grads[i][k] = grads[i];
                    else
                      for (unsigned int d = 0; d < dim; ++d)
                        data.shape_grads[i][k][d] = grads[i][d];
                  }
              else
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  {
                    Tensor<2, dim> add_grads;
                    for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
                      add_grads += inverse_node_matrix(j, i) * grads[j];
                    if constexpr (dim == spacedim)
                      data.shape_grads[i][k] = add_grads;
                    else
                      for (unsigned int d = 0; d < dim; ++d)
                        data.shape_grads[i][k][d] = add_grads[d];
                  }
            }
 
          if (update_flags & update_hessians)
            {
              if (inverse_node_matrix.n_cols() == 0)
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  {
                    if constexpr (dim == spacedim)
                      data.shape_grad_grads[i][k] = grad_grads[i];
                    else
                      for (unsigned int d = 0; d < dim; ++d)
                        for (unsigned int e = 0; e < dim; ++e)
                          data.shape_grad_grads[i][k][d][e] =
                            grad_grads[i][d][e];
                  }
              else
                for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
                  {
                    Tensor<3, dim> add_grad_grads;
                    for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
                      add_grad_grads +=
                        inverse_node_matrix(j, i) * grad_grads[j];
                    if constexpr (dim == spacedim)
                      data.shape_grad_grads[i][k] = add_grad_grads;
                    else
                      for (unsigned int d = 0; d < dim; ++d)
                        for (unsigned int e = 0; e < dim; ++e)
                          data.shape_grad_grads[i][k][d][e] =
                            add_grad_grads[d][e];
                  }
            }
        }
    return data_ptr;
  }
 
  virtual void
  fill_fe_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const CellSimilarity::Similarity                            cell_similarity,
    const Quadrature<dim>                                      &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  using FiniteElement<dim, spacedim>::fill_fe_face_values;
 
  virtual void
  fill_fe_face_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_no,
    const hp::QCollection<dim - 1>                             &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  virtual void
  fill_fe_subface_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_no,
    const unsigned int                                          sub_no,
    const Quadrature<dim - 1>                                  &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
  {
  public:
    Table<2, Tensor<1, dim>> shape_values;
 
    Table<2, DerivativeForm<1, dim, spacedim>> shape_grads;
 
    Table<2, DerivativeForm<2, dim, spacedim>> shape_grad_grads;
 
    mutable std::vector<double>              dof_sign_change;
    mutable std::vector<Tensor<1, spacedim>> transformed_shape_values;
    // for shape_gradient computations
    mutable std::vector<Tensor<2, spacedim>> transformed_shape_grads;
    mutable std::vector<Tensor<2, dim>>      untransformed_shape_grads;
    // for shape_hessian computations
    mutable std::vector<Tensor<3, spacedim>> transformed_shape_hessians;
    mutable std::vector<Tensor<3, dim>> untransformed_shape_hessian_tensors;
  };
 
 
 
  const std::unique_ptr<const TensorPolynomialsBase<dim>> poly_space;
 
  FullMatrix<double> inverse_node_matrix;
 
  mutable Threads::Mutex cache_mutex;
 
  mutable Point<dim> cached_point;
 
  mutable std::vector<Tensor<1, dim>> cached_values;
 
  mutable std::vector<Tensor<2, dim>> cached_grads;
 
  mutable std::vector<Tensor<3, dim>> cached_grad_grads;
};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif