fe_enriched.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2019 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_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_update_flags.h>

#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>

#include <map>
#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_domination(const FiniteElement<dim, spacedim> &fe_other,
                         const unsigned int codim = 0) const override final;



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



namespace ColorEnriched
{
  template <int dim, int spacedim = dim>
  using predicate_function = std::function<bool(
    const typename Triangulation<dim, spacedim>::cell_iterator &)>;

  namespace internal
  {
    template <int dim, int spacedim>
    bool
    find_connection_between_subdomains(
      const hp::DoFHandler<dim, spacedim> &    dof_handler,
      const predicate_function<dim, spacedim> &predicate_1,
      const predicate_function<dim, spacedim> &predicate_2);

    template <int dim, int spacedim>
    unsigned int
    color_predicates(
      const hp::DoFHandler<dim, spacedim> &                 dof_handler,
      const std::vector<predicate_function<dim, spacedim>> &predicates,
      std::vector<unsigned int> &                           predicate_colors);

    template <int dim, int spacedim>
    void
    set_cellwise_color_set_and_fe_index(
      hp::DoFHandler<dim, spacedim> &                       dof_handler,
      const std::vector<predicate_function<dim, spacedim>> &predicates,
      const std::vector<unsigned int> &                     predicate_colors,
      std::map<unsigned int, std::map<unsigned int, unsigned int>>
        &                                  cellwise_color_predicate_map,
      std::vector<std::set<unsigned int>> &fe_sets);

    template <int dim, int spacedim>
    void
    make_colorwise_enrichment_functions(
      const unsigned int                                      n_colors,
      const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments,
      const std::map<unsigned int, std::map<unsigned int, unsigned int>>
        &cellwise_color_predicate_map,
      std::vector<std::function<const Function<spacedim> *(
        const typename Triangulation<dim, spacedim>::cell_iterator &)>>
        &color_enrichments);


    template <int dim, int spacedim>
    void
    make_fe_collection_from_colored_enrichments(
      const unsigned int n_colors,
      const std::vector<std::set<unsigned int>>
        &fe_sets, // total list of color sets possible
      const std::vector<std::function<const Function<spacedim> *(
        const typename Triangulation<dim, spacedim>::cell_iterator &)>>
        &color_enrichments, // color wise enrichment functions
      const FiniteElement<dim, spacedim> &fe_base, // basic FE element
      const FiniteElement<dim, spacedim>
        &fe_enriched, // FE multiplied by enrichment function
      const FE_Nothing<dim, spacedim> &fe_nothing,
      hp::FECollection<dim, spacedim> &fe_collection);
  } // namespace internal



  template <int dim, int spacedim = dim>
  struct Helper
  {
    Helper(const FiniteElement<dim, spacedim> &                    fe_base,
           const FiniteElement<dim, spacedim> &                    fe_enriched,
           const std::vector<predicate_function<dim, spacedim>> &  predicates,
           const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments);

    const hp::FECollection<dim, spacedim> &
    build_fe_collection(hp::DoFHandler<dim, spacedim> &dof_handler);

  private:
    hp::FECollection<dim, spacedim> fe_collection;

    const FiniteElement<dim, spacedim> &fe_base;

    const FiniteElement<dim, spacedim> &fe_enriched;

    const FE_Nothing<dim, spacedim> fe_nothing;

    const std::vector<predicate_function<dim, spacedim>> predicates;

    const std::vector<std::shared_ptr<Function<spacedim>>> enrichments;

    using cell_iterator_function = std::function<const Function<spacedim> *(
      const typename Triangulation<dim, spacedim>::cell_iterator &)>;

    std::vector<cell_iterator_function> color_enrichments;

    std::vector<unsigned int> predicate_colors;

    unsigned int n_colors;

    std::map<unsigned int, std::map<unsigned int, unsigned int>>
      cellwise_color_predicate_map;

    std::vector<std::set<unsigned int>> fe_sets;
  };
} // namespace ColorEnriched

DEAL_II_NAMESPACE_CLOSE

#endif // dealii_fe_enriched_h