fe_enriched.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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_enriched_h
#define dealii__fe_enriched_h

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

// We require c++14 here even though we only need features that are part of
// c++11, but it turns out that gcc 4.6.x and 4.7.x don't support c++11
// features like delegating constructors.

#ifdef DEAL_II_WITH_CXX14
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/base/function.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/quadrature.h>

#include <vector>
#include <utility>
#include <numeric>

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim=dim>
class FE_Enriched : public FiniteElement<dim, spacedim>
{
public:

  FE_Enriched (const FiniteElement<dim,spacedim> &fe_base,
               const FiniteElement<dim,spacedim> &fe_enriched,
               const Function<spacedim>      *enrichment_function);

  FE_Enriched (const FiniteElement<dim,spacedim> &fe_base);

  FE_Enriched (const FiniteElement<dim,spacedim> *fe_base,
               const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,
               const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions);

private:
  FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
               const std::vector< unsigned int > &multiplicities,
               const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions);
public:

  virtual FiniteElement<dim,spacedim> *clone() const;

  virtual
  UpdateFlags
  requires_update_flags (const UpdateFlags update_flags) const;

  virtual std::string get_name () const;

  virtual const FiniteElement<dim,spacedim> &
  base_element (const unsigned int index) const;

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

  virtual const FullMatrix<double> &
  get_restriction_matrix (const unsigned int child,
                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;

  virtual const FullMatrix<double> &
  get_prolongation_matrix (const unsigned int child,
                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;


  virtual bool hp_constraints_are_implemented () const;

  virtual void
  get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                                 FullMatrix<double>       &matrix) const;

  virtual void
  get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                                    const unsigned int        subface,
                                    FullMatrix<double>       &matrix) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

  virtual
  FiniteElementDomination::Domination
  compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;


  const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > >
  get_enrichments() const;

  const FESystem<dim,spacedim> &
  get_fe_system() const;

protected:

  class InternalData : public FiniteElement<dim,spacedim>::InternalDataBase
  {
  public:
    struct EnrichmentValues
    {
      std::vector<double> values;
      std::vector<Tensor<1,spacedim> > gradients;
      std::vector<SymmetricTensor<2, spacedim> > hessians;
    };

    InternalData ( std::unique_ptr<typename FESystem<dim,spacedim>::InternalData> fesystem_data);

    typename FiniteElement<dim,spacedim>::InternalDataBase &
    get_fe_data (const unsigned int base_no) const;

    internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
    get_fe_output_object (const unsigned int base_no) const;

    std::unique_ptr<typename FESystem<dim,spacedim>::InternalData> fesystem_data;

    mutable std::vector<std::vector<EnrichmentValues> > enrichment;
  };

  std::vector<std::vector<std::vector<unsigned int> > > base_no_mult_local_enriched_dofs;

  const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > enrichments;

  const bool is_enriched;

  template <int dim_1>
  typename FiniteElement<dim,spacedim>::InternalDataBase *
  setup_data (std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fes_data,
              const UpdateFlags      flags,
              const Quadrature<dim_1> &quadrature) const;

  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
  get_data (const UpdateFlags      flags,
            const Mapping<dim,spacedim>    &mapping,
            const Quadrature<dim> &quadrature,
            ::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const;

  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
  get_face_data (const UpdateFlags      update_flags,
                 const Mapping<dim,spacedim>    &mapping,
                 const Quadrature<dim-1> &quadrature,
                 ::internal::FEValues::FiniteElementRelatedData< dim, spacedim >        &output_data) const;

  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
  get_subface_data (const UpdateFlags      update_flags,
                    const Mapping<dim,spacedim>    &mapping,
                    const Quadrature<dim-1> &quadrature,
                    ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const;

  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;

private:
  void initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
                   const std::vector<unsigned int> &multiplicities);

  FESystem<dim,spacedim> fe_system;

  template <int dim_1>
  void
  multiply_by_enrichment (const Quadrature<dim_1>                                       &quadrature,
                          const InternalData                                            &fe_data,
                          const internal::FEValues::MappingRelatedData<dim,spacedim>    &mapping_data,
                          const typename Triangulation< dim, spacedim >::cell_iterator  &cell,
                          internal::FEValues::FiniteElementRelatedData<dim,spacedim>    &output_data) const;
};

//}
DEAL_II_NAMESPACE_CLOSE

#endif // CXX14

#endif // dealii__fe_enriched_h