fe_poly.h

// ---------------------------------------------------------------------
//
// 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_fe_poly_h
#define dealii_fe_poly_h


#include <deal.II/fe/fe.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/std_cxx14/memory.h>

DEAL_II_NAMESPACE_OPEN


template <class PolynomialType, int dim=PolynomialType::dimension, int spacedim=dim>
class FE_Poly : public FiniteElement<dim,spacedim>
{
public:
  FE_Poly (const PolynomialType &poly_space,
           const FiniteElementData<dim> &fe_data,
           const std::vector<bool> &restriction_is_additive_flags,
           const std::vector<ComponentMask> &nonzero_components);

  unsigned int get_degree () const;

  // for documentation, see the FiniteElement base class
  virtual
  UpdateFlags
  requires_update_flags (const UpdateFlags update_flags) const;

  std::vector<unsigned int> get_poly_space_numbering() const;

  std::vector<unsigned int> get_poly_space_numbering_inverse() const;

  virtual double shape_value (const unsigned int i,
                              const Point<dim> &p) const;

  virtual double shape_value_component (const unsigned int i,
                                        const Point<dim> &p,
                                        const unsigned int component) const;

  virtual Tensor<1,dim> shape_grad (const unsigned int  i,
                                    const Point<dim>   &p) const;

  virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
                                              const Point<dim> &p,
                                              const unsigned int component) const;

  virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
                                         const Point<dim> &p) const;

  virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
                                                   const Point<dim> &p,
                                                   const unsigned int component) const;

  virtual Tensor<3,dim> shape_3rd_derivative (const unsigned int  i,
                                              const Point<dim>   &p) const;

  virtual Tensor<3,dim> shape_3rd_derivative_component (const unsigned int i,
                                                        const Point<dim>   &p,
                                                        const unsigned int component) const;

  virtual Tensor<4,dim> shape_4th_derivative (const unsigned int  i,
                                              const Point<dim>   &p) const;

  virtual Tensor<4,dim> shape_4th_derivative_component (const unsigned int i,
                                                        const Point<dim>   &p,
                                                        const unsigned int component) const;

protected:
  /*
   * 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
  {
    // generate a new data object and
    // initialize some fields
    auto data = std_cxx14::make_unique<InternalData>();
    data->update_each = requires_update_flags(update_flags);

    const unsigned int n_q_points = quadrature.size();

    // initialize some scratch arrays. we need them for the underlying
    // polynomial to put the values and derivatives of shape functions
    // to put there, depending on what the user requested
    std::vector<double> values(update_flags & update_values ?
                               this->dofs_per_cell : 0);
    std::vector<Tensor<1,dim> > grads(update_flags & update_gradients ?
                                      this->dofs_per_cell : 0);
    std::vector<Tensor<2,dim> > grad_grads(update_flags & update_hessians ?
                                           this->dofs_per_cell : 0);
    std::vector<Tensor<3,dim> > third_derivatives(update_flags & update_3rd_derivatives ?
                                                  this->dofs_per_cell : 0);
    std::vector<Tensor<4,dim> > fourth_derivatives;   // won't be needed, so leave empty

    // now also initialize fields the fields of this class's own
    // temporary storage, depending on what we need for the given
    // update flags.
    //
    // there is one exception from the rule: if we are dealing with
    // cells (i.e., if this function is not called via
    // get_(sub)face_data()), then we can already store things in the
    // final location where FEValues::reinit() later wants to see
    // things. we then don't need the intermediate space. we determine
    // whether we are on a cell by asking whether the number of
    // elements in the output array equals the number of quadrature
    // points (yes, it's a cell) or not (because in that case the
    // number of quadrature points we use here equals the number of
    // quadrature points summed over *all* faces or subfaces, whereas
    // the number of output slots equals the number of quadrature
    // points on only *one* face)
    if ((update_flags & update_values)
        &&
        !((output_data.shape_values.n_rows() > 0)
          &&
          (output_data.shape_values.n_cols() == n_q_points)))
      data->shape_values.reinit (this->dofs_per_cell, n_q_points);

    if (update_flags & update_gradients)
      data->shape_gradients.reinit (this->dofs_per_cell, n_q_points);

    if (update_flags & update_hessians)
      data->shape_hessians.reinit (this->dofs_per_cell, n_q_points);

    if (update_flags & update_3rd_derivatives)
      data->shape_3rd_derivatives.reinit (this->dofs_per_cell, n_q_points);

    // next already fill those fields of which we have information by
    // now. note that the shape gradients are only those on the unit
    // cell, and need to be transformed when visiting an actual cell
    if (update_flags & (update_values | update_gradients
                        | update_hessians | update_3rd_derivatives) )
      for (unsigned int i=0; i<n_q_points; ++i)
        {
          poly_space.compute(quadrature.point(i),
                             values, grads, grad_grads,
                             third_derivatives,
                             fourth_derivatives);

          // the values of shape functions at quadrature points don't change.
          // consequently, write these values right into the output array if
          // we can, i.e., if the output array has the correct size. this is
          // the case on cells. on faces, we already precompute data on *all*
          // faces and subfaces, but we later on copy only a portion of it
          // into the output object; in that case, copy the data from all
          // faces into the scratch object
          if (update_flags & update_values)
            if (output_data.shape_values.n_rows() > 0)
              {
                if (output_data.shape_values.n_cols() == n_q_points)
                  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
                    output_data.shape_values[k][i] = values[k];
                else
                  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
                    data->shape_values[k][i] = values[k];
              }

          // for everything else, derivatives need to be transformed,
          // so we write them into our scratch space and only later
          // copy stuff into where FEValues wants it
          if (update_flags & update_gradients)
            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
              data->shape_gradients[k][i] = grads[k];

          if (update_flags & update_hessians)
            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
              data->shape_hessians[k][i] = grad_grads[k];

          if (update_flags & update_3rd_derivatives)
            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
              data->shape_3rd_derivatives[k][i] = third_derivatives[k];
        }
    return std::move(data);
  }

  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;

  virtual
  void
  fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator           &cell,
                       const unsigned int                                                   face_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;

  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;

  class InternalData : public FiniteElement<dim,spacedim>::InternalDataBase
  {
  public:
    Table<2,double> shape_values;

    Table<2,Tensor<1,dim> > shape_gradients;

    Table<2,Tensor<2,dim> > shape_hessians;

    Table<2,Tensor<3,dim> > shape_3rd_derivatives;
  };

  void
  correct_third_derivatives (internal::FEValuesImplementation::FiniteElementRelatedData<dim,spacedim>       &output_data,
                             const internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>       &mapping_data,
                             const unsigned int                                                n_q_points,
                             const unsigned int                                                dof) const;

  PolynomialType poly_space;
};

DEAL_II_NAMESPACE_CLOSE

#endif