]> https://gitweb.dealii.org/ - dealii.git/commitdiff
added Helper function to create FECollection for overlapping sub-domains with individ...
authornivesh <nidhinivesh.d@gmail.com>
Tue, 29 May 2018 09:21:23 +0000 (11:21 +0200)
committernivesh <nidhinivesh.d@gmail.com>
Wed, 20 Jun 2018 18:34:25 +0000 (20:34 +0200)
20 files changed:
doc/doxygen/images/3source_fe_indices.png [new file with mode: 0644]
doc/news/changes/major/20180420NiveshD [new file with mode: 0644]
include/deal.II/fe/fe_enriched.h
include/deal.II/fe/fe_enriched.templates.h [new file with mode: 0644]
source/fe/fe_enriched.cc
source/fe/fe_enriched.inst.in
tests/fe/fe_enriched_color_01.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_01.output [new file with mode: 0644]
tests/fe/fe_enriched_color_02.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_02.with_mpi=true.output [new file with mode: 0644]
tests/fe/fe_enriched_color_03.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_03.with_mpi=true.output [new file with mode: 0644]
tests/fe/fe_enriched_color_04.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_04.with_mpi=true.output [new file with mode: 0644]
tests/fe/fe_enriched_color_05.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_05.with_mpi=true.output [new file with mode: 0644]
tests/fe/fe_enriched_color_06.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output [new file with mode: 0644]
tests/fe/fe_enriched_color_07.cc [new file with mode: 0644]
tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/3source_fe_indices.png b/doc/doxygen/images/3source_fe_indices.png
new file mode 100644 (file)
index 0000000..f07c0b9
Binary files /dev/null and b/doc/doxygen/images/3source_fe_indices.png differ
diff --git a/doc/news/changes/major/20180420NiveshD b/doc/news/changes/major/20180420NiveshD
new file mode 100644 (file)
index 0000000..4f63595
--- /dev/null
@@ -0,0 +1,9 @@
+New: The new class ColorEnriched::Helper constructs an hp::FECollection, a 
+collection of finite element objects used by hp::DoFHandler, in a domain 
+with multiple, possibly overlapping, sub-domains with individual
+enrichment functions. Note that the overlapping regions may have
+multiple enrichment functions associated with them. This is implemented
+using a general constructor of FE_Enriched which allows different
+enrichment functions.
+<br>
+(Nivesh Dommaraju, 2018/04/20)
index 55e66f67201e0ec328445d8485837aee5c73f974..67335a223e1fc62be9d023e01228425fa9f47603 100644 (file)
 
 #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>
@@ -691,6 +696,315 @@ private:
       &output_data) const;
 };
 
+
+
+/**
+ * This namespace consists of a class needed to create a collection
+ * of FE_Enriched finite elements (hp::FECollection)
+ * to be used with hp::DoFHandler in a domain with multiple, possibly
+ * overlapping, sub-domains with individual enrichment functions.
+ *
+ * To create hp::FECollection, a graph coloring algorithm is used to assign
+ * colors to enrichment functions before creating hp::FECollection. Hence the
+ * name.
+ */
+namespace ColorEnriched
+{
+  /**
+   * An alias template for predicate function which returns a
+   * boolean for a Triangulation<dim,spacedim>::cell_iterator object.
+   *
+   * This is used by helper functions and in the implementation of
+   * ColorEnriched::Helper class.
+   */
+  template <int dim, int spacedim = dim>
+  using predicate_function = std::function<bool(
+    const typename Triangulation<dim, spacedim>::cell_iterator &)>;
+
+  /**
+   * ColorEnriched::Helper class creates a collection of FE_Enriched finite
+   * elements (hp::FECollection) to be used with hp::DoFHandler in a domain
+   * with multiple, possibly overlapping, sub-domains with individual
+   * enrichment functions. Note that the overlapping regions may have
+   * multiple enrichment functions associated with them. This is implemented
+   * using a general constructor of FE_Enriched object which allows different
+   * enrichment functions.
+   *
+   * Consider a domain with multiple enriched sub-domains
+   * which are disjoint i.e. not connected with each other.
+   * To ensure $C^0$ continuity at the interface between
+   * the enriched sub-domain (characterized by a single enrichment
+   * function) and the non-enriched domain, we can use an FE_Enriched
+   * object in the enriched sub-domain and in the non-enriched domain
+   * a standard finite element (eg: FE_Q) wrapped into an FE_Enriched
+   * object (which internally uses a dominating FE_Nothing object).
+   * Refer to the documentation on FE_Enriched for more
+   * information on this. It is to be noted that an FE_Enriched
+   * object is constructed using a base FE
+   * (FiniteElement objects) and one or more
+   * enriched FEs. FE_Nothing is a dummy enriched FE.
+   *
+   * The situation becomes more
+   * complicated when two enriched sub-domains
+   * share an interface. When the number of enrichment functions are
+   * same for the sub-domains, FE_Enriched object of one sub-domain
+   * is constructed such that each enriched FE is paired (figuratively) with a
+   * FE_Nothing in the FE_Enriched object of the other sub-domain.
+   * For example, let the FEs fe_enr1 and fe_enr2, which will be
+   * used with enrichment functions, correspond
+   * to the two sub-domains. Then the FE_Enriched objects of the two
+   * sub-domains are built using
+   * [fe_base, fe_enr1, fe_nothing] and
+   * [fe_base, fe_nothing, fe_enr2] respectively.
+   * Note that the size of the vector of enriched FEs
+   * (used in FE_Enriched constructor) is equal to 2, the
+   * same as the number of enrichment functions. When the number of enrichment
+   * functions is not the same, additional enriched FEs are paired
+   * with FE_Nothing. This ensures that the enriched DOF's at the interface
+   * are set to zero by the DoFTools::make_hanging_node_constraints() function.
+   * Using these two strategies, we construct the appropriate FE_Enriched
+   * using the general constructor. Note that this is
+   * done on a mesh without hanging nodes.
+   *
+   * Now consider a domain with multiple sub-domains which may share
+   * an interface with each other. As discussed previously,
+   * the number of enriched FEs in the FE_Enriched object of each
+   * sub-domain needs to be equal to the number of sub-domains. This is because
+   * we are not using the information of how the domains are connected
+   * and any sub-domain may share interface with any other sub-domain (not
+   * considering overlaps for now!). However, in general, a given sub-domain
+   * shares an interface only with a few sub-domains. This warrants
+   * the use of a graph coloring algorithm to reduce
+   * the size of the vector of enriched FEs
+   * (used in the FE_Enriched constructor). By giving the sub-domains
+   * that share no interface the same color, a single 'std::function'
+   * that returns different enrichment functions for each
+   * sub-domain can be constructed. Then the size of the vector of enriched
+   * FEs is equal to the number of different colors
+   * used for predicates (or sub-domains).
+   *
+   * @note The graph coloring function, SparsityTools::color_sparsity_pattern,
+   * used for assigning colors to the sub-domains
+   * needs MPI (use Utilities::MPI::MPI_InitFinalize to initialize MPI
+   * and the necessary Zoltan setup).
+   * The coloring function, based on Zoltan, is a parallel coloring
+   * algorithm but is used in serial by SparsityTools::color_sparsity_pattern.
+   *
+   * Construction of the Helper class needs a base FiniteElement @p fe_base,
+   * an enriched FiniteElement @p fe_enriched (used for all the
+   * enrichment functions), a vector of predicate
+   * functions (used to define sub-domains) as well as the corresponding
+   * enrichment functions. The FECollection object, a collection of FE_Enriched
+   * objects to be used with an hp::DoFHandler object, can be retrieved
+   * using the member function build_fe_collection which also modifies the
+   * active FE indices of the hp::DoFHandler object (provided as an argument
+   * to the build_fe_collection function).
+   *
+   * <h3>Simple example</h3>
+   * Consider a domain with three sub-domains defined by predicate functions.
+   * Different cells are associated with FE indices as shown in the following
+   * image. The three equal-sized square-shaped sub-domains 'a', 'b'
+   * and 'c' can be seen. The predicates associated with these sub-domains
+   * are also labeled 'a', 'b' and 'c'.
+   * The subdomains 'a' and 'b' intersect with cell labeled with FE
+   * index 3. The cells in 'c' are labeled with FE
+   * index 1. As can be seen, connections exist between 'a' and 'b',
+   * 'b' and 'c' but 'a' and 'c' are not connected.
+   *
+   * \htmlonly <style>div.image
+   * img[src="3source_fe_indices.png"]{width:25%;}</style> \endhtmlonly
+   * @image html 3source_fe_indices.png "Active FE indices"
+   *
+   * As discussed before, the colors of predicates are alloted using
+   * the graph coloring algorithm. Each predicate is a node in the graph and if
+   * two sub-domains share an interface, the corresponding predicates
+   * should be given different colors.
+   * (The predicate colors are different from what is shown
+   * in the image. The colors in the image are as per FE indices).
+   * Predicates 'a' and 'c' can be given the same color since they
+   * are not connected but the color given to 'b' has to be different from
+   * 'a' and 'c'.
+   *
+   * The name of finite element at an index (i) of @p fe_collection
+   * (hp::FECollection) can be obtained by
+   * <code>fe_collection[index].get_name()</code> and is
+   * show in the table below. Note that all the FE_Enriched elements
+   * are of the same size and FE_Nothing<2>(dominating) is used as
+   * discussed before.
+   *
+   * <table>
+   * <tr>
+   * <th>Active FE index</th>
+   * <th>FiniteElement name</code> </th>
+   * </tr>
+   * <tr>
+   * <td>0</td>
+   * <td><code>FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)]</code></td>
+   * </tr>
+   * <tr>
+   * <td>1</td>
+   * <td><code>FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)]</code></td>
+   * </tr>
+   * <tr>
+   * <td>2</td>
+   * <td><code>FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)]</code></td>
+   * </tr>
+   * <tr>
+   * <td>3</td>
+   * <td><code>FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)]</code></td>
+   * </tr>
+   * </table>
+   *
+   * The internal data members used by this class need to be available when the
+   * problem is solved. This can be ensured by declaring the object static,
+   * which is deallocated only when the program terminates. An alternative
+   * would be to use it as a data member of the containing class. Since vector
+   * of predicates and enrichment functions may not be available while
+   * constructing the Helper, a 'std::shared_ptr' to Helper object can be used
+   * and constructed when the predicates and enrichment functions are
+   * available.
+   *
+   * <h3>Example usage:</h3>
+   * @code
+   * FE_Q<dim> fe_base(2);
+   * FE_Q<dim> fe_enriched(1);
+   * std::vector< predicate_function<dim> > predicates;
+   * std::vector< std::shared_ptr<Function<dim>> > enrichments;
+   *
+   * Triangulation<dim>  triangulation;
+   * hp::DoFHandler<dim> dof_handler(triangulation);
+   *
+   * static ColorEnriched::Helper<dim> FE_helper(fe_base,
+   *                                             fe_enriched,
+   *                                             predicates,
+   *                                             enrichments);
+   * const hp::FECollection<dim>&
+   * fe_collection(FE_helper.build_fe_collection(dof_handler));
+   * @endcode
+   */
+  template <int dim, int spacedim = dim>
+  struct Helper
+  {
+    /**
+     * Constructor for Helper class.
+     *
+     * @param fe_base A base FiniteElement
+     * @param fe_enriched An enriched FiniteElement
+     * @param predicates std::vector of predicates defining the sub-domains.
+     * <code>@p predicates[i]</code> returns true for a cell if it
+     * belongs to a sub-domain with index (i).
+     * @param enrichments std::vector of enrichment functions
+     */
+    Helper(const FE_Q<dim, spacedim> &                             fe_base,
+           const FE_Q<dim, spacedim> &                             fe_enriched,
+           const std::vector<predicate_function<dim, spacedim>> &  predicates,
+           const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments);
+
+    /**
+     * Prepares an hp::DoFHandler object. The active FE indices of
+     * mesh cells are initialized to work with
+     * ColorEnriched::Helper<dim,spacedim>::fe_collection.
+     *
+     * Returns an hp::FECollection object.
+     *
+     * @param dof_handler an hp::DoFHandler object
+     * @return hp::FECollection, a collection of
+     * finite elements needed by @p dof_handler.
+     */
+    const hp::FECollection<dim, spacedim> &
+    build_fe_collection(hp::DoFHandler<dim, spacedim> &dof_handler);
+
+  private:
+    /**
+     * Contains a collection of FiniteElement objects needed by an
+     * hp::DoFHandler object.
+     */
+    hp::FECollection<dim, spacedim> fe_collection;
+
+    /**
+     * A base FiniteElement used for constructing FE_Enriched
+     * object required by ColorEnriched::Helper<dim,spacedim>::fe_collection.
+     */
+    const FE_Q<dim, spacedim> &fe_base;
+
+    /**
+     * An enriched FiniteElement used for constructing FE_Enriched
+     * object required by ColorEnriched::Helper<dim,spacedim>::fe_collection.
+     */
+    const FE_Q<dim, spacedim> &fe_enriched;
+
+    /**
+     * A finite element with zero degrees of freedom used for
+     * constructing FE_Enriched object required by
+     * ColorEnriched::Helper<dim,spacedim>::fe_collection
+     */
+    const FE_Nothing<dim, spacedim> fe_nothing;
+
+    /**
+     * std::vector of predicates defining the sub-domains.
+     * <code>@p predicates[i]</code> returns true for a cell if it
+     * belongs to a sub-domain with index (i).
+     */
+    const std::vector<predicate_function<dim, spacedim>> predicates;
+
+    /**
+     * std::vector of enrichment functions corresponding
+     * to the predicates. These are needed while constructing
+     * ColorEnriched::Helper<dim,spacedim>::fe_collection.
+     */
+    const std::vector<std::shared_ptr<Function<spacedim>>> enrichments;
+
+    /**
+     * An alias template for any callable target such as functions, lambda
+     * expressions, function objects that take a
+     * Triangulation<dim,spacedim>::cell_iterator
+     * and return a pointer to Function<dim>. This is used to define
+     * Helper<dim,spacedim>::color_enrichments
+     * which returns an enrichment function
+     * for a cell in Triangulation<dim,spacedim>.
+     */
+    using cell_iterator_function = std::function<const Function<spacedim> *(
+      const typename Triangulation<dim, spacedim>::cell_iterator &)>;
+
+    /**
+     * std::vector of functions that take in a cell
+     * and return a function pointer. These are needed while constructing
+     * fe_collection.
+     *
+     * color_enrichments[i](cell_iterator) calls the correct enrichment function
+     * (i.e. whose corresponding predicate has the color i) for the cell.
+     * Note that this call to cell_iterator returns the enrichment function
+     * which is a pointer to class derived from Function<dim>.
+     */
+    std::vector<cell_iterator_function> color_enrichments;
+
+    /**
+     * std::vector of colors (unsigned int) associated
+     * with each sub-domain. No two connected sub-domains (i.e. sub-domains that
+     * share a vertex) have the same color.
+     */
+    std::vector<unsigned int> predicate_colors;
+
+    /**
+     * Total number of different colors in predicate_colors
+     */
+    unsigned int num_colors;
+
+    /**
+     * A map of maps used to associate
+     * a cell with a map that in turn associates colors of active predicates in
+     * the cell with corresponding predicate ids.
+     */
+    std::map<unsigned int, std::map<unsigned int, unsigned int>>
+      cellwise_color_predicate_map;
+
+    /**
+     * A vector of different possible color sets for a given set of
+     * predicates and hp::DoFHandler object
+     */
+    std::vector<std::set<unsigned int>> fe_sets;
+  };
+} // namespace ColorEnriched
 //}
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/include/deal.II/fe/fe_enriched.templates.h b/include/deal.II/fe/fe_enriched.templates.h
new file mode 100644 (file)
index 0000000..674aebe
--- /dev/null
@@ -0,0 +1,212 @@
+// ---------------------------------------------------------------------
+//
+// 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_templates_h
+#define dealii_fe_enriched_templates_h
+
+#include <deal.II/fe/fe_enriched.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ColorEnriched
+{
+  namespace internal
+  {
+    /**
+     * Returns true if there is a connection between subdomains in the mesh
+     * associated with @p hp::DoFHandler i.e., if the subdomains share at least
+     * a vertex. The two subdomains are defined by predicates provided by
+     * @p predicate_1 and @p predicate_2. A predicate is a function (or
+     * object of a type with an operator()) which takes in a cell iterator and
+     * gives a boolean. It is said to be active in a cell if it returns true.
+     *
+     * An example of a custom predicate is one that checks the distance from a
+     * fixed point. Note that the operator() takes in a cell iterator. Using the
+     * constructor, the fixed point and the distance can be chosen.
+     * @code
+     * <int dim>
+     * struct predicate
+     * {
+     *     predicate(const Point<dim> p, const int radius)
+     *     :p(p),radius(radius){}
+     *
+     *     template <class Iterator>
+     *     bool operator () (const Iterator &i)
+     *     {
+     *         return ( (i->center() - p).norm() < radius);
+     *     }
+     *
+     * private:
+     *     Point<dim> p;
+     *     int radius;
+     *
+     * };
+     * @endcode
+     * and then the function can be used as follows to find if the subdomains
+     * are connected.
+     * @code
+     * find_connection_between_subdomains
+     * (dof_handler,
+     *  predicate<dim>(Point<dim>(0,0), 1)
+     *  predicate<dim>(Point<dim>(2,2), 1));
+     * @endcode
+     *
+     * @param[in] hp::DoFHandler object
+     * @param[in] predicate_1 A function (or object of a type with an
+     * operator()) defining the subdomain 1. The function takes in a cell and
+     * returns a boolean.
+     * @param[in] predicate_2 Same as @p predicate_1 but defines subdomain 2.
+     * @return A boolean "true" if the subdomains share atleast a vertex.
+     */
+    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);
+
+    /**
+     * Assign colors to subdomains using Graph coloring algorithm where each
+     * subdomain is considered as a graph node. Subdomains which are
+     * connected i.e share atleast a vertex have different color. Each subdomain
+     * is defined using a predicate function of @p predicates.
+     *
+     * @param[in] dof_handler a hp::DoFHandler object
+     * @param[in] predicates predicates defining the subdomains
+     * @param[out] predicate_colors Colors (unsigned int) associated with each
+     * subdomain.
+     */
+    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);
+
+    /**
+     * Used to construct data members @p cellwise_color_predicate_map and
+     * @p fe_sets of Helper class. Inputs are hp::DoFHandler object,
+     * vector of predicates and colors associated with them. Before calling
+     * this function, colors can be assigned to predicates (i.e subdomains)
+     * using the function color_predicates.
+     *
+     * Each active FE index has a set of colors associated with it.
+     * A cell with an active FE index i has a set of colors given by
+     * <code>fe_sets[i]</code>. An active FE index with color {a,b}
+     * means that the cell has two active predicates (i.e they return true
+     * for the cell) of color a and b.
+     *
+     * Eg: fe_sets = { {}, {1}, {2}, {1,2} } means
+     * Cells with active FE index 0 have no predicates associated.
+     * Cells with index 1 have a active predicate with color 1.
+     * Cells with index 2 have a active predicate with color 2.
+     * Cells with index 3 have active predicates with color 1 and color 2.
+     *
+     * A map of maps cellwise_color_predicate_map is used to associate
+     * predicate colors in cells with predicate ids. For this purpose, each
+     * cell is given a unique id which is stored in material id for now.
+     * When the grid is refined, material id is inherited to the children, so
+     * map which associates material id with color map will still be relevant.
+     *
+     * Now the color map can be explained with an example. If the cell with
+     * material id 100 has active predicates 4 (color = 1) and 5 (color = 2),
+     * the map will insert pairs (1, 4) and (2, 5) at key 100 (i.e unique id
+     * of cell is mapped with a map which associates color with predicate id).
+     *
+     * @param[in] dof_handler hp::DoFHandler object
+     * @param[in] predicates vector of predicates defining the subdomains.
+     * <code>@p predicates[i]</code> returns true for a cell if it
+     * belongs to subdomain with index i.
+     * @param[in] predicate_colors vector of colors (unsigned int) associated
+     * with each subdomain.
+     * @param[out] cellwise_color_predicate_map A map of maps used to associate
+     * predicate colors in cells with predicate ids.
+     * @param[out] fe_sets a vector of color lists
+     */
+    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);
+
+    /**
+     * A function that returns a vector of enrichment functions corresponding
+     * to a color. The size of the vector is equal to total number of different
+     * colors associated with predicates (i.e subdomains).
+     *
+     * Assume that a cell has a active predicates with ids 4 (color = 1) and
+     * 5 (color = 2). cellwise_color_predicate_map has this information
+     * provided we know the material id.
+     *
+     * Now a call to color_enrichment[1](cell) should in turn call
+     * enrichments[4](cell).
+     *
+     * @param[in] num_colors number of colors for predicates
+     * @param[in] enrichments vector of enrichment functions
+     * @param[in] cellwise_color_predicate_map A map of maps used to associate
+     * predicate colors in cells with predicate ids.
+     * @param[out] color_enrichments A vector of functions that take in cell
+     * and return a function pointer.
+     */
+    template <int dim, int spacedim>
+    void
+    make_colorwise_enrichment_functions(
+      const unsigned int &                                    num_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);
+
+
+    /**
+     * Creates a hp::FECollection object constructed using FE_Enriched
+     * elements which itself is constructed using color enrichment functions
+     * and is of size equal to number of colors.
+     *
+     * @param[in] fe_sets a vector of color lists
+     * @param[in] color_enrichments A vector of functions that take in cell
+     * and return a function pointer.
+     * @param[in] fe_base base FiniteElement
+     * @param[in] fe_enriched enriched FiniteElements
+     * @param[in] fe_nothing a finite element with zero degrees of freedom
+     * @param[out] fe_collection a collection of
+     * finite elements
+     */
+    template <int dim, int spacedim>
+    void
+    make_fe_collection_from_colored_enrichments(
+      const unsigned int &num_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 FE_Q<dim, spacedim> &fe_base, // basic FE element
+      const FE_Q<dim, spacedim>
+        &fe_enriched, // FE element multiplied by enrichment function
+      const FE_Nothing<dim, spacedim> &fe_nothing,
+      hp::FECollection<dim, spacedim> &fe_collection);
+  } // namespace internal
+} // namespace ColorEnriched
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // dealii_fe_enriched_templates_h
index c08b19a20d9a5880d7c11d179f0abddef824ce30..179805213842861b8276fdaaa3351709317c1e44 100644 (file)
 #include <deal.II/base/std_cxx14/memory.h>
 
 #include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_enriched.templates.h>
 #include <deal.II/fe/fe_tools.h>
 
+#include <deal.II/lac/sparsity_tools.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1048,6 +1050,374 @@ FE_Enriched<dim, spacedim>::InternalData::get_fe_output_object(
 }
 
 
+namespace ColorEnriched
+{
+  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)
+    {
+      // Use a vector to mark vertices
+      std::vector<bool> vertices_subdomain_1(
+        dof_handler.get_triangulation().n_vertices(), false);
+
+      // Mark vertices that belong to cells in subdomain 1
+      for (auto cell : dof_handler.active_cell_iterators())
+        if (predicate_1(cell)) // True ==> part of subdomain 1
+          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
+               ++v)
+            vertices_subdomain_1[cell->vertex_index(v)] = true;
+
+      // Find if cells in subdomain 2 and subdomain 1 share vertices.
+      for (auto cell : dof_handler.active_cell_iterators())
+        if (predicate_2(cell)) // True ==> part of subdomain 2
+          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
+               ++v)
+            if (vertices_subdomain_1[cell->vertex_index(v)] == true)
+              {
+                return true;
+              }
+      return false;
+    }
+
+
+
+    template <int dim, int spacedim>
+    unsigned int
+    color_predicates(
+      const hp::DoFHandler<dim, spacedim> &                 mesh,
+      const std::vector<predicate_function<dim, spacedim>> &predicates,
+      std::vector<unsigned int> &                           predicate_colors)
+    {
+      const unsigned int num_indices = predicates.size();
+
+      // Use sparsity pattern to represent connections between subdomains.
+      // Each predicate (i.e a subdomain) is a node in the graph.
+      DynamicSparsityPattern dsp;
+      dsp.reinit(num_indices, num_indices);
+
+      /*
+       * Find connections between subdomains taken two at a time.
+       * If the connection exists, add it to a graph object represented
+       * by dynamic sparsity pattern.
+       */
+      for (unsigned int i = 0; i < num_indices; ++i)
+        for (unsigned int j = i + 1; j < num_indices; ++j)
+          if (internal::find_connection_between_subdomains(
+                mesh, predicates[i], predicates[j]))
+            dsp.add(i, j);
+
+      dsp.symmetrize();
+
+      // Copy dynamic sparsity pattern to sparsity pattern needed by
+      // coloring function
+      SparsityPattern sp_graph;
+      sp_graph.copy_from(dsp);
+      predicate_colors.resize(num_indices);
+
+      // Assign each predicate with a color and return number of colors
+      return SparsityTools::color_sparsity_pattern(sp_graph, 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)
+    {
+      // clear output variables first
+      fe_sets.clear();
+      cellwise_color_predicate_map.clear();
+
+      /*
+       * Add first element of fe_sets which is empty by default. This means that
+       * the default, FE index = 0 is associated with an empty set, since no
+       * predicate is active in these regions.
+       */
+      fe_sets.resize(1);
+
+      /*
+       * Loop through cells and find set of predicate colors associated
+       * with the cell. As an example, a cell with an FE index associated
+       * with colors {a,b} means that predicates active in the cell have
+       * colors a or b.
+       *
+       * Create new active FE index in case of the color
+       * set is not already listed in fe_sets. If the set already exists,
+       * find index of the set in fe_sets. In either case, use the id in
+       * fe_sets to modify cell->active_fe_index.
+       *
+       * Associate each cell_id with a set of pairs. The pair represents
+       * predicate color and the active predicate with that color.
+       * Each color can only correspond to a single predicate since
+       * predicates with the same color correspond to disjoint domains.
+       * This is what the graph coloring in color_predicates
+       * function ensures. The number of pairs is equal to the number
+       * of predicates active in the given cell.
+       */
+      unsigned int map_index(0);
+      auto         cell = dof_handler.begin_active();
+      auto         endc = dof_handler.end();
+      for (; cell != endc; ++cell)
+        {
+          // set default FE index ==> no enrichment and no active predicates
+          cell->set_active_fe_index(0);
+
+          // Give each cell a unique id, which the cellwise_color_predicate_map
+          // can later use to associate a colors of active predicates with
+          // the actual predicate id.
+          //
+          // When the grid is refined, material id is inherited to
+          // children cells. So, the cellwise_color_predicate_map stays
+          // relevant.
+          cell->set_material_id(map_index);
+          std::set<unsigned int> color_list;
+
+          // loop through active predicates for the cell and insert map.
+          // Eg: if the cell with material id 100 has active
+          // predicates 4 (color = 1) and 5 (color = 2), the map will insert
+          // pairs (1, 4) and (2, 5) at key 100 (i.e unique id of cell is
+          // mapped with a map which associates color with predicate id)
+          // Note that color list for the cell would be {1,2}.
+          for (unsigned int i = 0; i < predicates.size(); ++i)
+            {
+              if (predicates[i](cell))
+                {
+                  /*
+                   * create a pair predicate color and predicate id and add it
+                   * to a map associated with each enriched cell
+                   */
+                  auto ret = cellwise_color_predicate_map[map_index].insert(
+                    std::pair<unsigned int, unsigned int>(predicate_colors[i],
+                                                          i));
+
+                  AssertThrow(
+                    ret.second == 1,
+                    ExcMessage("Only one enrichment function per color"));
+
+                  color_list.insert(predicate_colors[i]);
+                }
+            }
+
+
+          /*
+           * check if color combination is already added.
+           * If already added, set the active FE index based on
+           * its index in the fe_sets. If the combination doesn't
+           * exist, add the set to fe_sets and once again set the
+           * active FE index as last index in fe_sets.
+           *
+           * Eg: if cell has color list {1,2} associated and
+           * fe_sets = { {}, {1}, {2} } for now, we need to add a new set {1,2}
+           * to fe_sets and a new active FE index 3 because 0 to 2 FE indices
+           * are already taken by existing sets in fe_sets.
+           */
+          bool found = false;
+          if (!color_list.empty())
+            {
+              for (unsigned int j = 0; j < fe_sets.size(); ++j)
+                {
+                  if (fe_sets[j] == color_list)
+                    {
+                      found = true;
+                      cell->set_active_fe_index(j);
+                      break;
+                    }
+                }
+
+              if (!found)
+                {
+                  fe_sets.push_back(color_list);
+                  cell->set_active_fe_index(fe_sets.size() - 1);
+                }
+            }
+          ++map_index; // map_index should be unique to cells
+        }
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    make_colorwise_enrichment_functions(
+      const unsigned int &                                    num_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)
+    {
+      color_enrichments.clear();
+
+      // Each color should be associated with a single enrichment function
+      // called color enrichment function which calls the correct enrichment
+      // function for a given cell.
+      //
+      // Assume that a cell has a active predicates with ids 4 (color = 1) and
+      // 5 (color = 2). cellwise_color_predicate_map has this information
+      // provided we know the material id.
+      //
+      // Now a call to color_enrichment[1](cell) should in turn call
+      // enrichments[4](cell). That is just what this lambda is doing.
+      color_enrichments.resize(num_colors);
+      for (unsigned int i = 0; i < num_colors; ++i)
+        {
+          color_enrichments[i] =
+            [&, i](const typename Triangulation<dim, spacedim>::cell_iterator
+                     &cell) {
+              unsigned int id = cell->material_id();
+
+              /*
+               * i'th color_enrichment function corresponds to (i+1)'th color.
+               * Since FE_Enriched takes function pointers, we return a
+               * function pointer.
+               */
+              return enrichments[cellwise_color_predicate_map.at(id).at(i + 1)]
+                .get();
+            };
+        }
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    make_fe_collection_from_colored_enrichments(
+      const unsigned int &                       num_colors,
+      const std::vector<std::set<unsigned int>> &fe_sets,
+      const std::vector<std::function<const Function<spacedim> *(
+        const typename Triangulation<dim, spacedim>::cell_iterator &)>>
+        &                              color_enrichments,
+      const FE_Q<dim, spacedim> &      fe_base,
+      const FE_Q<dim, spacedim> &      fe_enriched,
+      const FE_Nothing<dim, spacedim> &fe_nothing,
+      hp::FECollection<dim, spacedim> &fe_collection)
+    {
+      // define dummy function which is associated with FE_Nothing
+      using cell_function = std::function<const Function<spacedim> *(
+        const typename Triangulation<dim, spacedim>::cell_iterator &)>;
+      cell_function dummy_function;
+      dummy_function =
+        [=](const typename Triangulation<dim, spacedim>::cell_iterator &)
+        -> const Function<spacedim> * {
+        AssertThrow(false,
+                    ExcMessage("Called enrichment function for FE_Nothing"));
+        return nullptr;
+      };
+
+      using EnrichmentFunctions_2DVector =
+        std::vector<std::vector<std::function<const Function<spacedim> *(
+          const typename Triangulation<dim, spacedim>::cell_iterator &)>>>;
+
+      // loop through color sets and create FE_enriched element for each
+      // of them provided before calling this function, we have color
+      // enrichment function associated with each color.
+      for (unsigned int color_set_id = 0; color_set_id != fe_sets.size();
+           ++color_set_id)
+        {
+          std::vector<const FiniteElement<dim, spacedim> *> vec_fe_enriched(
+            num_colors, &fe_nothing);
+          EnrichmentFunctions_2DVector functions(num_colors, {dummy_function});
+
+          for (auto it = fe_sets[color_set_id].begin();
+               it != fe_sets[color_set_id].end();
+               ++it)
+            {
+              // Given a color id ( = *it), corresponding color enrichment
+              // function is at index id-1 because color_enrichments is
+              // is indexed from zero.
+              const unsigned int ind = *it - 1;
+
+              AssertIndexRange(ind, vec_fe_enriched.size());
+              AssertIndexRange(ind, functions.size());
+              AssertIndexRange(ind, color_enrichments.size());
+
+              // Assume a active predicate colors {1,2} for a cell.
+              // We then need to create a vector of FE enriched elements
+              // with vec_fe_enriched[0] = vec_fe_enriched[1] = &fe_enriched
+              // which can later be associated with enrichment functions.
+              vec_fe_enriched[ind] = &fe_enriched;
+
+              // color_set_id'th color function is (color_set_id-1)
+              // element of color wise enrichments
+              functions[ind].assign(1, color_enrichments[ind]);
+            }
+
+          AssertDimension(vec_fe_enriched.size(), functions.size());
+
+          FE_Enriched<dim, spacedim> fe_component(
+            &fe_base, vec_fe_enriched, functions);
+          fe_collection.push_back(fe_component);
+        }
+    }
+  } // namespace internal
+
+
+
+  template <int dim, int spacedim>
+  Helper<dim, spacedim>::Helper(
+    const FE_Q<dim, spacedim> &                             fe_base,
+    const FE_Q<dim, spacedim> &                             fe_enriched,
+    const std::vector<predicate_function<dim, spacedim>> &  predicates,
+    const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments) :
+    fe_base(fe_base),
+    fe_enriched(fe_enriched),
+    fe_nothing(1, true),
+    predicates(predicates),
+    enrichments(enrichments)
+  {}
+
+
+
+  template <int dim, int spacedim>
+  const hp::FECollection<dim, spacedim> &
+  Helper<dim, spacedim>::build_fe_collection(
+    hp::DoFHandler<dim, spacedim> &dof_handler)
+  {
+    // color the predicates based on connections between corresponding
+    // subdomains
+    if (predicates.size() != 0)
+      num_colors =
+        internal::color_predicates(dof_handler, predicates, predicate_colors);
+    else
+      num_colors = 0;
+
+    // create color maps and color list for each cell
+    internal::set_cellwise_color_set_and_fe_index(dof_handler,
+                                                  predicates,
+                                                  predicate_colors,
+                                                  cellwise_color_predicate_map,
+                                                  fe_sets);
+    // setup color wise enrichment functions
+    // i'th function corresponds to (i+1) color!
+    internal::make_colorwise_enrichment_functions<dim, spacedim>(
+      num_colors, enrichments, cellwise_color_predicate_map, color_enrichments);
+
+    // make FE_Collection
+    internal::make_fe_collection_from_colored_enrichments(num_colors,
+                                                          fe_sets,
+                                                          color_enrichments,
+                                                          fe_base,
+                                                          fe_enriched,
+                                                          fe_nothing,
+                                                          fe_collection);
+
+    return fe_collection;
+  }
+} // namespace ColorEnriched
+
+
 // explicit instantiations
 #include "fe_enriched.inst"
 
index f80860e9a058ee7c084b1c79825ed2b94ba6127f..41f5a674299ddbc54f4e305569c3f649c4346c2d 100644 (file)
@@ -18,5 +18,60 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
     template class FE_Enriched<deal_II_dimension, deal_II_space_dimension>;
+
+    namespace ColorEnriched \{
+        namespace internal \{
+            template
+            bool find_connection_between_subdomains
+            (const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension>  &dof_handler,
+             const predicate_function<deal_II_dimension, deal_II_space_dimension> &predicate_1,
+             const predicate_function<deal_II_dimension, deal_II_space_dimension> &predicate_2);
+
+
+            template
+            unsigned int color_predicates
+            (const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &dof_handler,
+             const std::vector<predicate_function<deal_II_dimension, deal_II_space_dimension>> &,
+             std::vector<unsigned int> &);
+
+
+            template
+            void set_cellwise_color_set_and_fe_index
+            (hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &dof_handler,
+             const std::vector<predicate_function <deal_II_dimension, deal_II_space_dimension>>
+             &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
+            void make_colorwise_enrichment_functions<deal_II_dimension, deal_II_space_dimension>
+            (const unsigned int &num_colors,
+             const std::vector< std::shared_ptr <Function<deal_II_space_dimension>> >
+             &enrichments,
+             const std::map<unsigned int,std::map<unsigned int, unsigned int> >
+             &cellwise_color_predicate_map,
+             std::vector<std::function<const Function<deal_II_space_dimension>*
+             (const typename Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator &)>>
+             &color_enrichments);
+
+
+            template
+            void make_fe_collection_from_colored_enrichments
+            (const unsigned int &num_colors,
+             const std::vector <std::set<unsigned int>> &fe_sets,
+             const std::vector<std::function<const Function<deal_II_space_dimension>*
+             (const typename Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator &)>>
+             &color_enrichments,
+             const FE_Q<deal_II_dimension, deal_II_space_dimension> &fe_base,
+             const FE_Q<deal_II_dimension, deal_II_space_dimension> &fe_enriched,
+             const FE_Nothing<deal_II_dimension, deal_II_space_dimension> &fe_nothing,
+             hp::FECollection<deal_II_dimension, deal_II_space_dimension> &fe_collection);
+        \}
+        template struct Helper<deal_II_dimension, deal_II_space_dimension>;
+    \}
+
 #endif
   }
diff --git a/tests/fe/fe_enriched_color_01.cc b/tests/fe/fe_enriched_color_01.cc
new file mode 100644 (file)
index 0000000..ddee396
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test ColorEnriched::internal::find_connection_between_subdomains(...)
+// function. Check if the function correctly finds if two subdomains
+// share an edge/node.
+
+#include <deal.II/fe/fe_enriched.templates.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+/*
+ * Construct a class template which finds if a cell
+ * is within region or not, based on distance of cell
+ * center from region center.
+ */
+template <int dim>
+struct predicate_template
+{
+  predicate_template(const Point<dim> p, const int radius) :
+    p(p),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i)
+  {
+    return ((i->center() - p).norm() < radius);
+  }
+
+private:
+  Point<dim> p;
+  int        radius;
+};
+
+
+
+template <int dim>
+void
+test()
+{
+  deallog << "dim = " << dim << std::endl;
+
+  // Construct grid
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  GridGenerator::hyper_cube(triangulation, -20, 20);
+  triangulation.refine_global(4);
+
+  // Construct vector of predicates for 2 and 3 dimensions
+  Assert(dim == 2 || dim == 3, ExcDimensionMismatch2(dim, 2, 3));
+  typedef std::function<bool(
+    const typename Triangulation<dim>::active_cell_iterator &)>
+                                  predicate_function;
+  std::vector<predicate_function> predicates;
+  predicates.resize(5);
+  if (dim == 2)
+    {
+      // Radius set such that every region has 2^2 = 4 cells
+      predicates[1] = predicate_template<dim>(Point<dim>(7.5, 7.5), 2);
+      predicates[2] = predicate_template<dim>(Point<dim>(5, 5), 2);
+      predicates[0] = predicate_template<dim>(Point<dim>(), 2);
+      predicates[3] = predicate_template<dim>(Point<dim>(-5, -5), 2);
+      predicates[4] = predicate_template<dim>(Point<dim>(-10, -10), 2);
+    }
+  else if (dim == 3)
+    {
+      // Radius set such that every region has 2^3 = 8 cells
+      predicates[1] = predicate_template<dim>(Point<dim>(7.5, 7.5, 7.5), 3);
+      predicates[2] = predicate_template<dim>(Point<dim>(5, 5, 5), 3);
+      predicates[0] = predicate_template<dim>(Point<dim>(), 3);
+      predicates[3] = predicate_template<dim>(Point<dim>(-5, -5, -5), 3);
+      predicates[4] = predicate_template<dim>(Point<dim>(-10, -10, -10), 3);
+    }
+
+  // Check pair-wise connections between predicate regions
+  for (int i = 0; i < 5; ++i)
+    for (int j = 0; j < 5; ++j)
+      {
+        deallog << i << ":" << j << "="
+                << ColorEnriched::internal::find_connection_between_subdomains(
+                     dof_handler, predicates[i], predicates[j])
+                << std::endl;
+      }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  try
+    {
+      test<2>();
+      test<3>();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_01.output b/tests/fe/fe_enriched_color_01.output
new file mode 100644 (file)
index 0000000..919b93e
--- /dev/null
@@ -0,0 +1,53 @@
+
+DEAL::dim = 2
+DEAL::0:0=1
+DEAL::0:1=0
+DEAL::0:2=1
+DEAL::0:3=1
+DEAL::0:4=0
+DEAL::1:0=0
+DEAL::1:1=1
+DEAL::1:2=1
+DEAL::1:3=0
+DEAL::1:4=0
+DEAL::2:0=1
+DEAL::2:1=1
+DEAL::2:2=1
+DEAL::2:3=0
+DEAL::2:4=0
+DEAL::3:0=1
+DEAL::3:1=0
+DEAL::3:2=0
+DEAL::3:3=1
+DEAL::3:4=1
+DEAL::4:0=0
+DEAL::4:1=0
+DEAL::4:2=0
+DEAL::4:3=1
+DEAL::4:4=1
+DEAL::dim = 3
+DEAL::0:0=1
+DEAL::0:1=0
+DEAL::0:2=1
+DEAL::0:3=1
+DEAL::0:4=0
+DEAL::1:0=0
+DEAL::1:1=1
+DEAL::1:2=1
+DEAL::1:3=0
+DEAL::1:4=0
+DEAL::2:0=1
+DEAL::2:1=1
+DEAL::2:2=1
+DEAL::2:3=0
+DEAL::2:4=0
+DEAL::3:0=1
+DEAL::3:1=0
+DEAL::3:2=0
+DEAL::3:3=1
+DEAL::3:4=1
+DEAL::4:0=0
+DEAL::4:1=0
+DEAL::4:2=0
+DEAL::4:3=1
+DEAL::4:4=1
diff --git a/tests/fe/fe_enriched_color_02.cc b/tests/fe/fe_enriched_color_02.cc
new file mode 100644 (file)
index 0000000..f0657fd
--- /dev/null
@@ -0,0 +1,163 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the function ColorEnriched::internal::color_predicates.
+ * Check if the function correctly colors vector of predicates
+ * using graph coloring.
+ *
+ * Two predicates are said to be connected if cells belonging to
+ * different predicates touch each other.
+ */
+
+#include <deal.II/fe/fe_enriched.templates.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include "../tests.h"
+
+/*
+ * Predicate function needed by ColorEnriched::internal::color_predicates
+ * implemented using a struct.
+ */
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+
+/*
+ * Type used to defined vector of predicates needed by the function
+ * ColorEnriched::internal::color_predicates.
+ */
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+int
+main(int argc, char **argv)
+{
+  // Intialize MPI as required by Zoltan library used for graph coloring by this
+  // test.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+  MPILogInitAll                    all;
+
+  // Make basic grid
+  const unsigned int dim = 2;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -20, 20);
+  triangulation.refine_global(4);
+  hp::DoFHandler<dim> dof_handler(triangulation);
+
+  // check the coloring function on different set of predicates.
+  std::vector<predicate_function<dim>> vec_predicates;
+  std::vector<unsigned int>            predicate_colors;
+  {
+    // case 1: predicates are not connected
+    vec_predicates.clear();
+    vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-10, 10), 2));
+    vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 0), 2));
+
+    predicate_colors.resize(vec_predicates.size());
+
+    ColorEnriched::internal::color_predicates<dim>(
+      dof_handler, vec_predicates, predicate_colors);
+
+    deallog << "Case 1" << std::endl;
+    for (auto i : predicate_colors)
+      {
+        deallog << i << std::endl;
+      }
+  }
+
+  {
+    // case 2: Two predicates that are connected.
+    vec_predicates.clear();
+    vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-10, 10), 2));
+    vec_predicates.push_back(
+      EnrichmentPredicate<dim>(Point<dim>(-7.5, 7.5), 2));
+
+    predicate_colors.resize(vec_predicates.size());
+
+    ColorEnriched::internal::color_predicates(
+      dof_handler, vec_predicates, predicate_colors);
+
+    deallog << "Case 2" << std::endl;
+    for (auto i : predicate_colors)
+      {
+        deallog << i << std::endl;
+      }
+  }
+
+  {
+    // case 3: connections between predicates is as follows -
+    // 0-1 (overlap connection),
+    // 3-4 (edge connection)
+    vec_predicates.clear();
+    vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-10, 10), 2));
+    vec_predicates.push_back(
+      EnrichmentPredicate<dim>(Point<dim>(-7.5, 7.5), 2));
+    vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 0), 2));
+    vec_predicates.push_back(
+      EnrichmentPredicate<dim>(Point<dim>(7.5, -7.5), 2));
+    vec_predicates.push_back(
+      EnrichmentPredicate<dim>(Point<dim>(12.5, -12.5), 2));
+
+    predicate_colors.resize(vec_predicates.size());
+
+    ColorEnriched::internal::color_predicates(
+      dof_handler, vec_predicates, predicate_colors);
+
+    deallog << "Case 3" << std::endl;
+    for (auto i : predicate_colors)
+      {
+        deallog << i << std::endl;
+      }
+  }
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_02.with_mpi=true.output b/tests/fe/fe_enriched_color_02.with_mpi=true.output
new file mode 100644 (file)
index 0000000..64cb19d
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:0::Case 1
+DEAL:0::1
+DEAL:0::1
+DEAL:0::Case 2
+DEAL:0::1
+DEAL:0::2
+DEAL:0::Case 3
+DEAL:0::1
+DEAL:0::2
+DEAL:0::1
+DEAL:0::1
+DEAL:0::2
diff --git a/tests/fe/fe_enriched_color_03.cc b/tests/fe/fe_enriched_color_03.cc
new file mode 100644 (file)
index 0000000..24b69ad
--- /dev/null
@@ -0,0 +1,175 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test function - ColorEnriched::internal::set_cellwise_color_set_and_fe_index
+ * for a set of predicates.
+ * Check for each cell, if appropriate fe index and color-index map is set.
+ * Color-index map associates different colors of different enrichment
+ * functions with corresponding enrichment function index.
+ */
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_enriched.templates.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include <map>
+
+#include "../tests.h"
+
+/*
+ * Predicate function needed by ColorEnriched::internal::color_predicates
+ * implemented using a struct.
+ */
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+
+/*
+ * Type used to defined vector of predicates needed by the function
+ * ColorEnriched::internal::color_predicates.
+ */
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+int
+main(int argc, char **argv)
+{
+  // Intialize MPI as required by Zoltan library used for graph coloring by this
+  // test.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+  MPILogInitAll                    all;
+
+  // Make basic grid
+  const unsigned int  dim = 2;
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  GridGenerator::hyper_cube(triangulation, -2, 2);
+  triangulation.refine_global(2);
+
+  // Make predicates. Predicate 0 and 1 overlap.
+  // Predicate 2 is connected to 0.
+  std::vector<predicate_function<dim>> vec_predicates;
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-1, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(1.5, -1.5), 1));
+
+  // Do manual coloring since we are not testing coloring function here!
+  std::vector<unsigned int> predicate_colors;
+  predicate_colors.resize(vec_predicates.size());
+  ColorEnriched::internal::color_predicates(
+    dof_handler, vec_predicates, predicate_colors);
+
+  // Make required objects to call function set_cellwise_color_set_and_fe_index
+  std::map<unsigned int, std::map<unsigned int, unsigned int>>
+                                      cellwise_color_predicate_map;
+  std::vector<std::set<unsigned int>> fe_sets;
+  ColorEnriched::internal::set_cellwise_color_set_and_fe_index(
+    dof_handler,
+    vec_predicates,
+    predicate_colors,
+    cellwise_color_predicate_map,
+    fe_sets);
+
+  /*
+   * Run through active cells to check fe index, colors of
+   * enrichment functions associated with it.
+   *
+   * A unique color set corresponds to an fe index.
+   *
+   * Eg: If an fe index 1 corresponds to color set {2,3},
+   * means that a cell with fe index 1 has enrichment functions
+   * which are colored 2 and 3. Here different enrichment function
+   * have same color 2 but for a given cell only one of them would
+   * be relevant. So all additional information we need is which
+   * enrichment function is relevant that has color 2. We need to
+   * do the same thing with color 2 as well.
+   *
+   * Each cell is assigned unique material id by the function
+   * set_cellwise_color_set_and_fe_index. Now using material id,
+   * each cell is associated with a map which assigns a color to a
+   * particular enrichment function id.
+   *
+   */
+  auto cell = dof_handler.begin_active();
+  auto endc = dof_handler.end();
+  for (unsigned int cell_index = 0; cell != endc; ++cell, ++cell_index)
+    {
+      // print true predicates for a cell as a binary code
+      // 0 1 0 indicates predicate with index 2 is true in the cell.
+      deallog << cell->index() << ":predicates=";
+      for (auto predicate : vec_predicates)
+        deallog << predicate(cell) << ":";
+
+      // print color predicate pairs for a cell
+      // Note that here material id is used to identify cells
+      // Here (1,2) indicates predicate 2 of color 1 is relavant for cell.
+      deallog << "(color, enrichment_func_id):";
+      for (auto color_predicate_pair :
+           cellwise_color_predicate_map[cell->material_id()])
+        {
+          deallog << "(" << color_predicate_pair.first << ","
+                  << color_predicate_pair.second << "):";
+        }
+
+      // For a cell, print fe active index and corresponding fe set.
+      //{1,2} indicates 2 enrichment functions of color 1 and 2 are relevant.
+      deallog << ":fe_active_index:" << cell->active_fe_index() << ":fe_set:";
+      for (auto fe_set_element : fe_sets[cell->active_fe_index()])
+        deallog << fe_set_element << ":";
+      deallog << std::endl;
+    }
+
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_03.with_mpi=true.output b/tests/fe/fe_enriched_color_03.with_mpi=true.output
new file mode 100644 (file)
index 0000000..38128ac
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL:0::0:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::1:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::2:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::3:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::4:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::5:predicates=0:0:1:(color, enrichment_func_id):(1,2)::fe_active_index:1:fe_set:1:
+DEAL:0::6:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::7:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::8:predicates=0:1:0:(color, enrichment_func_id):(2,1)::fe_active_index:2:fe_set:2:
+DEAL:0::9:predicates=1:1:0:(color, enrichment_func_id):(1,0):(2,1)::fe_active_index:3:fe_set:1:2:
+DEAL:0::10:predicates=0:1:0:(color, enrichment_func_id):(2,1)::fe_active_index:2:fe_set:2:
+DEAL:0::11:predicates=1:1:0:(color, enrichment_func_id):(1,0):(2,1)::fe_active_index:3:fe_set:1:2:
+DEAL:0::12:predicates=1:0:0:(color, enrichment_func_id):(1,0)::fe_active_index:1:fe_set:1:
+DEAL:0::13:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
+DEAL:0::14:predicates=1:0:0:(color, enrichment_func_id):(1,0)::fe_active_index:1:fe_set:1:
+DEAL:0::15:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set:
diff --git a/tests/fe/fe_enriched_color_04.cc b/tests/fe/fe_enriched_color_04.cc
new file mode 100644 (file)
index 0000000..ffd8c5a
--- /dev/null
@@ -0,0 +1,176 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test function: ColorEnriched::internal::make_colorwise_enrichment_functions
+ * for a set of predicates.
+ *
+ * The container of enrichment functions created by the function is
+ * called color_enrichments because they return the correct enrichment
+ * function if cell and color of the function are provided.
+ */
+
+#include <deal.II/fe/fe_enriched.templates.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include <map>
+
+#include "../tests.h"
+
+/*
+ * Predicate function needed by ColorEnriched::internal::color_predicates
+ * implemented using a struct.
+ */
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+
+/*
+ * Type used to defined vector of predicates needed by the function
+ * ColorEnriched::internal::color_predicates.
+ */
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+int
+main(int argc, char **argv)
+{
+  // Intialize MPI as required by Zoltan library used for graph coloring by this
+  // test.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+  MPILogInitAll                    all;
+
+  // Make basic grid
+  const unsigned int  dim = 2;
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  GridGenerator::hyper_cube(triangulation, -2, 2);
+  triangulation.refine_global(2);
+
+
+  // Make predicates. Predicate 0 and 1 overlap.
+  // Predicate 2 is connected to 0.
+  std::vector<predicate_function<dim>> vec_predicates;
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-1, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(1.5, -1.5), 1));
+
+  // Do manual coloring since we are not testing coloring function here!
+  std::vector<unsigned int> predicate_colors;
+  predicate_colors.resize(vec_predicates.size());
+  unsigned int num_colors = ColorEnriched::internal::color_predicates(
+    dof_handler, vec_predicates, predicate_colors);
+
+
+  // Make required objects to call function set_cellwise_color_set_and_fe_index
+  std::map<unsigned int, std::map<unsigned int, unsigned int>>
+                                      cellwise_color_predicate_map;
+  std::vector<std::set<unsigned int>> fe_sets;
+  ColorEnriched::internal::set_cellwise_color_set_and_fe_index(
+    dof_handler,
+    vec_predicates,
+    predicate_colors,
+    cellwise_color_predicate_map,
+    fe_sets);
+
+  // Construct vector of enrichment functions
+  std::vector<std::shared_ptr<Function<dim>>> vec_enrichments;
+  vec_enrichments.reserve(vec_predicates.size());
+  for (unsigned int i = 0; i < vec_predicates.size(); ++i)
+    {
+      // constant function is chosen as enrichment function
+      ConstantFunction<dim> func(i);
+      vec_enrichments.push_back(std::make_shared<ConstantFunction<dim>>(func));
+    }
+
+  // Construct container for color enrichment functions needed
+  // by function make_colorwise_enrichment_functions
+  std::vector<std::function<const Function<dim> *(
+    const typename Triangulation<dim>::cell_iterator &)>>
+    color_enrichments;
+
+  ColorEnriched::internal::make_colorwise_enrichment_functions<dim, dim>(
+    num_colors,
+    vec_enrichments,
+    cellwise_color_predicate_map,
+    color_enrichments);
+
+
+  deallog << "color wise enrichment functions:" << std::endl;
+  auto cell = dof_handler.begin_active();
+  auto endc = dof_handler.end();
+  for (unsigned int cell_index = 0; cell != endc; ++cell, ++cell_index)
+    {
+      // Print ids of predicates active in the cell
+      unsigned int cell_id = cell->index();
+      deallog << cell_id << ":predicates=";
+      for (auto predicate : vec_predicates)
+        deallog << predicate(cell) << ":";
+
+      /*
+       * Check if a color and enrichment index map exists for the cell.
+       * If so print the color and corresponding enrichment functions
+       * value at the cell center. Note that indices of color enrichment
+       * functions starts with zero not one unlike color indices.
+       */
+      if (cellwise_color_predicate_map.count(cell->material_id()) == 1)
+        for (unsigned int color = 1; color <= num_colors; ++color)
+          if (cellwise_color_predicate_map.at(cell->material_id())
+                .count(color) == 1)
+            deallog << ":color:" << color << ":func_value:"
+                    << color_enrichments[color - 1](cell)->value(
+                         cell->center());
+
+      deallog << std::endl;
+    }
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_04.with_mpi=true.output b/tests/fe/fe_enriched_color_04.with_mpi=true.output
new file mode 100644 (file)
index 0000000..96671dd
--- /dev/null
@@ -0,0 +1,18 @@
+
+DEAL:0::color wise enrichment functions:
+DEAL:0::0:predicates=0:0:0:
+DEAL:0::1:predicates=0:0:0:
+DEAL:0::2:predicates=0:0:0:
+DEAL:0::3:predicates=0:0:0:
+DEAL:0::4:predicates=0:0:0:
+DEAL:0::5:predicates=0:0:1::color:1:func_value:2.00000
+DEAL:0::6:predicates=0:0:0:
+DEAL:0::7:predicates=0:0:0:
+DEAL:0::8:predicates=0:1:0::color:2:func_value:1.00000
+DEAL:0::9:predicates=1:1:0::color:1:func_value:0.00000:color:2:func_value:1.00000
+DEAL:0::10:predicates=0:1:0::color:2:func_value:1.00000
+DEAL:0::11:predicates=1:1:0::color:1:func_value:0.00000:color:2:func_value:1.00000
+DEAL:0::12:predicates=1:0:0::color:1:func_value:0.00000
+DEAL:0::13:predicates=0:0:0:
+DEAL:0::14:predicates=1:0:0::color:1:func_value:0.00000
+DEAL:0::15:predicates=0:0:0:
diff --git a/tests/fe/fe_enriched_color_05.cc b/tests/fe/fe_enriched_color_05.cc
new file mode 100644 (file)
index 0000000..339e4a3
--- /dev/null
@@ -0,0 +1,264 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test function: ColorEnriched::internal
+ * ::make_fe_collection_from_colored_enrichments for a set of predicates.
+ *
+ * The function return FE_Collection which is then printed to test.
+ */
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_enriched.templates.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/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/data_postprocessor.h>
+
+#include <map>
+
+#include "../tests.h"
+
+// uncomment when debugging
+//#define DATA_OUT_FE_ENRICHED
+
+/*
+ * Predicate function needed by ColorEnriched::internal::color_predicates
+ * implemented using a struct.
+ */
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+
+/*
+ * Type used to defined vector of predicates needed by the function
+ * ColorEnriched::internal::color_predicates.
+ */
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+template <int dim>
+void
+plot_shape_function(hp::DoFHandler<dim> &dof_handler, unsigned int patches = 5)
+{
+  std::cout << "n_cells: " << dof_handler.get_triangulation().n_active_cells()
+            << std::endl;
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  constraints.close();
+
+  // output to check if all is good:
+  std::vector<Vector<double>> shape_functions;
+  std::vector<std::string>    names;
+  for (unsigned int s = 0; s < dof_handler.n_dofs(); s++)
+    {
+      Vector<double> shape_function;
+      shape_function.reinit(dof_handler.n_dofs());
+      shape_function[s] = 1.0;
+
+      // if the dof is constrained, first output unconstrained vector
+      if (constraints.is_constrained(s))
+        {
+          names.push_back(std::string("UN_") +
+                          dealii::Utilities::int_to_string(s, 2));
+          shape_functions.push_back(shape_function);
+        }
+
+      names.push_back(std::string("N_") +
+                      dealii::Utilities::int_to_string(s, 2));
+
+      // make continuous/constrain:
+      constraints.distribute(shape_function);
+      shape_functions.push_back(shape_function);
+    }
+
+  DataOut<dim, hp::DoFHandler<dim>> data_out;
+  data_out.attach_dof_handler(dof_handler);
+
+  // get material ids:
+  Vector<float> fe_index(dof_handler.get_triangulation().n_active_cells());
+  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler
+                                                              .begin_active(),
+                                                     endc = dof_handler.end();
+  for (unsigned int index = 0; cell != endc; ++cell, ++index)
+    {
+      fe_index[index] = cell->active_fe_index();
+    }
+  data_out.add_data_vector(fe_index, "fe_index");
+
+  for (unsigned int i = 0; i < shape_functions.size(); i++)
+    data_out.add_data_vector(shape_functions[i], names[i]);
+
+  data_out.build_patches(patches);
+
+  std::string filename =
+    "hp-shape_functions_" + dealii::Utilities::int_to_string(dim) + "D.vtu";
+  std::ofstream output(filename.c_str());
+  data_out.write_vtu(output);
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  // Intialize MPI as required by Zoltan library used for graph coloring by this
+  // test.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+  MPILogInitAll                    all;
+
+  // Make basic grid
+  const unsigned int  dim = 2;
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  GridGenerator::hyper_cube(triangulation, -2, 2);
+  triangulation.refine_global(2);
+
+  // Make predicates. Predicate 0 and 1 overlap.
+  std::vector<predicate_function<dim>> vec_predicates;
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-1, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 1), 1));
+
+  // find colors for predicates
+  std::vector<unsigned int> predicate_colors;
+  predicate_colors.resize(vec_predicates.size());
+  unsigned int num_colors = ColorEnriched::internal::color_predicates(
+    dof_handler, vec_predicates, predicate_colors);
+
+  // Make required objects to call function set_cellwise_color_set_and_fe_index
+  std::map<unsigned int, std::map<unsigned int, unsigned int>>
+                                      cellwise_color_predicate_map;
+  std::vector<std::set<unsigned int>> fe_sets;
+
+  ColorEnriched::internal::set_cellwise_color_set_and_fe_index(
+    dof_handler,
+    vec_predicates,
+    predicate_colors,
+    cellwise_color_predicate_map,
+    fe_sets);
+
+  // Construct vector of enrichment functions
+  std::vector<std::shared_ptr<Function<dim>>> vec_enrichments;
+  vec_enrichments.reserve(vec_predicates.size());
+  for (unsigned int i = 0; i < vec_predicates.size(); ++i)
+    {
+      // constant function.
+      ConstantFunction<dim> func(10 + i); // constant function
+      vec_enrichments.push_back(std::make_shared<ConstantFunction<dim>>(func));
+    }
+
+  // Construct container for color enrichment functions needed
+  // by function make_colorwise_enrichment_functions
+  std::vector<std::function<const Function<dim> *(
+    const typename Triangulation<dim>::cell_iterator &)>>
+    color_enrichments;
+
+  ColorEnriched::internal::make_colorwise_enrichment_functions<dim, dim>(
+    num_colors,      // needs number of colors
+    vec_enrichments, // enrichment functions based on predicate id
+    cellwise_color_predicate_map,
+    color_enrichments);
+
+  // Construct object needed to call make_fe_collection_from_colored_enrichments
+  FE_Q<dim>             fe_base(2);
+  FE_Q<dim>             fe_enriched(1);
+  FE_Nothing<dim>       fe_nothing(1, true);
+  hp::FECollection<dim> fe_collection;
+  ColorEnriched::internal::make_fe_collection_from_colored_enrichments(
+    num_colors,
+    fe_sets,           // total list of color sets possible
+    color_enrichments, // color wise enrichment functions
+    fe_base,           // basic fe element
+    fe_enriched,       // fe element multiplied by enrichment function
+    fe_nothing,
+    fe_collection);
+
+  // print all the different fe sets needed by different cells
+  deallog << "fe sets:" << std::endl;
+  for (auto fe_set : fe_sets)
+    {
+      deallog << "color:";
+      for (auto color : fe_set)
+        deallog << ":" << color;
+      deallog << std::endl;
+    }
+
+  // check if fe_collection is correctly constructed by function
+  deallog << "fe_collection[index] mapping:" << std::endl;
+  for (unsigned int index = 0; index != fe_collection.size(); ++index)
+    {
+      deallog << "name:" << fe_collection[index].get_name() << std::endl;
+      deallog << "n_blocks:" << fe_collection[index].n_blocks() << std::endl;
+      deallog << "n_comp:" << fe_collection[index].n_components() << std::endl;
+      deallog << "n_dofs:" << fe_collection[index].n_dofs_per_cell()
+              << std::endl;
+    }
+
+#ifdef DATA_OUT_FE_ENRICHED
+  GridTools::partition_triangulation(
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD), triangulation);
+  dof_handler.distribute_dofs(*fe_collection);
+
+  plot_shape_function<dim>(dof_handler, 5);
+#endif
+
+  dof_handler.clear();
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_05.with_mpi=true.output b/tests/fe/fe_enriched_color_05.with_mpi=true.output
new file mode 100644 (file)
index 0000000..6258268
--- /dev/null
@@ -0,0 +1,23 @@
+
+DEAL:0::fe sets:
+DEAL:0::color:
+DEAL:0::color::1
+DEAL:0::color::1:2
+DEAL:0::color::2
+DEAL:0::fe_collection[index] mapping:
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:9
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:13
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:17
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:13
diff --git a/tests/fe/fe_enriched_color_06.cc b/tests/fe/fe_enriched_color_06.cc
new file mode 100644 (file)
index 0000000..3b48a1b
--- /dev/null
@@ -0,0 +1,142 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test function: ColorEnriched::internal
+ * ::make_fe_collection_from_colored_enrichments for a set of predicates.
+ *
+ * The function return FE_Collection which is then printed to test.
+ */
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_enriched.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/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/data_postprocessor.h>
+
+#include <map>
+
+#include "../tests.h"
+
+/*
+ * Testing helper class
+ */
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+
+/*
+ * Type used to defined vector of predicates needed by the function
+ * ColorEnriched::internal::color_predicates.
+ */
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+int
+main(int argc, char **argv)
+{
+  // Intialize MPI as required by Zoltan library used for graph coloring by this
+  // test.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+  MPILogInitAll                    all;
+
+  // Make basic grid
+  const unsigned int  dim = 2;
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  GridGenerator::hyper_cube(triangulation, -2, 2);
+  triangulation.refine_global(2);
+
+  // Make predicates. Predicate 0 and 1 overlap.
+  std::vector<predicate_function<dim>> vec_predicates;
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(-1, 1), 1));
+  vec_predicates.push_back(EnrichmentPredicate<dim>(Point<dim>(0, 1), 1));
+
+  // Construct vector of enrichment functions
+  std::vector<std::shared_ptr<Function<dim>>> vec_enrichments;
+  vec_enrichments.reserve(vec_predicates.size());
+  for (unsigned int i = 0; i < vec_predicates.size(); ++i)
+    {
+      // constant function.
+      ConstantFunction<dim> func(10 + i); // constant function
+      vec_enrichments.push_back(std::make_shared<ConstantFunction<dim>>(func));
+    }
+
+  // Construct helper class to construct fe collection
+  FE_Q<dim>                         fe_base(2);
+  FE_Q<dim>                         fe_enriched(1);
+  static ColorEnriched::Helper<dim> fe_space(
+    fe_base, fe_enriched, vec_predicates, vec_enrichments);
+  const hp::FECollection<dim> &fe_collection(
+    fe_space.build_fe_collection(dof_handler));
+
+  // check if fe_collection is correctly constructed by function
+  deallog << "fe_collection[index] mapping:" << std::endl;
+  for (unsigned int index = 0; index != fe_collection.size(); ++index)
+    {
+      deallog << "name:" << fe_collection[index].get_name() << std::endl;
+      deallog << "n_blocks:" << fe_collection[index].n_blocks() << std::endl;
+      deallog << "n_comp:" << fe_collection[index].n_components() << std::endl;
+      deallog << "n_dofs:" << fe_collection[index].n_dofs_per_cell()
+              << std::endl;
+    }
+
+  dof_handler.clear();
+  return 0;
+}
diff --git a/tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output b/tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output
new file mode 100644 (file)
index 0000000..bcfe1f9
--- /dev/null
@@ -0,0 +1,18 @@
+
+DEAL:0::fe_collection[index] mapping:
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:9
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:13
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:17
+DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)]
+DEAL:0::n_blocks:3
+DEAL:0::n_comp:1
+DEAL:0::n_dofs:13
diff --git a/tests/fe/fe_enriched_color_07.cc b/tests/fe/fe_enriched_color_07.cc
new file mode 100644 (file)
index 0000000..083c61c
--- /dev/null
@@ -0,0 +1,1988 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Testing the ColorEnriched::Helper class by solving problems
+ */
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_cspline.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/parsed_function.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_enriched.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_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/slepc_solver.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <math.h>
+
+#include <map>
+#include <set>
+#include <vector>
+
+#include "../tests.h"
+
+template <int dim>
+class SigmaFunction : public Function<dim>
+{
+  Point<dim>          center;
+  FunctionParser<dim> func;
+
+public:
+  SigmaFunction() : Function<dim>(), func(1)
+  {}
+
+  // to help with resize function. doesn't copy function parser(func)!
+  SigmaFunction(SigmaFunction &&other) : center(other.center), func(1)
+  {}
+
+  void
+  initialize(const Point<dim> & center,
+             const double &     sigma,
+             const std::string &func_expr);
+  double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+  Tensor<1, dim>
+  gradient(const Point<dim> &p, const unsigned int component = 0) const;
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          value_list) const;
+};
+
+template <int dim>
+void
+SigmaFunction<dim>::initialize(const Point<dim> & _center,
+                               const double &     sigma,
+                               const std::string &func_expr)
+{
+  center = _center;
+  std::string                   variables;
+  std::map<std::string, double> constants = {{"sigma", sigma},
+                                             {"pi", numbers::PI}};
+
+  AssertThrow(dim == 1 || dim == 2 || dim == 3,
+              ExcMessage("Dimension not implemented"));
+  switch (dim)
+    {
+      case 1:
+        variables = "x";
+        break;
+      case 2:
+        variables = "x, y";
+        break;
+      case 3:
+        variables = "x, y, z";
+        break;
+    }
+
+  func.initialize(variables, func_expr, constants);
+}
+
+template <int dim>
+inline double
+SigmaFunction<dim>::value(const Point<dim> & p,
+                          const unsigned int component) const
+{
+  const Point<dim> d(p - center);
+  return func.value(d, component);
+}
+
+
+template <int dim>
+inline Tensor<1, dim>
+SigmaFunction<dim>::gradient(const Point<dim> & p,
+                             const unsigned int component) const
+{
+  const Point<dim> d(p - center);
+  return func.gradient(d, component);
+}
+
+template <int dim>
+void
+SigmaFunction<dim>::value_list(const std::vector<Point<dim>> &points,
+                               std::vector<double> &          value_list) const
+{
+  const unsigned int n_points = points.size();
+
+  AssertDimension(points.size(), value_list.size());
+
+  for (unsigned int p = 0; p < n_points; ++p)
+    value_list[p] = value(points[p]);
+}
+
+
+
+template <int dim>
+struct EnrichmentPredicate
+{
+  EnrichmentPredicate(const Point<dim> origin, const double radius) :
+    origin(origin),
+    radius(radius)
+  {}
+
+  template <class Iterator>
+  bool
+  operator()(const Iterator &i) const
+  {
+    return ((i->center() - origin).norm_square() < radius * radius);
+  }
+
+  const Point<dim> &
+  get_origin()
+  {
+    return origin;
+  }
+  const double &
+  get_radius()
+  {
+    return radius;
+  }
+
+private:
+  const Point<dim> origin;
+  const double     radius;
+};
+
+
+template <int dim>
+class SplineEnrichmentFunction : public Function<dim>
+{
+public:
+  SplineEnrichmentFunction(const Point<dim> &         origin,
+                           const std::vector<double> &interpolation_points_1d,
+                           const std::vector<double> &interpolation_values_1d) :
+    Function<dim>(1),
+    origin(origin),
+    interpolation_points(interpolation_points_1d),
+    interpolation_values(interpolation_values_1d),
+    cspline(interpolation_points, interpolation_values)
+  {}
+
+  SplineEnrichmentFunction(SplineEnrichmentFunction &&other) :
+    Function<dim>(1),
+    origin(other.origin),
+    interpolation_points(other.interpolation_points),
+    interpolation_values(other.interpolation_values),
+    cspline(interpolation_points, interpolation_values)
+  {}
+
+  SplineEnrichmentFunction(const SplineEnrichmentFunction &other) :
+    Function<dim>(1),
+    origin(other.origin),
+    interpolation_points(other.interpolation_points),
+    interpolation_values(other.interpolation_values),
+    cspline(interpolation_points, interpolation_values)
+  {}
+
+
+
+  virtual double
+  value(const Point<dim> &point, const unsigned int component = 0) const
+  {
+    Tensor<1, dim> dist = point - origin;
+    const double   r    = dist.norm();
+    return cspline.value(Point<1>(r), component);
+  }
+
+  virtual Tensor<1, dim>
+  gradient(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    Tensor<1, dim> dist = p - origin;
+    const double   r    = dist.norm();
+    Assert(r > 0., ExcDivideByZero());
+    dist /= r;
+    Assert(component == 0, ExcMessage("Not implemented"));
+    return cspline.gradient(Point<1>(r))[0] * dist;
+  }
+
+private:
+  /**
+   * origin
+   */
+  const Point<dim>    origin;
+  std::vector<double> interpolation_points;
+  std::vector<double> interpolation_values;
+  // enrichment function as CSpline based on radius
+  Functions::CSpline<1> cspline;
+};
+
+
+
+struct ParameterCollection
+{
+  ParameterCollection(const std::string &file_name);
+
+  ParameterCollection(const int &                dim,
+                      const double &             size,
+                      const unsigned int &       shape,
+                      const unsigned int &       global_refinement,
+                      const unsigned int &       cycles,
+                      const unsigned int &       fe_base_degree,
+                      const unsigned int &       fe_enriched_degree,
+                      const unsigned int &       max_iterations,
+                      const double &             tolerance,
+                      const std::string &        rhs_value_expr,
+                      const std::string &        boundary_value_expr,
+                      const std::string &        rhs_radial_problem,
+                      const std::string &        boundary_radial_problem,
+                      const std::string &        exact_soln_expr,
+                      const unsigned int &       patches,
+                      const unsigned int &       debug_level,
+                      const unsigned int &       n_enrichments,
+                      const std::vector<double> &points_enrichments,
+                      const std::vector<double> &radii_predicates,
+                      const std::vector<double> &sigmas);
+
+  void
+  print();
+
+  void set_enrichment_point(Point<2> &p, const unsigned int i)
+  {
+    AssertDimension(dim, 2);
+    p(0) = points_enrichments[2 * i];
+    p(1) = points_enrichments[2 * i + 1];
+  }
+  void set_enrichment_point(Point<3> &p, const unsigned int i)
+  {
+    AssertDimension(dim, 3);
+    p(0) = points_enrichments[3 * i];
+    p(1) = points_enrichments[3 * i + 1];
+    p(2) = points_enrichments[3 * i + 2];
+  }
+
+  int          dim;
+  double       size;
+  unsigned int shape; // 0 = ball, 1 = cube
+  unsigned int global_refinement;
+  unsigned int cycles;
+  unsigned int fe_base_degree;
+  unsigned int fe_enriched_degree;
+  unsigned int max_iterations;
+  double       tolerance;
+
+  // parameters related to exact solution
+  std::string rhs_value_expr;
+  std::string boundary_value_expr;
+
+  // value = true ==> estimate exact solution from radial problem
+  std::string rhs_radial_problem;
+  std::string boundary_radial_problem;
+
+  std::string exact_soln_expr;
+
+  unsigned int patches;
+  // debug level = 0(output nothing),
+  // 1 (print statements)
+  // 2 (output solution)
+  // 3 (+ output grid data as well)
+  // 9 (+ shape functions as well)
+  unsigned int        debug_level;
+  unsigned int        n_enrichments;
+  std::vector<double> points_enrichments;
+  std::vector<double> radii_predicates;
+  std::vector<double> sigmas;
+};
+
+
+
+ParameterCollection::ParameterCollection(const std::string &file_name)
+{
+  // std::cout << "...reading parameters" << std::endl;
+
+  ParameterHandler prm;
+
+  // declare parameters
+  prm.enter_subsection("geometry");
+  prm.declare_entry("dim", "2", Patterns::Integer());
+  prm.declare_entry("size", "1", Patterns::Double(0));
+  prm.declare_entry("shape", "1", Patterns::Integer(0));
+  prm.declare_entry("Global refinement", "1", Patterns::Integer(1));
+  prm.declare_entry("cycles", "0", Patterns::Integer(0));
+  prm.leave_subsection();
+
+
+  prm.enter_subsection("solver");
+  prm.declare_entry("fe base degree", "1", Patterns::Integer(1));
+  prm.declare_entry("fe enriched degree", "1", Patterns::Integer(1));
+  prm.declare_entry("max iterations", "1000", Patterns::Integer(1));
+  prm.declare_entry("tolerance", "1e-8", Patterns::Double(0));
+  prm.leave_subsection();
+
+
+  prm.enter_subsection("expressions");
+  prm.declare_entry("rhs value", "0", Patterns::Anything());
+  prm.declare_entry("boundary value", "0", Patterns::Anything());
+  prm.declare_entry("rhs value radial problem", "0", Patterns::Anything());
+  prm.declare_entry("boundary value radial problem", "0", Patterns::Anything());
+  prm.declare_entry("exact solution expression", "", Patterns::Anything());
+  prm.declare_entry("estimate exact solution", "false", Patterns::Bool());
+  prm.leave_subsection();
+
+
+  prm.enter_subsection("output");
+  prm.declare_entry("patches", "1", Patterns::Integer(1));
+  prm.declare_entry("debug level", "0", Patterns::Integer(0, 9));
+  prm.leave_subsection();
+
+  // parse parameter file
+  prm.parse_input(file_name, "#end-of-dealii parser");
+
+
+  // get parameters
+  prm.enter_subsection("geometry");
+  dim               = prm.get_integer("dim");
+  size              = prm.get_double("size");
+  shape             = prm.get_integer("shape");
+  global_refinement = prm.get_integer("Global refinement");
+  cycles            = prm.get_integer("cycles");
+  prm.leave_subsection();
+
+  prm.enter_subsection("solver");
+  fe_base_degree     = prm.get_integer("fe base degree");
+  fe_enriched_degree = prm.get_integer("fe enriched degree");
+  max_iterations     = prm.get_integer("max iterations");
+  tolerance          = prm.get_double("tolerance");
+  prm.leave_subsection();
+
+  prm.enter_subsection("expressions");
+  rhs_value_expr          = prm.get("rhs value");
+  boundary_value_expr     = prm.get("boundary value");
+  rhs_radial_problem      = prm.get("rhs value radial problem");
+  boundary_radial_problem = prm.get("boundary value radial problem");
+  exact_soln_expr         = prm.get("exact solution expression");
+  prm.leave_subsection();
+
+  prm.enter_subsection("output");
+  patches     = prm.get_integer("patches");
+  debug_level = prm.get_integer("debug level");
+  prm.leave_subsection();
+
+
+  // manual parsing
+  // open parameter file
+  std::ifstream prm_file(file_name);
+
+  // read lines until "#end-of-dealii parser" is reached
+  std::string line;
+  while (getline(prm_file, line))
+    if (line == "#end-of-dealii parser")
+      break;
+
+  AssertThrow(
+    line == "#end-of-dealii parser",
+    ExcMessage("line missing in parameter file = \'#end-of-dealii parser\' "));
+
+  // function to read next line not starting with # or empty
+  auto read_next_proper_line = [&](std::string &line) {
+    while (getline(prm_file, line))
+      {
+        if (line.size() == 0 || line[0] == '#' || line[0] == ' ')
+          continue;
+        else
+          break;
+      }
+  };
+
+  std::stringstream s_stream;
+
+  // read num of enrichement points
+  read_next_proper_line(line);
+  s_stream.str(line);
+  s_stream >> n_enrichments;
+
+  // note vector of points
+  for (unsigned int i = 0; i != n_enrichments; ++i)
+    {
+      read_next_proper_line(line);
+      s_stream.clear();
+      s_stream.str(line);
+
+      points_enrichments.resize(dim * n_enrichments);
+
+      if (dim == 2)
+        {
+          double x, y;
+          s_stream >> x >> y;
+          points_enrichments[2 * i]     = x;
+          points_enrichments[2 * i + 1] = y;
+        }
+      else if (dim == 3)
+        {
+          double x, y, z;
+          s_stream >> x >> y >> z;
+          points_enrichments[3 * i]     = x;
+          points_enrichments[3 * i + 1] = y;
+          points_enrichments[3 * i + 2] = z;
+        }
+      else
+        AssertThrow(false, ExcMessage("Dimension not implemented"));
+    }
+
+  // note vector of radii for predicates
+  for (unsigned int i = 0; i != n_enrichments; ++i)
+    {
+      read_next_proper_line(line);
+      s_stream.clear();
+      s_stream.str(line);
+
+      double r;
+      s_stream >> r;
+      radii_predicates.push_back(r);
+    }
+
+  // note vector of sigmas for rhs
+  for (unsigned int i = 0; i != n_enrichments; ++i)
+    {
+      read_next_proper_line(line);
+      s_stream.clear();
+      s_stream.str(line);
+
+      double r;
+      s_stream >> r;
+      sigmas.push_back(r);
+    }
+}
+
+
+
+ParameterCollection::ParameterCollection(
+  const int &                dim,
+  const double &             size,
+  const unsigned int &       shape,
+  const unsigned int &       global_refinement,
+  const unsigned int &       cycles,
+  const unsigned int &       fe_base_degree,
+  const unsigned int &       fe_enriched_degree,
+  const unsigned int &       max_iterations,
+  const double &             tolerance,
+  const std::string &        rhs_value_expr,
+  const std::string &        boundary_value_expr,
+  const std::string &        rhs_radial_problem,
+  const std::string &        boundary_radial_problem,
+  const std::string &        exact_soln_expr,
+  const unsigned int &       patches,
+  const unsigned int &       debug_level,
+  const unsigned int &       n_enrichments,
+  const std::vector<double> &points_enrichments,
+  const std::vector<double> &radii_predicates,
+  const std::vector<double> &sigmas) :
+  dim(dim),
+  size(size),
+  shape(shape),
+  global_refinement(global_refinement),
+  cycles(cycles),
+  fe_base_degree(fe_base_degree),
+  fe_enriched_degree(fe_enriched_degree),
+  max_iterations(max_iterations),
+  tolerance(tolerance),
+  rhs_value_expr(rhs_value_expr),
+  boundary_value_expr(boundary_value_expr),
+  rhs_radial_problem(rhs_radial_problem),
+  boundary_radial_problem(boundary_radial_problem),
+  exact_soln_expr(exact_soln_expr),
+  patches(patches),
+  debug_level(debug_level),
+  n_enrichments(n_enrichments),
+  points_enrichments(points_enrichments),
+  radii_predicates(radii_predicates),
+  sigmas(sigmas)
+{}
+
+
+
+void
+ParameterCollection::print()
+{
+  std::cout << "Dim : " << dim << std::endl
+            << "Size : " << size << std::endl
+            << "Shape : " << shape << std::endl
+            << "Global refinement : " << global_refinement << std::endl
+            << "Cycles : " << cycles << std::endl
+            << "FE base degree : " << fe_base_degree << std::endl
+            << "FE enriched degree : " << fe_enriched_degree << std::endl
+            << "Max Iterations : " << max_iterations << std::endl
+            << "Tolerance : " << tolerance << std::endl
+            << "rhs - main problem : " << rhs_value_expr << std::endl
+            << "boundary value - main problem : " << boundary_value_expr
+            << std::endl
+            << "rhs of radial problem : " << rhs_radial_problem << std::endl
+            << "boundary value of radial problem : " << boundary_radial_problem
+            << std::endl
+            << "exact solution expr : " << exact_soln_expr << std::endl
+            << "estimate exact solution using radial problem : "
+            << "Patches used for output: " << patches << std::endl
+            << "Debug level: " << debug_level << std::endl
+            << "Number of enrichments: " << n_enrichments << std::endl;
+
+  std::cout << "Enrichment points : " << std::endl;
+  for (unsigned int i = 0; i < points_enrichments.size(); i = i + dim)
+    {
+      for (int d = 0; d < dim; ++d)
+        std::cout << points_enrichments[i + d] << " ";
+
+      std::cout << std::endl;
+    }
+
+  std::cout << "Enrichment radii : " << std::endl;
+  for (auto r : radii_predicates)
+    std::cout << r << std::endl;
+
+  std::cout << "Sigma values of different sources : " << std::endl;
+  for (auto r : sigmas)
+    std::cout << r << std::endl;
+}
+
+
+
+/*
+ * EstimateEnrichmentFunction is used to estimate enrichment function by
+ * solveing a 1D poisson problem with right hand side and boundary
+ * expression provided as a string of single variable 'x', to be
+ * interpreted as distance from @par center.
+ *
+ * Eg: For a 2D poisson problem with right hand side expression R and boundary
+ * expression B given by functions R( x^2 + y^2) and B( x^2 + y^2 ), an
+ * equivalent radial problem can be solved using this class by setting
+ * @par rhs_expr = R( x^2)
+ * @par boundary_expr = B (x^2)
+ * Note that the original poisson problem is defined by right hand side
+ * and boundary expression dependent on square of distance (x^2 + y^2)
+ * from center.
+ */
+class EstimateEnrichmentFunction
+{
+public:
+  EstimateEnrichmentFunction(const Point<1> &   center,
+                             const double &     domain_size,
+                             const double &     sigma,
+                             const std::string &rhs_expr,
+                             const std::string &boundary_expr,
+                             const double &     refinement = 11);
+  EstimateEnrichmentFunction(const Point<1> &   center,
+                             const double &     left_bound,
+                             const double &     right_bound,
+                             const double &     sigma,
+                             const std::string &rhs_expr,
+                             const std::string &boundary_expr,
+                             const double &     refinement = 11);
+  ~EstimateEnrichmentFunction();
+  void
+  run();
+  void
+  evaluate_at_x_values(std::vector<double> &interpolation_points,
+                       std::vector<double> &interpolation_values);
+  double
+  value(const Point<1> &p, const unsigned int &component = 0);
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve();
+  void
+  refine_grid();
+  void
+                      output_results() const;
+  Point<1>            center;
+  double              domain_size;
+  double              left_bound, right_bound;
+  double              sigma;
+  std::string         rhs_expr;
+  std::string         boundary_expr;
+  std::vector<double> rhs_values;
+
+public:
+  unsigned int debug_level;
+
+private:
+  Triangulation<1>     triangulation;
+  unsigned int         refinement;
+  FE_Q<1>              fe;
+  DoFHandler<1>        dof_handler;
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+};
+
+EstimateEnrichmentFunction::EstimateEnrichmentFunction(
+  const Point<1> &   center,
+  const double &     domain_size,
+  const double &     sigma,
+  const std::string &rhs_expr,
+  const std::string &boundary_expr,
+  const double &     refinement) :
+  center(center),
+  domain_size(domain_size),
+  sigma(sigma),
+  rhs_expr(rhs_expr),
+  boundary_expr(boundary_expr),
+  debug_level(0),
+  refinement(refinement),
+  fe(1),
+  dof_handler(triangulation)
+{
+  left_bound  = center[0] - domain_size / 2;
+  right_bound = center[0] + domain_size / 2;
+}
+
+
+EstimateEnrichmentFunction::EstimateEnrichmentFunction(
+  const Point<1> &   center,
+  const double &     left_bound,
+  const double &     right_bound,
+  const double &     sigma,
+  const std::string &rhs_expr,
+  const std::string &boundary_expr,
+  const double &     refinement) :
+  center(center),
+  left_bound(left_bound),
+  right_bound(right_bound),
+  sigma(sigma),
+  rhs_expr(rhs_expr),
+  boundary_expr(boundary_expr),
+  debug_level(0),
+  refinement(refinement),
+  fe(1),
+  dof_handler(triangulation)
+{
+  domain_size = right_bound - left_bound;
+}
+
+
+void
+EstimateEnrichmentFunction::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, left_bound, right_bound);
+  triangulation.refine_global(refinement);
+}
+
+
+void
+EstimateEnrichmentFunction::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+  sparsity_pattern.copy_from(dsp);
+  system_matrix.reinit(sparsity_pattern);
+  solution.reinit(dof_handler.n_dofs());
+  system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+void
+EstimateEnrichmentFunction::assemble_system()
+{
+  QGauss<1>        quadrature_formula(2);
+  SigmaFunction<1> rhs;
+  rhs.initialize(center, sigma, rhs_expr);
+  FEValues<1>        fe_values(fe,
+                        quadrature_formula,
+                        update_values | update_gradients |
+                          update_quadrature_points | update_JxW_values);
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+      cell_matrix = 0;
+      cell_rhs    = 0;
+      rhs_values.resize(n_q_points);
+      rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
+
+      for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
+        {
+          double radius = center.distance(fe_values.quadrature_point(q_index));
+
+          //-1/r (r*u_r) = f form converts to
+          // r(u_r, v_r) = (r*f,v)
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                cell_matrix(i, j) +=
+                  (radius * fe_values.shape_grad(i, q_index) *
+                   fe_values.shape_grad(j, q_index)) *
+                  fe_values.JxW(q_index);
+              cell_rhs(i) +=
+                radius * (fe_values.shape_value(i, q_index) *
+                          rhs_values[q_index] * fe_values.JxW(q_index));
+            }
+        }
+      cell->get_dof_indices(local_dof_indices);
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        {
+          for (unsigned int j = 0; j < dofs_per_cell; ++j)
+            system_matrix.add(
+              local_dof_indices[i], local_dof_indices[j], cell_matrix(i, j));
+          system_rhs(local_dof_indices[i]) += cell_rhs(i);
+        }
+    }
+  std::map<types::global_dof_index, double> boundary_values;
+  SigmaFunction<1>                          boundary_func;
+  boundary_func.initialize(center, sigma, boundary_expr);
+
+  VectorTools::interpolate_boundary_values(
+    dof_handler, 0, boundary_func, boundary_values);
+  VectorTools::interpolate_boundary_values(
+    dof_handler, 1, boundary_func, boundary_values);
+
+  MatrixTools::apply_boundary_values(
+    boundary_values, system_matrix, solution, system_rhs);
+}
+
+
+void
+EstimateEnrichmentFunction::solve()
+{
+  SolverControl solver_control(50000, 1e-12, false, false);
+  SolverCG<>    solver(solver_control);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+}
+
+
+void
+EstimateEnrichmentFunction::refine_grid()
+{
+  Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+  KellyErrorEstimator<1>::estimate(dof_handler,
+                                   QGauss<1 - 1>(3),
+                                   typename FunctionMap<1>::type(),
+                                   solution,
+                                   estimated_error_per_cell);
+  GridRefinement::refine_and_coarsen_fixed_number(
+    triangulation, estimated_error_per_cell, 0.2, 0.01);
+  triangulation.execute_coarsening_and_refinement();
+}
+
+
+void
+EstimateEnrichmentFunction::output_results() const
+{
+  DataOut<1> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution, "solution");
+  data_out.build_patches();
+  std::ofstream output("radial_solution.vtk");
+  data_out.write_vtk(output);
+}
+
+
+void
+EstimateEnrichmentFunction::run()
+{
+  if (debug_level >= 1)
+    std::cout << "Solving problem in 1.: " << 1 << " with center: " << center
+              << ", size: " << domain_size << ", sigma: " << sigma << std::endl;
+
+  make_grid();
+
+  double old_value = 0, value = 1, relative_change = 1;
+  bool   start = true;
+  do
+    {
+      if (!start)
+        {
+          refine_grid();
+          ++refinement;
+        }
+
+      if (debug_level >= 1)
+        std::cout << "Refinement level: " << refinement << std::endl;
+
+      setup_system();
+      assemble_system();
+      solve();
+
+      value = VectorTools::point_value(dof_handler, solution, center);
+      if (!start)
+        {
+          relative_change = fabs((old_value - value) / old_value);
+        }
+      start     = false;
+      old_value = value;
+    }
+  while (relative_change > 0.005);
+
+  if (debug_level >= 1)
+    std::cout << "Radial solution at origin = " << value
+              << " after global refinement " << refinement << std::endl;
+
+  if (debug_level >= 1)
+    output_results();
+}
+
+
+void
+EstimateEnrichmentFunction::evaluate_at_x_values(
+  std::vector<double> &interpolation_points,
+  std::vector<double> &interpolation_values)
+{
+  if (interpolation_values.size() != interpolation_points.size())
+    interpolation_values.resize(interpolation_points.size());
+
+  // x varies from 0 to 2*sigma.
+  // factor 2 because once a cell is decided to be enriched based on its center,
+  // its quadrature points can cause x to be twice!
+  for (unsigned int i = 0; i != interpolation_values.size(); ++i)
+    {
+      double value = VectorTools::point_value(
+        dof_handler, solution, Point<1>(interpolation_points[i]));
+      interpolation_values[i] = value;
+    }
+}
+
+
+double
+EstimateEnrichmentFunction::value(const Point<1> &    p,
+                                  const unsigned int &component)
+{
+  return VectorTools::point_value(dof_handler, solution, p);
+}
+
+
+EstimateEnrichmentFunction::~EstimateEnrichmentFunction()
+{
+  triangulation.clear();
+}
+
+
+
+template <int dim>
+void
+plot_shape_function(hp::DoFHandler<dim> &dof_handler, unsigned int patches = 5)
+{
+  std::cout << "...start plotting shape function" << std::endl;
+  std::cout << "Patches for output: " << patches << std::endl;
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  constraints.close();
+
+  // find set of dofs which belong to enriched cells
+  std::set<unsigned int> enriched_cell_dofs;
+  for (auto cell : dof_handler.active_cell_iterators())
+    if (cell->active_fe_index() != 0)
+      {
+        unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+        std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+        cell->get_dof_indices(local_dof_indices);
+        enriched_cell_dofs.insert(local_dof_indices.begin(),
+                                  local_dof_indices.end());
+      }
+
+  // output to check if all is good:
+  std::vector<Vector<double>> shape_functions;
+  std::vector<std::string>    names;
+  for (auto dof : enriched_cell_dofs)
+    {
+      Vector<double> shape_function;
+      shape_function.reinit(dof_handler.n_dofs());
+      shape_function[dof] = 1.0;
+
+      // if the dof is constrained, first output unconstrained vector
+      names.push_back(std::string("C_") +
+                      dealii::Utilities::int_to_string(dof, 2));
+      shape_functions.push_back(shape_function);
+
+      //      names.push_back(std::string("UC_") +
+      //                      dealii::Utilities::int_to_string(s,2));
+
+      //      // make continuous/constraint:
+      //      constraints.distribute(shape_function);
+      //      shape_functions.push_back(shape_function);
+    }
+
+  if (dof_handler.n_dofs() < 100)
+    {
+      std::cout << "...start printing support points" << std::endl;
+
+      std::map<types::global_dof_index, Point<dim>> support_points;
+      MappingQ1<dim>                                mapping;
+      hp::MappingCollection<dim>                    hp_mapping;
+      for (unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i)
+        hp_mapping.push_back(mapping);
+      DoFTools::map_dofs_to_support_points(
+        hp_mapping, dof_handler, support_points);
+
+      const std::string base_filename =
+        "DOFs" + dealii::Utilities::int_to_string(dim) + "_p" +
+        dealii::Utilities::int_to_string(0);
+
+      const std::string filename = base_filename + ".gp";
+      std::ofstream     f(filename.c_str());
+
+      f << "set terminal png size 400,410 enhanced font \"Helvetica,8\""
+        << std::endl
+        << "set output \"" << base_filename << ".png\"" << std::endl
+        << "set size square" << std::endl
+        << "set view equal xy" << std::endl
+        << "unset xtics                                                                                   "
+        << std::endl
+        << "unset ytics" << std::endl
+        << "unset grid" << std::endl
+        << "unset border" << std::endl
+        << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle"
+        << std::endl;
+      GridOut grid_out;
+      grid_out.write_gnuplot(dof_handler.get_triangulation(), f);
+      f << "e" << std::endl;
+
+      DoFTools::write_gnuplot_dof_support_point_info(f, support_points);
+
+      f << "e" << std::endl;
+
+      std::cout << "...finished printing support points" << std::endl;
+    }
+
+  DataOut<dim, hp::DoFHandler<dim>> data_out;
+  data_out.attach_dof_handler(dof_handler);
+
+  // get material ids:
+  Vector<float> fe_index(dof_handler.get_triangulation().n_active_cells());
+  for (auto cell : dof_handler.active_cell_iterators())
+    {
+      fe_index[cell->active_cell_index()] = cell->active_fe_index();
+    }
+  data_out.add_data_vector(fe_index, "fe_index");
+
+  for (unsigned int i = 0; i < shape_functions.size(); i++)
+    data_out.add_data_vector(shape_functions[i], names[i]);
+
+  data_out.build_patches(patches);
+
+  std::string   filename = "shape_functions.vtu";
+  std::ofstream output(filename.c_str());
+  data_out.write_vtu(output);
+
+  std::cout << "...finished plotting shape functions" << std::endl;
+}
+
+
+
+template <int dim>
+using predicate_function =
+  std::function<bool(const typename Triangulation<dim>::cell_iterator &)>;
+
+
+
+/**
+ * Main class
+ */
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem();
+  LaplaceProblem(const ParameterCollection &prm);
+  virtual ~LaplaceProblem();
+  void
+  run();
+
+protected:
+  void
+  initialize();
+  void
+  build_fe_space();
+  virtual void
+  make_enrichment_functions();
+  void
+  setup_system();
+
+private:
+  void
+  build_tables();
+  void
+  assemble_system();
+  unsigned int
+  solve();
+  void
+  refine_grid();
+  void
+  output_results(const unsigned int cycle);
+  void
+  process_solution();
+
+protected:
+  ParameterCollection prm;
+  unsigned int        n_enriched_cells;
+
+  Triangulation<dim>  triangulation;
+  hp::DoFHandler<dim> dof_handler;
+
+  std::shared_ptr<const hp::FECollection<dim>> fe_collection;
+  hp::QCollection<dim>                         q_collection;
+
+  FE_Q<dim>       fe_base;
+  FE_Q<dim>       fe_enriched;
+  FE_Nothing<dim> fe_nothing;
+
+  IndexSet locally_owned_dofs;
+  IndexSet locally_relevant_dofs;
+
+  PETScWrappers::MPI::SparseMatrix system_matrix;
+  PETScWrappers::MPI::Vector       solution;
+  Vector<double>                   localized_solution;
+  PETScWrappers::MPI::Vector       system_rhs;
+
+  ConstraintMatrix constraints;
+
+  MPI_Comm           mpi_communicator;
+  const unsigned int n_mpi_processes;
+  const unsigned int this_mpi_process;
+
+  ConditionalOStream pcout;
+
+  std::vector<SigmaFunction<dim>> vec_rhs;
+
+  using cell_iterator_function = std::function<Function<dim> *(
+    const typename hp::DoFHandler<dim>::active_cell_iterator &)>;
+
+  std::vector<std::shared_ptr<Function<dim>>> vec_enrichments;
+  std::vector<predicate_function<dim>>        vec_predicates;
+  std::vector<cell_iterator_function>         color_enrichments;
+
+  // output vectors. size triangulation.n_active_cells()
+  // change to Vector
+  std::vector<Vector<float>> predicate_output;
+  Vector<float>              color_output;
+  Vector<float>              vec_fe_index;
+  Vector<float>              mat_id;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem() :
+  prm(),
+  n_enriched_cells(0),
+  dof_handler(triangulation),
+  fe_base(prm.fe_base_degree),
+  fe_enriched(prm.fe_enriched_degree),
+  fe_nothing(1, true),
+  mpi_communicator(MPI_COMM_WORLD),
+  n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator)),
+  this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator)),
+  pcout(std::cout, (this_mpi_process == 0))
+{
+  prm.print();
+
+  pcout << "...default parameters set" << std::endl;
+}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem(const ParameterCollection &_par) :
+  prm(_par),
+  n_enriched_cells(0),
+  dof_handler(triangulation),
+  fe_base(prm.fe_base_degree),
+  fe_enriched(prm.fe_enriched_degree),
+  fe_nothing(1, true),
+  mpi_communicator(MPI_COMM_WORLD),
+  n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator)),
+  this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator)),
+  pcout(std::cout, (this_mpi_process == 0) && (prm.debug_level >= 1))
+{
+  AssertThrow(prm.dim == dim, ExcMessage("parameter file dim != problem dim"));
+  prm.print();
+  pcout << "...parameters set" << std::endl;
+}
+
+
+/*
+ * Construct basic grid, vector of predicate functions and
+ * right hand side function of the problem.
+ */
+template <int dim>
+void
+LaplaceProblem<dim>::initialize()
+{
+  pcout << "...Start initializing" << std::endl;
+
+  /*
+   * set up basic grid which is a hyper cube or hyper ball based on
+   * parameter file. Refine as per the global refinement value in the
+   * parameter file.
+   *
+   */
+  if (prm.shape == 1)
+    GridGenerator::hyper_cube(triangulation, -prm.size / 2.0, prm.size / 2.0);
+  else if (prm.shape == 0)
+    {
+      Point<dim> center = Point<dim>();
+      GridGenerator::hyper_ball(triangulation, center, prm.size / 2.0);
+      triangulation.set_all_manifold_ids_on_boundary(0);
+      static SphericalManifold<dim> spherical_manifold(center);
+      triangulation.set_manifold(0, spherical_manifold);
+    }
+  else
+    AssertThrow(false, ExcMessage("Shape not implemented."));
+  triangulation.refine_global(prm.global_refinement);
+
+  /*
+   * Ensure that num of radii, sigma and focal points of enrichment domains
+   * provided matches the number of predicates. Since focal points of
+   * enrichment domains are stored as vector of numbers, the dimension of
+   * points_enrichments is dim times the n_enrichments in prm file.
+   */
+  Assert(
+    prm.points_enrichments.size() / dim == prm.n_enrichments &&
+      prm.radii_predicates.size() == prm.n_enrichments &&
+      prm.sigmas.size() == prm.n_enrichments,
+    ExcMessage(
+      "Number of enrichment points, predicate radii and sigmas should be equal"));
+
+  /*
+   * Construct vector of predicate functions, where a function at index i
+   * returns true for a given cell if it belongs to enrichment domain i.
+   * The decision is based on cell center being at a distance within radius
+   * of the domain from focal point of the domain.
+   */
+  for (unsigned int i = 0; i != prm.n_enrichments; ++i)
+    {
+      Point<dim> p;
+      prm.set_enrichment_point(p, i);
+      vec_predicates.push_back(
+        EnrichmentPredicate<dim>(p, prm.radii_predicates[i]));
+    }
+
+  /*
+   * Construct a vector of right hand side functions where the actual right hand
+   * side of problem is evaluated through summation of contributions from each
+   * of the functions in the vector. Each function i at a given point is
+   * evaluated with respect to point's position {x_i, y_i, z_i} relative to
+   * focal point of enrichment domain i. The value is then a function of x_i,
+   * y_i, z_i and sigmas[i] given by the rhs_value_expr[i] in parameter file.
+   */
+  vec_rhs.resize(prm.n_enrichments);
+  for (unsigned int i = 0; i != prm.n_enrichments; ++i)
+    {
+      Point<dim> p;
+      prm.set_enrichment_point(p, i);
+      vec_rhs[i].initialize(p, prm.sigmas[i], prm.rhs_value_expr);
+    }
+
+  pcout << "...finish initializing" << std::endl;
+}
+
+
+
+/*
+ * Create a enrichment function associated with each enrichment domain.
+ * Here the approximate solution is assumed to be only dependent on
+ * distance (r) from the focal point of enrichment domain. This is
+ * true for poisson problem (a linear PDE) provided the right hand side is a
+ * radial function and drops exponentially for points farther from focal
+ * point of enrichment domain.
+ */
+template <int dim>
+void
+LaplaceProblem<dim>::make_enrichment_functions()
+{
+  pcout << "!!! Make enrichment function called" << std::endl;
+
+  for (unsigned int i = 0; i < vec_predicates.size(); ++i)
+    {
+      /*
+       * Formulate a 1d/radial problem with center and size appropriate
+       * to cover the enrichment domain. The function determining
+       * right hand side and boundary value of the problem is provided in
+       * parameter file. The sigma governing this function is the same as
+       * sigma provided for the corresponding enrichment domain and predicate
+       * function for which the enrichment function is to be estimated.
+       *
+       * The center is 0 since the function is evaluated with respect to
+       * relative position from focal point anyway.
+       *
+       * For hexahedral cells, dimension can extend upto sqrt(3) < 2 times!
+       * So take a factor of 4 as size of the problem. This ensures that
+       * enrichment function can be evaluated at all points in the enrichment
+       * domain
+       */
+      double center = 0;
+      double sigma  = prm.sigmas[i];
+      double size   = prm.radii_predicates[i] * 4;
+
+      /*
+       * Radial problem is solved only when enrichment domain has a positive
+       * radius and hence non-empty domain.
+       */
+      if (prm.radii_predicates[i] != 0)
+        {
+          EstimateEnrichmentFunction radial_problem(
+            Point<1>(center),
+            size,
+            sigma,
+            prm.rhs_radial_problem,
+            prm.boundary_radial_problem);
+          radial_problem.debug_level = prm.debug_level; // print output
+          radial_problem.run();
+          pcout << "solved problem with "
+                << "x and sigma : " << center << ", " << sigma << std::endl;
+
+          // make points at which solution needs to interpolated
+          std::vector<double> interpolation_points, interpolation_values;
+          double              cut_point = 3 * sigma;
+          unsigned int        n1 = 15, n2 = 15;
+          double              radius      = size / 2;
+          double              right_bound = center + radius;
+          double h1 = cut_point / n1, h2 = (radius - cut_point) / n2;
+          for (double p = center; p < center + cut_point; p += h1)
+            interpolation_points.push_back(p);
+          for (double p = center + cut_point; p < right_bound; p += h2)
+            interpolation_points.push_back(p);
+          interpolation_points.push_back(right_bound);
+
+          // add enrichment function only when predicate radius is non-zero
+          radial_problem.evaluate_at_x_values(interpolation_points,
+                                              interpolation_values);
+
+
+          // construct enrichment function and push
+          Point<dim> p;
+          prm.set_enrichment_point(p, i);
+          SplineEnrichmentFunction<dim> func(
+            p, interpolation_points, interpolation_values);
+          vec_enrichments.push_back(
+            std::make_shared<SplineEnrichmentFunction<dim>>(func));
+        }
+      else
+        {
+          pcout << "Dummy function added at " << i << std::endl;
+          ConstantFunction<dim> func(0);
+          vec_enrichments.push_back(
+            std::make_shared<ConstantFunction<dim>>(func));
+        }
+    }
+}
+
+
+
+/*
+ * Since each enrichment domain has different enrichment function
+ * associated with it and the cells common to different enrichment
+ * domains need to treated differently, we use helper function in
+ * ColorEnriched namespace to construct finite element space and necessary
+ * data structures required by it. The helper function is also used to
+ * set FE indices of dof handler cells. We also set the cell with a unique
+ * material id for now, which is used to map cells with a pairs of
+ * color and corresponding enrichment function. All this is internally
+ * used to figure out the correct set of enrichment functions to be used
+ * for a cell.
+ *
+ * The quadrature point collection, with size equal to FE collection
+ * is also constructed here.
+ */
+template <int dim>
+void
+LaplaceProblem<dim>::build_fe_space()
+{
+  pcout << "...building fe space" << std::endl;
+
+  make_enrichment_functions();
+  static ColorEnriched::Helper<dim> fe_space(
+    fe_base, fe_enriched, vec_predicates, vec_enrichments);
+  fe_collection = std::make_shared<const hp::FECollection<dim>>(
+    fe_space.build_fe_collection(dof_handler));
+  pcout << "size of fe collection: " << fe_collection->size() << std::endl;
+
+  if (prm.debug_level == 9)
+    {
+      if (triangulation.n_active_cells() < 100)
+        {
+          pcout << "...start print fe indices" << std::endl;
+
+          // print fe index
+          const std::string base_filename =
+            "fe_indices" + dealii::Utilities::int_to_string(dim) + "_p" +
+            dealii::Utilities::int_to_string(0);
+          const std::string filename = base_filename + ".gp";
+          std::ofstream     f(filename.c_str());
+
+          f << "set terminal png size 400,410 enhanced font \"Helvetica,8\""
+            << std::endl
+            << "set output \"" << base_filename << ".png\"" << std::endl
+            << "set size square" << std::endl
+            << "set view equal xy" << std::endl
+            << "unset xtics" << std::endl
+            << "unset ytics" << std::endl
+            << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle"
+            << std::endl;
+          GridOut().write_gnuplot(triangulation, f);
+          f << "e" << std::endl;
+
+          for (auto it : dof_handler.active_cell_iterators())
+            f << it->center() << " \"" << it->active_fe_index() << "\"\n";
+
+          f << std::flush << "e" << std::endl;
+          pcout << "...finished print fe indices" << std::endl;
+        }
+
+      if (triangulation.n_active_cells() < 100)
+        {
+          pcout << "...start print cell indices" << std::endl;
+
+          // print cell ids
+          const std::string base_filename =
+            "cell_id" + dealii::Utilities::int_to_string(dim) + "_p" +
+            dealii::Utilities::int_to_string(0);
+          const std::string filename = base_filename + ".gp";
+          std::ofstream     f(filename.c_str());
+
+          f << "set terminal png size 400,410 enhanced font \"Helvetica,8\""
+            << std::endl
+            << "set output \"" << base_filename << ".png\"" << std::endl
+            << "set size square" << std::endl
+            << "set view equal xy" << std::endl
+            << "unset xtics" << std::endl
+            << "unset ytics" << std::endl
+            << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle"
+            << std::endl;
+          GridOut().write_gnuplot(triangulation, f);
+          f << "e" << std::endl;
+
+          for (auto it : dof_handler.active_cell_iterators())
+            f << it->center() << " \"" << it->index() << "\"\n";
+
+          f << std::flush << "e" << std::endl;
+
+          pcout << "...end print cell indices" << std::endl;
+        }
+    }
+
+  // q collections the same size as different material identities
+  q_collection.push_back(QGauss<dim>(4));
+  for (unsigned int i = 1; i < fe_collection->size(); ++i)
+    q_collection.push_back(QGauss<dim>(10));
+
+  pcout << "...building fe space" << std::endl;
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::setup_system()
+{
+  pcout << "...start setup system" << std::endl;
+
+  GridTools::partition_triangulation(n_mpi_processes, triangulation);
+
+  dof_handler.distribute_dofs(*fe_collection);
+
+  DoFRenumbering::subdomain_wise(dof_handler);
+  std::vector<IndexSet> locally_owned_dofs_per_proc =
+    DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
+  locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
+  locally_relevant_dofs.clear();
+  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+  constraints.clear();
+  constraints.reinit(locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+  SigmaFunction<dim> boundary_value_func;
+  Point<dim>         p;
+  prm.set_enrichment_point(p, 0);
+  boundary_value_func.initialize(p, prm.sigmas[0], prm.boundary_value_expr);
+
+
+  VectorTools::interpolate_boundary_values(
+    dof_handler, 0, boundary_value_func, constraints);
+  constraints.close();
+
+  // Initialise the stiffness and mass matrices
+  DynamicSparsityPattern dsp(locally_relevant_dofs);
+  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+  std::vector<types::global_dof_index> n_locally_owned_dofs(n_mpi_processes);
+  for (unsigned int i = 0; i < n_mpi_processes; ++i)
+    n_locally_owned_dofs[i] = locally_owned_dofs_per_proc[i].n_elements();
+
+  SparsityTools::distribute_sparsity_pattern(
+    dsp, n_locally_owned_dofs, mpi_communicator, locally_relevant_dofs);
+
+  system_matrix.reinit(
+    locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
+
+  solution.reinit(locally_owned_dofs, mpi_communicator);
+  system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+  pcout << "...finished setup system" << std::endl;
+}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::assemble_system()
+{
+  pcout << "...assemble system" << std::endl;
+
+  system_matrix = 0;
+  system_rhs    = 0;
+
+  FullMatrix<double> cell_system_matrix;
+  Vector<double>     cell_rhs;
+
+  std::vector<types::global_dof_index> local_dof_indices;
+
+  std::vector<double> rhs_value, tmp_rhs_value;
+
+  hp::FEValues<dim> fe_values_hp(*fe_collection,
+                                 q_collection,
+                                 update_values | update_gradients |
+                                   update_quadrature_points |
+                                   update_JxW_values);
+
+
+  for (auto cell : dof_handler.active_cell_iterators())
+    if (cell->subdomain_id() == this_mpi_process)
+      {
+        fe_values_hp.reinit(cell);
+        const FEValues<dim> &fe_values = fe_values_hp.get_present_fe_values();
+
+        const unsigned int &dofs_per_cell = cell->get_fe().dofs_per_cell;
+        const unsigned int &n_q_points    = fe_values.n_quadrature_points;
+
+        /*
+         * Initialize rhs values vector to zero. Add values calculated
+         * from each of different rhs functions (vec_rhs).
+         */
+        rhs_value.assign(n_q_points, 0);
+        tmp_rhs_value.assign(n_q_points, 0);
+        for (unsigned int i = 0; i < vec_rhs.size(); ++i)
+          {
+            vec_rhs[i].value_list(fe_values.get_quadrature_points(),
+                                  tmp_rhs_value);
+
+            // add tmp to the total one at quadrature points
+            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+              {
+                rhs_value[q_point] += tmp_rhs_value[q_point];
+              }
+          }
+
+        local_dof_indices.resize(dofs_per_cell);
+        cell_system_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_rhs.reinit(dofs_per_cell);
+
+        cell_system_matrix = 0;
+        cell_rhs           = 0;
+
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              for (unsigned int j = i; j < dofs_per_cell; ++j)
+                cell_system_matrix(i, j) +=
+                  (fe_values.shape_grad(i, q_point) *
+                   fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
+
+              cell_rhs(i) +=
+                (rhs_value[q_point] * fe_values.shape_value(i, q_point) *
+                 fe_values.JxW(q_point));
+            }
+
+        // exploit symmetry
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          for (unsigned int j = i; j < dofs_per_cell; ++j)
+            cell_system_matrix(j, i) = cell_system_matrix(i, j);
+
+        cell->get_dof_indices(local_dof_indices);
+
+        constraints.distribute_local_to_global(cell_system_matrix,
+                                               cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix,
+                                               system_rhs);
+      }
+
+  system_matrix.compress(VectorOperation::add);
+  system_rhs.compress(VectorOperation::add);
+
+  pcout << "...finished assemble system" << std::endl;
+}
+
+template <int dim>
+unsigned int
+LaplaceProblem<dim>::solve()
+{
+  pcout << "...solving" << std::endl;
+  SolverControl solver_control(prm.max_iterations, prm.tolerance, false, false);
+  PETScWrappers::SolverCG cg(solver_control, mpi_communicator);
+
+  // choose preconditioner
+#define amg
+#ifdef amg
+  PETScWrappers::PreconditionBoomerAMG::AdditionalData additional_data;
+  additional_data.symmetric_operator = true;
+  PETScWrappers::PreconditionBoomerAMG preconditioner(system_matrix,
+                                                      additional_data);
+#else
+  PETScWrappers::PreconditionJacobi preconditioner(system_matrix);
+#endif
+
+  cg.solve(system_matrix, solution, system_rhs, preconditioner);
+
+  Vector<double> local_soln(solution);
+
+  constraints.distribute(local_soln);
+  solution = local_soln;
+  pcout << "...finished solving" << std::endl;
+  return solver_control.last_step();
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::refine_grid()
+{
+  const Vector<double> localized_solution(solution);
+  Vector<float>        local_error_per_cell(triangulation.n_active_cells());
+
+  hp::QCollection<dim - 1> q_collection_face;
+  for (unsigned int i = 0; i < q_collection.size(); ++i)
+    q_collection_face.push_back(QGauss<dim - 1>(1));
+
+  KellyErrorEstimator<dim>::estimate(dof_handler,
+                                     q_collection_face,
+                                     typename FunctionMap<dim>::type(),
+                                     localized_solution,
+                                     local_error_per_cell,
+                                     ComponentMask(),
+                                     nullptr,
+                                     n_mpi_processes,
+                                     this_mpi_process);
+  const unsigned int n_local_cells =
+    GridTools::count_cells_with_subdomain_association(triangulation,
+                                                      this_mpi_process);
+  PETScWrappers::MPI::Vector distributed_all_errors(
+    mpi_communicator, triangulation.n_active_cells(), n_local_cells);
+  for (unsigned int i = 0; i < local_error_per_cell.size(); ++i)
+    if (local_error_per_cell(i) != 0)
+      distributed_all_errors(i) = local_error_per_cell(i);
+  distributed_all_errors.compress(VectorOperation::insert);
+  const Vector<float> localized_all_errors(distributed_all_errors);
+  GridRefinement::refine_and_coarsen_fixed_fraction(
+    triangulation, localized_all_errors, 0.85, 0);
+  triangulation.execute_coarsening_and_refinement();
+  ++prm.global_refinement;
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::output_results(const unsigned int cycle)
+{
+  pcout << "...output results" << std::endl;
+  pcout << "Patches used: " << prm.patches << std::endl;
+
+  Vector<double>     exact_soln_vector, error_vector;
+  SigmaFunction<dim> exact_solution;
+
+  if (prm.exact_soln_expr != "")
+    {
+      // create exact solution vector
+      exact_soln_vector.reinit(dof_handler.n_dofs());
+      exact_solution.initialize(
+        Point<dim>(), prm.sigmas[0], prm.exact_soln_expr);
+      VectorTools::project(dof_handler,
+                           constraints,
+                           q_collection,
+                           exact_solution,
+                           exact_soln_vector);
+
+      // create error vector
+      error_vector.reinit(dof_handler.n_dofs());
+      Vector<double> full_solution(localized_solution);
+      error_vector += full_solution;
+      error_vector -= exact_soln_vector;
+    }
+
+  Assert(cycle < 10, ExcNotImplemented());
+  if (this_mpi_process == 0)
+    {
+      std::string filename = "solution-";
+      filename += Utilities::to_string(cycle);
+      filename += ".vtk";
+      std::ofstream output(filename.c_str());
+
+      DataOut<dim, hp::DoFHandler<dim>> data_out;
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(localized_solution, "solution");
+      if (prm.exact_soln_expr != "")
+        {
+          data_out.add_data_vector(exact_soln_vector, "exact_solution");
+          data_out.add_data_vector(error_vector, "error_vector");
+        }
+      data_out.build_patches(prm.patches);
+      data_out.write_vtk(output);
+      output.close();
+    }
+  pcout << "...finished output results" << std::endl;
+}
+
+
+
+// use this only when exact solution is known
+template <int dim>
+void
+LaplaceProblem<dim>::process_solution()
+{
+  Vector<float> difference_per_cell(triangulation.n_active_cells());
+  double        L2_error, H1_error;
+
+  if (!prm.exact_soln_expr.empty())
+    {
+      pcout << "...using exact solution for error calculation" << std::endl;
+
+      SigmaFunction<dim> exact_solution;
+      exact_solution.initialize(
+        Point<dim>(), prm.sigmas[0], prm.exact_soln_expr);
+
+      VectorTools::integrate_difference(dof_handler,
+                                        localized_solution,
+                                        exact_solution,
+                                        difference_per_cell,
+                                        q_collection,
+                                        VectorTools::L2_norm);
+      L2_error = VectorTools::compute_global_error(
+        triangulation, difference_per_cell, VectorTools::L2_norm);
+
+      VectorTools::integrate_difference(dof_handler,
+                                        localized_solution,
+                                        exact_solution,
+                                        difference_per_cell,
+                                        q_collection,
+                                        VectorTools::H1_norm);
+      H1_error = VectorTools::compute_global_error(
+        triangulation, difference_per_cell, VectorTools::H1_norm);
+    }
+
+  pcout << "refinement h_smallest Dofs L2_norm H1_norm" << std::endl;
+  pcout << prm.global_refinement << " "
+        << prm.size / std::pow(2.0, prm.global_refinement) << " "
+        << dof_handler.n_dofs() << " " << L2_error << " " << H1_error
+        << std::endl;
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::run()
+{
+  pcout << "...run problem" << std::endl;
+  double norm_soln_old(0), norm_rel_change_old(1);
+
+  // Run making grids and building fe space only once.
+  initialize();
+  build_fe_space();
+
+
+  if (this_mpi_process == 0)
+    deallog << "Solving problem with number of sources: " << prm.n_enrichments
+            << std::endl;
+
+  for (unsigned int cycle = 0; cycle <= prm.cycles; ++cycle)
+    {
+      pcout << "Cycle " << cycle << std::endl;
+
+      setup_system();
+
+      pcout << "Number of active cells:       "
+            << triangulation.n_active_cells() << std::endl
+            << "Number of degrees of freedom: " << dof_handler.n_dofs()
+            << std::endl;
+
+      if (prm.debug_level == 9 && this_mpi_process == 0)
+        plot_shape_function<dim>(dof_handler);
+
+      assemble_system();
+      auto n_iterations = solve();
+      pcout << "Number of iterations: " << n_iterations << std::endl;
+      localized_solution.reinit(dof_handler.n_dofs());
+      localized_solution = solution;
+      double value =
+        VectorTools::point_value(dof_handler, localized_solution, Point<dim>());
+      pcout << "Solution at origin:   " << value << std::endl;
+
+
+      // calculate L2 norm of solution
+      if (this_mpi_process == 0)
+        {
+          pcout << "calculating L2 norm of soln" << std::endl;
+          double        norm_soln_new(0), norm_rel_change_new(0);
+          Vector<float> difference_per_cell(triangulation.n_active_cells());
+          VectorTools::integrate_difference(dof_handler,
+                                            localized_solution,
+                                            ZeroFunction<dim>(),
+                                            difference_per_cell,
+                                            q_collection,
+                                            VectorTools::H1_norm);
+          norm_soln_new = VectorTools::compute_global_error(
+            triangulation, difference_per_cell, VectorTools::H1_norm);
+          // relative change can only be calculated for cycle > 0
+          if (cycle > 0)
+            {
+              norm_rel_change_new =
+                std::abs((norm_soln_new - norm_soln_old) / norm_soln_old);
+              pcout << "relative change of solution norm "
+                    << norm_rel_change_new << std::endl;
+            }
+
+          // moniter relative change of norm in later stages
+          if (cycle > 1)
+            {
+              deallog << (norm_rel_change_new < norm_rel_change_old)
+                      << std::endl;
+            }
+
+          norm_soln_old = norm_soln_new;
+
+          // first sample of relative change of norm comes only cycle = 1
+          if (cycle > 0)
+            norm_rel_change_old = norm_rel_change_new;
+
+          pcout << "End of L2 calculation" << std::endl;
+        }
+
+      if (prm.debug_level >= 2 && this_mpi_process == 0)
+        output_results(cycle);
+
+      // Donot refine if loop is at the end
+      if (cycle != prm.cycles)
+        refine_grid();
+
+      pcout << "...step run complete" << std::endl;
+    }
+  pcout << "...finished run problem" << std::endl;
+}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem()
+{
+  dof_handler.clear();
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+  {
+    /*
+     * Case 1: single source with known solution
+     */
+    {
+      ParameterCollection prm(
+        2,     // dimension
+        2,     // domain size
+        1,     // cube shape
+        3,     // global refinement
+        5,     // num of cycles grid is refined and solved again
+        1,     // fe base degree
+        1,     // fe enriched degree
+        50000, // max iterations
+        1e-9,  // tolerance
+        // rhs value
+        "-(exp(-(x*x + y*y)/(2*sigma*sigma))*(- 2*sigma*sigma + x*x + y*y))/(2*sigma*sigma*sigma*sigma*sigma*sigma*pi)",
+        // boundary value
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))",
+        // rhs value for radial problem solved to find enrichment function
+        "-(exp(-(x*x)/(2*sigma*sigma))*(- 2*sigma*sigma + x*x))/(2*sigma*sigma*sigma*sigma*sigma*sigma*pi)",
+        // boundary value for radial problem
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))",
+        // exact solution expression. If null nothing is done
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))",
+        1, // patches
+        1, // debug level
+        1, // num enrichments
+        // enrichment points interpreted 2 at a time if dimension is 2
+        {0, 0},
+        // radii defining different predicates
+        {0.4},
+        // sigmas defining different predicates
+        {0.1});
+
+
+      if (prm.dim == 2)
+        {
+          LaplaceProblem<2> problem(prm);
+          problem.run();
+        }
+      else if (prm.dim == 3)
+        {
+          LaplaceProblem<3> problem(prm);
+          problem.run();
+        }
+      else
+        AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3"));
+    }
+
+
+
+    /*
+     * Case 2: 3 sources
+     */
+    {
+      ParameterCollection prm(
+        2,     // dimension
+        4,     // domain size
+        1,     // cube shape
+        3,     // global refinement
+        4,     // num of cycles grid is refined and solved again
+        1,     // fe base degree
+        1,     // fe enriched degree
+        50000, // max iterations
+        1e-9,  // tolerance
+        // rhs value
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))",
+        // boundary value
+        "0",
+        // rhs value for radial problem solved to find enrichment function
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))",
+        // boundary value for radial problem
+        "0",
+        // exact solution expression. If null nothing is done
+        "",
+        1, // patches
+        1, // debug level
+        3, // num enrichments
+        // enrichment points interpreted 2 at a time if dimension is 2
+        {0.5, 0.5, 0, 0, -1, -1},
+        // radii defining different predicates
+        {0.4, 0.4, 0.4},
+        // sigmas defining different predicates
+        {0.1, 0.1, 0.1});
+
+
+      if (prm.dim == 2)
+        {
+          LaplaceProblem<2> problem(prm);
+          problem.run();
+        }
+      else if (prm.dim == 3)
+        {
+          LaplaceProblem<3> problem(prm);
+          problem.run();
+        }
+      else
+        AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3"));
+    }
+
+
+    /*
+     * Case 3: five sources 3d
+     */
+    {
+      ParameterCollection prm(
+        3,     // dimension
+        8,     // domain size
+        1,     // cube shape
+        4,     // global refinement
+        2,     // num of cycles grid is refined and solved again
+        1,     // fe base degree
+        1,     // fe enriched degree
+        50000, // max iterations
+        1e-9,  // tolerance
+        // rhs value
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y + z*z)/(2*sigma*sigma))",
+        // boundary value
+        "0",
+        // rhs value for radial problem solved to find enrichment function
+        "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))",
+        // boundary value for radial problem
+        "0",
+        // exact solution expression. If null nothing is done
+        "",
+        1, // patches
+        1, // debug level
+        5, // num enrichments
+        // enrichment points interpreted 3 at a time if dimension is 3
+        {1.5, 1.5, 1.5, 1, 1, 1, 0, 0, 0, -1, -1, -1, 1, -1, -1},
+        // radii defining different predicates
+        {0.45, 0.45, 0.45, 0.45, 0.45},
+        // sigmas defining different predicates
+        {0.1, 0.1, 0.1, 0.1, 0.1});
+
+
+      if (prm.dim == 2)
+        {
+          LaplaceProblem<2> problem(prm);
+          problem.run();
+        }
+      else if (prm.dim == 3)
+        {
+          LaplaceProblem<3> problem(prm);
+          problem.run();
+        }
+      else
+        AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3"));
+    }
+  }
+}
diff --git a/tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output b/tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..c5bd2f2
--- /dev/null
@@ -0,0 +1,14 @@
+
+DEAL:0::Solving problem with number of sources: 1
+DEAL:0::1
+DEAL:0::1
+DEAL:0::1
+DEAL:0::1
+DEAL:0::Solving problem with number of sources: 3
+DEAL:0::1
+DEAL:0::1
+DEAL:0::1
+DEAL:0::Solving problem with number of sources: 5
+DEAL:0::1
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.