fe_enriched.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_fe_enriched_h
#define dealii_fe_enriched_h

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

#include <deal.II/base/function.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/symmetric_tensor.h>

#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_update_flags.h>

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

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 std::unique_ptr<FiniteElement<dim, spacedim>>
  clone() const override;

  virtual UpdateFlags
  requires_update_flags(const UpdateFlags update_flags) const override;

  virtual std::string
  get_name() const override;

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

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

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

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


  virtual bool
  hp_constraints_are_implemented() const override;

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

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

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

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

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

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


  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::FEValuesImplementation::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>
  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
  setup_data(
    std::unique_ptr<typename FESystem<dim, spacedim>::InternalData> fes_data,
    const UpdateFlags                                               flags,
    const Quadrature<dim_1> &quadrature) const;

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

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

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

  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;

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

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

  const std::unique_ptr<const FESystem<dim, spacedim>> fe_system;

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

//}
DEAL_II_NAMESPACE_CLOSE

#endif // dealii_fe_enriched_h