--- /dev/null
+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)
#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>
&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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
#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
}
+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"
{
#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
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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:
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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:
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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"));
+ }
+ }
+}
--- /dev/null
+
+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
+
+
projects / dealii.git / commitdiff
