fe_poly_tensor.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2017 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_tensor_h
#define dealii__fe_poly_tensor_h


#include <deal.II/base/config.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/fe/fe.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/thread_management.h>


DEAL_II_NAMESPACE_OPEN

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

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

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

  // documentation inherited from the base class
  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;

  // 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;

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

  // 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;

protected:
  MappingType mapping_type;


  /* 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
  typename FiniteElement<dim,spacedim>::InternalDataBase *
  get_data(const UpdateFlags                                                    update_flags,
           const Mapping<dim,spacedim>                                         &/*mapping*/,
           const Quadrature<dim>                                               &quadrature,
           ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &/*output_data*/) const
  {
    // generate a new data object and
    // initialize some fields
    InternalData *data = new InternalData;
    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->sign_change.resize (this->dofs_per_cell);

    // initialize fields only if really
    // necessary. otherwise, don't
    // allocate memory
    if (update_flags & update_values)
      {
        values.resize (this->dofs_per_cell);
        data->shape_values.reinit (this->dofs_per_cell, n_q_points);
        if (mapping_type != mapping_none)
          data->transformed_shape_values.resize (n_q_points);
      }

    if (update_flags & update_gradients)
      {
        grads.resize (this->dofs_per_cell);
        data->shape_grads.reinit (this->dofs_per_cell, n_q_points);
        data->transformed_shape_grads.resize (n_q_points);

        if ( (mapping_type == mapping_raviart_thomas)
             ||
             (mapping_type == mapping_piola)
             ||
             (mapping_type == mapping_nedelec)
             ||
             (mapping_type == mapping_contravariant))
          data->untransformed_shape_grads.resize(n_q_points);
      }

    if (update_flags & update_hessians)
      {
        grad_grads.resize (this->dofs_per_cell);
        data->shape_grad_grads.reinit (this->dofs_per_cell, n_q_points);
        data->transformed_shape_hessians.resize (n_q_points);
        if ( mapping_type != mapping_none )
          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.compute(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->dofs_per_cell; ++i)
                  data->shape_values[i][k] = values[i];
              else
                for (unsigned int i=0; i<this->dofs_per_cell; ++i)
                  {
                    Tensor<1,dim> add_values;
                    for (unsigned int j=0; j<this->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->dofs_per_cell; ++i)
                  data->shape_grads[i][k] = grads[i];
              else
                for (unsigned int i=0; i<this->dofs_per_cell; ++i)
                  {
                    Tensor<2,dim> add_grads;
                    for (unsigned int j=0; j<this->dofs_per_cell; ++j)
                      add_grads += inverse_node_matrix(j,i) * grads[j];
                    data->shape_grads[i][k] = add_grads;
                  }
            }

          if (update_flags & update_hessians)
            {
              if (inverse_node_matrix.n_cols() == 0)
                for (unsigned int i=0; i<this->dofs_per_cell; ++i)
                  data->shape_grad_grads[i][k] = grad_grads[i];
              else
                for (unsigned int i=0; i<this->dofs_per_cell; ++i)
                  {
                    Tensor<3,dim> add_grad_grads;
                    for (unsigned int j=0; j<this->dofs_per_cell; ++j)
                      add_grad_grads += inverse_node_matrix(j,i) * grad_grads[j];
                    data->shape_grad_grads[i][k] = add_grad_grads;
                  }
            }

        }
    return 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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                  const typename FiniteElement<dim,spacedim>::InternalDataBase        &fe_internal,
                  ::internal::FEValues::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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                       const typename FiniteElement<dim,spacedim>::InternalDataBase        &fe_internal,
                       ::internal::FEValues::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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                          const typename FiniteElement<dim,spacedim>::InternalDataBase        &fe_internal,
                          ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const;

  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>                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;
  };



  PolynomialType 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