]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add FE_Enriched finite element 3115/head
authorDenis Davydov <davydden@gmail.com>
Wed, 28 Sep 2016 09:52:07 +0000 (11:52 +0200)
committerDenis Davydov <davydden@gmail.com>
Wed, 28 Sep 2016 13:35:19 +0000 (15:35 +0200)
35 files changed:
doc/doxygen/images/fe_enriched_1d.png [new file with mode: 0644]
doc/doxygen/images/fe_enriched_h-refinement.png [new file with mode: 0644]
doc/doxygen/stylesheet.css
doc/news/changes.h
include/deal.II/fe/fe_enriched.h [new file with mode: 0644]
include/deal.II/fe/fe_system.h
source/fe/CMakeLists.txt
source/fe/fe_enriched.cc [new file with mode: 0644]
source/fe/fe_enriched.inst.in [new file with mode: 0644]
tests/fe/fe_enriched_01.cc [new file with mode: 0644]
tests/fe/fe_enriched_01.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_02.cc [new file with mode: 0644]
tests/fe/fe_enriched_02.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_03.cc [new file with mode: 0644]
tests/fe/fe_enriched_03.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_04.cc [new file with mode: 0644]
tests/fe/fe_enriched_04.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_05.cc [new file with mode: 0644]
tests/fe/fe_enriched_05.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_06.cc [new file with mode: 0644]
tests/fe/fe_enriched_06.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_07.cc [new file with mode: 0644]
tests/fe/fe_enriched_07.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_08.cc [new file with mode: 0644]
tests/fe/fe_enriched_08.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_08a.cc [new file with mode: 0644]
tests/fe/fe_enriched_08a.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_compare_to_fe_system.cc [new file with mode: 0644]
tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_compare_to_fe_system_2.cc [new file with mode: 0644]
tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_step-36.cc [new file with mode: 0644]
tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output [new file with mode: 0644]
tests/fe/fe_enriched_step-36b.cc [new file with mode: 0644]
tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/fe_enriched_1d.png b/doc/doxygen/images/fe_enriched_1d.png
new file mode 100644 (file)
index 0000000..6d1af31
Binary files /dev/null and b/doc/doxygen/images/fe_enriched_1d.png differ
diff --git a/doc/doxygen/images/fe_enriched_h-refinement.png b/doc/doxygen/images/fe_enriched_h-refinement.png
new file mode 100644 (file)
index 0000000..909cee3
Binary files /dev/null and b/doc/doxygen/images/fe_enriched_h-refinement.png differ
index 1562747f8c819618e836144abd32d9c7b815f5b9..36c1fa4fed4231afd7debfdc1b29934f5db24a37 100644 (file)
@@ -26,3 +26,11 @@ div.doxygen-generated-exception-message {
     padding: 5px;
     font-size: 95%;
 }
+
+div.image img[src="fe_enriched_1d.png"] {
+    width:500px;
+}
+
+div.image img[src="fe_enriched_h-refinement.png"] {
+    width:500px;
+}
index de2cc3909372ee3da623012dc13ad18dfc1397ae..9859b2b5d2d352dbdd86f1664b721ad86b05a244 100644 (file)
@@ -206,6 +206,12 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
+ <li> New: FE_Enriched finite element class implements the partition of unitity
+ method which allows to enrich the finite element space based on a priori
+ knowledge about solution.
+ <br>
+ (Denis Davydov, 2016/09/28)
+ </li>
 
  <li> Improved: The doxygen documentation now contains nicely formatted
  boxes containing the text message of each exception. Several messages
diff --git a/include/deal.II/fe/fe_enriched.h b/include/deal.II/fe/fe_enriched.h
new file mode 100644 (file)
index 0000000..b3131e4
--- /dev/null
@@ -0,0 +1,613 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__fe_enriched_h
+#define dealii__fe_enriched_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_CXX11
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/quadrature.h>
+
+#include <vector>
+#include <utility>
+#include <numeric>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Implementation of a partition of unity finite element method (PUM) by
+ * Babuska and Melenk which enriches a standard
+ * finite element with an enrichment function multiplied with another (usually
+ * linear) finite element:
+ * \f[
+ * U(\mathbf x) = \sum_i N_i(\mathbf x) U_i + \sum_j N_j(\mathbf x) \sum_k F_k(\mathbf x) U_{jk}
+ * \f]
+ * where $ N_i(\mathbf x) $ and $ N_j(\mathbf x) $ are the underlying finite elements (including
+ * the mapping from the isoparametric element to the real element); $ F_k(\mathbf x) $
+ * are the scalar enrichment functions in real space (e.g. $ 1/r $, $ \exp(-r) $, etc);
+ * $ U_i $ and $ U_{jk} $ are the standard and enriched DoFs. This allows to
+ * include in the finite element space a priori knowledge about the partial
+ * differential equation being solved which in turn improves the local
+ * approximation properties of the spaces. This can be useful for highly oscillatory
+ * solutions, problems with domain corners or on unbounded domains or sudden
+ * changes of boundary conditions. PUM method uses finite element spaces which
+ * satisfy the partition of unity property (e.g. FE_Q). Among other properties
+ * this makes the resulting space to reproduce enrichment functions exactly.
+ *
+ * The simplest constructor of this class takes two finite element objects and an
+ * enrichment function to be used. For example
+ *
+ * @code
+ * FE_Enriched<dim> fe(FE_Q<dim>(2),
+ *                     FE_Q<dim>(1),
+ *                     function)
+ * @endcode
+ *
+ * In this case, standard DoFs are distributed by <code>FE_Q<dim>(2)</code>,
+ * whereas enriched DoFs are coming from a single finite element
+ * <code>FE_Q<dim>(1)</code> used with a single enrichment function
+ * <code>function</code>. In this case, the total number of DoFs on the
+ * enriched element is the sum of DoFs from <code>FE_Q<dim>(2)</code> and
+ * <code>FE_Q<dim>(1)</code>.
+ *
+ * As an example of an enrichment function, consider $ \exp(-x) $, which
+ * leads to the following shape functions on the unit element:
+ * <table>
+ *   <tr>
+ *     <td align="center">
+ *       @image html fe_enriched_1d.png
+ *     </td>
+ *     <td align="center">
+ *       @image html fe_enriched_h-refinement.png
+ *     </td>
+ *   </tr>
+ *   <tr>
+ *     <td align="center">
+ *       1d element, base and enriched shape functions.
+ *     </td>
+ *     <td align="center">
+ *       enriched shape function corresponding to the central vertex.
+ *     </td>
+ *   </tr>
+ * </table>
+ *
+ * Note that evaluation of gradients (hessians) of the enriched shape functions
+ * or the finite element field requires evaluation of gradients (gradients and hessians)
+ * of the enrichment functions:
+ * @f{align*}
+ *   U(\mathbf x)
+ *     &= \sum_i N_i(\mathbf x) U_i
+ *     + \sum_{j,k} N_j(\mathbf x) F_k(\mathbf x) U_{jk} \\
+ *   \mathbf \nabla U(\mathbf x)
+ *     &= \sum_i \mathbf \nabla N_i(\mathbf x) U_i
+ *     + \sum_{j,k} \left[\mathbf \nabla N_j(\mathbf x) F_k(\mathbf x) +
+ *                        N_j(\mathbf x) \mathbf \nabla F_k(\mathbf x) \right] U_{jk} \\
+ *   \mathbf \nabla \mathbf \nabla U(\mathbf x)
+ *     &= \sum_i \mathbf \nabla \mathbf \nabla N_i(\mathbf x) U_i
+ *     + \sum_{j,k} \left[\mathbf \nabla \mathbf \nabla N_j(\mathbf x) F_k(\mathbf x) +
+ *                        \mathbf \nabla F_k(\mathbf x) \mathbf \nabla N_j(\mathbf x) +
+ *                        \mathbf \nabla N_j(\mathbf x) \mathbf \nabla F_k(\mathbf x) +
+ *                        N_j(\mathbf x) \mathbf \nabla \mathbf \nabla F_k(\mathbf x) \right] U_{jk}
+ * @f}
+ *
+ * <h3>Using enriched and non-enriched FEs together</h3>
+ *
+ * In most applications it is beneficial to introduce enrichments only in
+ * some part of the domain (e.g. around a crack tip) and use standard FE (e.g. FE_Q)
+ * elsewhere.
+ * This can be achieved by using the hp finite element framework in deal.II
+ * that allows for the use of different elements on different cells. To make
+ * the resulting space $C^0$ continuous, it is then necessary for the DoFHandler
+ * class and DoFTools::make_hanging_node_constraints() function to be able to
+ * figure out what to do at the interface between enriched and non-enriched
+ * cells. Specifically, we want the degrees of freedom corresponding to
+ * enriched shape functions to be zero at these interfaces. These classes and
+ * functions can not to do this automatically, but the effect can be achieved
+ * by using not just a regular FE_Q on cells without enrichment, but to wrap
+ * the FE_Q into an FE_Enriched object <i>without actually enriching it</i>.
+ * This can be done as follows:
+ * @code
+ *   FE_Enriched<dim> fe_non_enriched(FE_Q<dim>(1));
+ * @endcode
+ * This constructor is equivalent to calling
+ * @code
+ *   FE_Enriched<dim> fe_non_enriched(FE_Q<dim>(1),
+ *                                    FE_Nothing<dim>(1,true),
+ *                                    NULL);
+ * @endcode
+ * and will result in the correct constraints for enriched DoFs attributed to
+ * support points on the interface between the two regions.
+ *
+ * <h3>References</h3>
+ *
+ * When using this class, please cite
+ * @code{.bib}
+ @Article{Davydov2016,
+   Title                    = {On the h-adaptive PUM and hp-adaptive FEM approaches applied to PDEs in quantum mechanics.},
+   Author                   = {Davydov, D and Gerasimov, T and Pelteret, J.-P. and Steinmann, P.},
+   Journal                  = {Journal of Computational Physics. Submitted},
+   Year                     = {2016},
+ }
+ * @endcode
+ * The PUM was introduced in
+ * @code{.bib}
+ @Article{Melenk1996,
+  Title                    = {The partition of unity finite element method: Basic theory and applications },
+  Author                   = {Melenk, J.M. and Babu\v{s}ka, I.},
+  Journal                  = {Computer Methods in Applied Mechanics and Engineering},
+  Year                     = {1996},
+  Number                   = {1--4},
+  Pages                    = {289 -- 314},
+  Volume                   = {139},
+}
+@Article{Babuska1997,
+  Title                    = {The partition of unity method},
+  Author                   = {Babu\v{s}ka, I. and Melenk, J. M.},
+  Journal                  = {International Journal for Numerical Methods in Engineering},
+  Year                     = {1997},
+  Number                   = {4},
+  Pages                    = {727--758},
+  Volume                   = {40},
+}
+ * @endcode
+ *
+ * <h3>Implementation</h3>
+ *
+ * The implementation of the class is based on FESystem which is aggregated as
+ * a private member. The simplest constructor <code> FE_Enriched<dim> fe(FE_Q<dim>(2), FE_Q<dim>(1),function)</code>
+ * will internally initialize FESystem as
+ *
+ * @code
+ * FESystem<dim> fe_system(FE_Q<dim>(2),1,
+ *                         FE_Q<dim>(1),1);
+ * @endcode
+ *
+ * Note that it would not be wise to have this class derived
+ * from FESystem as the latter concatenates the given elements into different
+ * components of a vector element, whereas the current class combines the given
+ * elements into the same components. For instance, if two scalar elements are
+ * given, the resulting element will be scalar rather than have two components
+ * when doing the same with an FESystem.
+ *
+ * The ordering of the shape function, @p interface_constrains, the @p prolongation (embedding)
+ * and the @p restriction matrices are taken from the FESystem class.
+ *
+ * @note This class is only available when deal.II is compiled with C++11.
+ *
+ * @ingroup fe
+ *
+ * @author Denis Davydov, 2016.
+ */
+template <int dim, int spacedim=dim>
+class FE_Enriched : public FiniteElement<dim, spacedim>
+{
+public:
+
+  /**
+   * Constructor which takes base FiniteElement @p fe_base and the enrichment
+   * FiniteElement @p fe_enriched which will be multiplied by the @p enrichment_function.
+   *
+   * In case @p fe_enriched is other than FE_Nothing, the lifetime of the
+   * @p enrichment_function must be at least as long as the FE_Enriched object.
+   */
+  FE_Enriched (const FiniteElement<dim,spacedim> &fe_base,
+               const FiniteElement<dim,spacedim> &fe_enriched,
+               const Function<spacedim>      *enrichment_function);
+
+  /**
+   * Constructor which only wraps the base FE @p fe_base.
+   * As for the enriched finite element space, FE_Nothing is used.
+   * Continuity constraints will be automatically generated when
+   * this non-enriched element is used in conjunction with enriched finite element
+   * within the hp::DoFHandler.
+   *
+   * See the discussion in the class documentation on how to use this element
+   * in the context of hp finite element methods.
+   */
+  FE_Enriched (const FiniteElement<dim,spacedim> &fe_base);
+
+  /**
+   * Constructor which takes pointer to the base FiniteElement @p fe_base and
+   * a vector of enriched FiniteElement's @p fe_enriched . @p fe_enriched[i]
+   * finite element will be enriched with functions in @p functions[i].
+   *
+   * This is the most general public constructor which also allows to have
+   * different enrichment functions in different disjoint parts of the domain.
+   * To that end the last argument provides an association of cell iterator
+   * to a Function. This is done to simplify the usage of this class when the
+   * number of disjoint domains with different functions is more than a few.
+   * Otherwise one would have to use different instance of this class for each
+   * disjoint enriched domain.
+   *
+   * @note When using the same finite element for enrichment with N
+   * different functions, it is advised to have the second argument of size 1
+   * and the last argument of size 1 x N. The same can be achieved by providing
+   * N equivalent enrichment elements while keeping the last argument of size
+   * N x 1. However this will be much more computationally expensive.
+   *
+   * @note When using different enrichment functions on disjoint domains, no
+   * checks are done by this class that the domains are actually disjoint.
+   */
+  FE_Enriched (const FiniteElement<dim,spacedim> *fe_base,
+               const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,
+               const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions);
+
+private:
+  /**
+   * The most general private constructor. The first two input parameters are
+   * consistent with those in FESystem. It is used internally only with
+   * <code>multiplicities[0]=1</code>, which is a logical requirement for this finite element.
+   */
+  FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
+               const std::vector< unsigned int > &multiplicities,
+               const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions);
+public:
+
+  /**
+   * @p clone function instead of a copy constructor.
+   *
+   * This function is needed by the constructors of @p FESystem.
+   */
+  virtual FiniteElement<dim,spacedim> *clone() const;
+
+  virtual
+  UpdateFlags
+  requires_update_flags (const UpdateFlags update_flags) const;
+
+  /**
+   * Return a string that identifies a finite element.
+   */
+  virtual std::string get_name () const;
+
+  /**
+   * Access to a composing element. The index needs to be smaller than the
+   * number of base elements. In the context of this class, the number of
+   * base elements is always more than one: a non-enriched element plus an
+   * element to be enriched, which could be FE_Nothing.
+   */
+  virtual const FiniteElement<dim,spacedim> &
+  base_element (const unsigned int index) const;
+
+  /**
+   * Return the value of the @p ith shape function at the point @p p. @p p is a
+   * point on the reference element.
+   *
+   * This function returns meaningful values only for non-enriched element as
+   * real-space enrichment requires evaluation of the function at the point in
+   * real-space.
+   */
+  virtual double shape_value(const unsigned int      i,
+                             const Point< dim >     &p) const;
+
+
+  /**
+   * @name Functions to support hp
+   * @{
+   */
+
+  /**
+   * Return whether this element implements hp constraints.
+   *
+   * This function returns @p true if and only if all its base elements return @p true
+   * for this function.
+   */
+  virtual bool hp_constraints_are_implemented () const;
+
+  /**
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
+   *
+   * Base elements of this element will have to implement this function. They
+   * may only provide interpolation matrices for certain source finite
+   * elements, for example those from the same family. If they don't implement
+   * interpolation from a given element, then they must throw an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented, which
+   * will get propagated out from this element.
+   */
+  virtual void
+  get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                                 FullMatrix<double>       &matrix) const;
+
+  /**
+   * Return the matrix interpolating from a face of of one element to the
+   * subface of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
+   *
+   * Base elements of this element will have to implement this function. They
+   * may only provide interpolation matrices for certain source finite
+   * elements, for example those from the same family. If they don't implement
+   * interpolation from a given element, then they must throw an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented, which
+   * will get propagated out from this element.
+   */
+  virtual void
+  get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                                    const unsigned int        subface,
+                                    FullMatrix<double>       &matrix) const;
+
+  /**
+   * If, on a vertex, several finite elements are active, the hp code first
+   * assigns the degrees of freedom of each of these FEs different global
+   * indices. It then calls this function to find out which of them should get
+   * identical values, and consequently can receive the same global DoF index.
+   * This function therefore returns a list of identities between DoFs of the
+   * present finite element object with the DoFs of @p fe_other, which is a
+   * reference to a finite element object representing one of the other finite
+   * elements active on this particular vertex. The function computes which of
+   * the degrees of freedom of the two finite element objects are equivalent,
+   * both numbered between zero and the corresponding value of dofs_per_vertex
+   * of the two finite elements. The first index of each pair denotes one of
+   * the vertex dofs of the present element, whereas the second is the
+   * corresponding index of the other finite element.
+   */
+  virtual
+  std::vector<std::pair<unsigned int, unsigned int> >
+  hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+  /**
+   * Same as hp_vertex_dof_indices(), except that the function treats degrees
+   * of freedom on lines.
+   */
+  virtual
+  std::vector<std::pair<unsigned int, unsigned int> >
+  hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+  /**
+   * Same as hp_vertex_dof_indices(), except that the function treats degrees
+   * of freedom on quads.
+   */
+  virtual
+  std::vector<std::pair<unsigned int, unsigned int> >
+  hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+  /**
+   * Return whether this element dominates the one given as argument when they
+   * meet at a common face, whether it is the other way around, whether
+   * neither dominates, or if either could dominate.
+   *
+   * For a definition of domination, see FiniteElementBase::Domination and in
+   * particular the
+   * @ref hp_paper "hp paper".
+   */
+  virtual
+  FiniteElementDomination::Domination
+  compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+  //@}
+
+
+  /**
+   * Return enrichment functions
+   */
+  const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > >
+  get_enrichments() const;
+
+  /**
+   * Return the underlying FESystem object.
+   */
+  const FESystem<dim,spacedim> &
+  get_fe_system() const;
+
+protected:
+
+  /**
+   * A class to hold internal data needed for evaluation of this FE at quadrature points.
+   */
+  class InternalData : public FiniteElement<dim,spacedim>::InternalDataBase
+  {
+  public:
+    /**
+     * For each Finite Element (base number) and each enrichment function (base_index)
+     * this struct will contain values, gradients and hessians of the enrichment functions.
+     */
+    struct EnrichmentValues
+    {
+      std::vector<double> values;
+      std::vector<Tensor<1,spacedim> > gradients;
+      std::vector<SymmetricTensor<2, spacedim> > hessians;
+    };
+
+    /**
+     * Constructor. Is used inside setup_data to wrap FESystem's internal
+     * data object. The former is called from get_data, get_subface_data and
+     * get_face_data which FE_Enriched has to implement.
+     *
+     * Since FESystem::get_data(), FESystem::get_face_data() and FESystem::get_subface_data()
+     * just create an object and return a pointer to it (i.e. they don't retain
+     * ownership), we store the cast result in a std::unique_ptr to indicate
+     * that InternalData owns the object.
+     */
+    InternalData ( std::unique_ptr<typename FESystem<dim,spacedim>::InternalData> fesystem_data);
+
+    /**
+     * Gives read-access to the pointer to a @p InternalData of the @p
+     * <code>base_no</code>th base element of FESystem's data.
+     */
+    typename FiniteElement<dim,spacedim>::InternalDataBase &
+    get_fe_data (const unsigned int base_no) const;
+
+    /**
+     * Gives read-access to the pointer to an object into which the
+     * <code>base_no</code>th base element will write its output when calling
+     * FiniteElement::fill_fe_values() and similar functions.
+     */
+    internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
+    get_fe_output_object (const unsigned int base_no) const;
+
+    /**
+     * Aggregate FESystem's internal data. It is used every time
+     * we call FESystem's fill_fe_values() and alike.
+     */
+    std::unique_ptr<typename FESystem<dim,spacedim>::InternalData> fesystem_data;
+
+    /**
+     * For each FE used in enrichment (base number <code>i</code>) and each enrichment function
+     * (base multiplicity <code>j</code>), <code>enrichment_values[i][j]</code> will be used to store
+     * possibly requested values, gradients and hessians of enrichment function <code>j</code>.
+     *
+     * The variable is made mutable as InternalData's provided to fill_fe_values and alike
+     * are const.
+     *
+     * @note We do not want to store this information in the finite element object itself,
+     * because this would mean that (i) only one FEValues object could use a finite element object at a time,
+     * and (ii) that these objects could not be used in a multithreaded context.
+     */
+    mutable std::vector<std::vector<EnrichmentValues> > enrichment;
+  };
+
+  /**
+   * For each finite element @p i used in enrichment and each enrichment function
+   * @p j associated with it (essentially its multiplicity),
+   * @p base_no_mult_local_enriched_dofs[i][j] contains the associated local DoFs
+   * on the FE_Enriched finite element.
+   */
+  std::vector<std::vector<std::vector<unsigned int> > > base_no_mult_local_enriched_dofs;
+
+  /**
+   * Enrichment functions.
+   * The size of the first vector is the same as the number of FiniteElement spaces used
+   * with enrichment. Whereas the size of the inner vector corresponds to the number
+   * of enrichment functions associated with a single FiniteElement.
+   */
+  const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > enrichments;
+
+  /**
+   * Auxiliary variable used to distinguish between the case when we do enrichment
+   * and when the class simply wraps another FiniteElement.
+   *
+   * This variable is initialized in the constructor by looping over a vector of
+   * enrichment elements and checking if all of them are FE_Nothing. If this is
+   * the case, then the value is set to <code>false</code>, otherwise it is
+   * <code>true</code>.
+   */
+  const bool is_enriched;
+
+  /**
+   * Auxiliary function called from get_data, get_face_data and
+   * get_subface_data. It take internal data of FESystem object in @p fes_data
+   * and the quadrature rule @p qudrature.
+   *
+   * This function essentially take the internal data from an instance of
+   * FESystem class and wraps it into our own InternalData class which
+   * additionally has objects to hold values/gradients/hessians of
+   * enrichment functions at each quadrature point depending on @p flags.
+   */
+  template <int dim_1>
+  typename FiniteElement<dim,spacedim>::InternalDataBase *
+  setup_data (std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fes_data,
+              const UpdateFlags      flags,
+              const Quadrature<dim_1> &quadrature) const;
+
+  /**
+   * Prepare internal data structures and fill in values independent of the
+   * cell. Returns a pointer to an object of which the caller of this function
+   * (FEValues) then has to assume ownership (which includes destruction when it is no
+   * more needed).
+   */
+  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
+  get_data (const UpdateFlags      flags,
+            const Mapping<dim,spacedim>    &mapping,
+            const Quadrature<dim> &quadrature,
+            dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const;
+
+  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
+  get_face_data (const UpdateFlags      update_flags,
+                 const Mapping<dim,spacedim>    &mapping,
+                 const Quadrature<dim-1> &quadrature,
+                 dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim >        &output_data) const;
+
+  virtual typename FiniteElement<dim,spacedim>::InternalDataBase *
+  get_subface_data (const UpdateFlags      update_flags,
+                    const Mapping<dim,spacedim>    &mapping,
+                    const Quadrature<dim-1> &quadrature,
+                    dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const;
+
+  virtual
+  void fill_fe_values (const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+                       const CellSimilarity::Similarity cell_similarity,
+                       const Quadrature<dim> &quadrature,
+                       const Mapping<dim, spacedim> &mapping,
+                       const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
+                       const dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
+                       const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+                       dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data
+                      ) const;
+
+  virtual
+  void
+  fill_fe_face_values ( const typename Triangulation<dim, spacedim>::cell_iterator   &cell,
+                        const unsigned int face_no,
+                        const Quadrature<dim-1> &quadrature,
+                        const Mapping<dim, spacedim> &mapping,
+                        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
+                        const dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
+                        const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+                        dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data
+                      ) const;
+
+  virtual
+  void
+  fill_fe_subface_values (const typename Triangulation< dim, spacedim >::cell_iterator &cell,
+                          const unsigned int face_no,
+                          const unsigned int sub_no,
+                          const Quadrature<dim-1> &quadrature,
+                          const Mapping<dim, spacedim> &mapping,
+                          const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,
+                          const dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
+                          const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+                          dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data
+                         ) const;
+
+private:
+  /**
+   * This function sets up the index table for the system as well as @p
+   * restriction and @p prolongation matrices.
+   */
+  void initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+                   const std::vector<unsigned int> &multiplicities);
+
+  /**
+   * The underlying FESystem object.
+   */
+  FESystem<dim,spacedim> fe_system;
+
+  /**
+   * After calling fill_fe_(face/subface_)values this function
+   * implements the chain rule to multiply stored shape values/gradient/hessians
+   * by those of enrichment function evaluated at quadrature points.
+   */
+  template <int dim_1>
+  void
+  multiply_by_enrichment (const Quadrature<dim_1>                                       &quadrature,
+                          const InternalData                                            &fe_data,
+                          const internal::FEValues::MappingRelatedData<dim,spacedim>    &mapping_data,
+                          const typename Triangulation< dim, spacedim >::cell_iterator  &cell,
+                          internal::FEValues::FiniteElementRelatedData<dim,spacedim>    &output_data) const;
+};
+
+//}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // CXX11
+
+#endif // dealii__fe_enriched_h
+
index 1f87be95e3f4db2428ab5e90859970667b919725..89b017bd1ecbb812e450ae106db45f59a1a41dce 100644 (file)
@@ -28,6 +28,8 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+template <int dim, int spacedim> class FE_Enriched;
+
 
 /**
  * This class provides an interface to group several elements together into
@@ -959,7 +961,6 @@ private:
       unsigned int> >
       base_elements;
 
-
   /**
    * Initialize the @p unit_support_points field of the FiniteElement class.
    * Called from the constructor.
@@ -1071,6 +1072,8 @@ private:
    * Mutex for protecting initialization of restriction and embedding matrix.
    */
   mutable Threads::Mutex mutex;
+
+  friend class FE_Enriched<dim,spacedim>;
 };
 
 
index 0e8f549797cc141fa10f5a1c0cf4d98d90e38372..f0a674c08ce10bf783781fc62ccb3479fc20f14f 100644 (file)
@@ -44,6 +44,7 @@ SET(_src
   fe_raviart_thomas_nodal.cc
   fe_series.cc
   fe_system.cc
+  fe_enriched.cc
   fe_tools.cc
   fe_tools_interpolate.cc
   fe_trace.cc
@@ -86,6 +87,7 @@ SET(_inst
   fe_raviart_thomas.inst.in
   fe_raviart_thomas_nodal.inst.in
   fe_system.inst.in
+  fe_enriched.inst.in
   fe_tools.inst.in
   fe_tools_interpolate.inst.in
   fe_trace.inst.in
diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc
new file mode 100644 (file)
index 0000000..7486d6b
--- /dev/null
@@ -0,0 +1,926 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+#include <deal.II/fe/fe_enriched.h>
+#ifdef DEAL_II_WITH_CXX11
+
+#include <deal.II/fe/fe_tools.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+  /**
+   * Auxiliary function to create multiplicity vector from input enrichment functions.
+   */
+  template <typename T>
+  std::vector<unsigned int>
+  build_multiplicities(const std::vector<std::vector<T > > &functions )
+  {
+    std::vector<unsigned int> multiplicities;
+    multiplicities.push_back(1); // the first one is non-enriched FE
+    for (unsigned int i = 0; i < functions.size(); i++)
+      multiplicities.push_back(functions[i].size());
+
+    return multiplicities;
+  }
+
+
+  /**
+   * Auxiliary function to build FiniteElement's vector
+   */
+  template <int dim, int spacedim>
+  std::vector< const FiniteElement< dim, spacedim > * >
+  build_fes(const FiniteElement<dim,spacedim> *fe_base,
+            const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)
+  {
+    std::vector< const FiniteElement< dim, spacedim > * > fes;
+    fes.push_back(fe_base);
+    for (unsigned int i = 0; i < fe_enriched.size(); i++)
+      fes.push_back(fe_enriched[i]);
+
+    return fes;
+  }
+
+
+  /**
+   * Auxiliary function which check consistency of the input parameters.
+   * Returns true if everything is ok.
+   */
+  template <int dim, int spacedim>
+  bool
+  consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
+                     const std::vector< unsigned int > &multiplicities,
+                     const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
+  {
+    AssertThrow(fes.size() > 0,
+                ExcMessage("FEs size should be >=1"));
+    AssertThrow(fes.size() == multiplicities.size(),
+                ExcMessage("FEs and multiplicities should have the same size"));
+
+    AssertThrow (functions.size() == fes.size() - 1,
+                 ExcDimensionMismatch(functions.size(), fes.size()-1));
+
+    AssertThrow(multiplicities[0] == 1,
+                ExcMessage("First multiplicity should be 1"));
+
+    const unsigned int n_comp_base = fes[0]->n_components();
+
+    // start from fe=1 as 0th is always non-enriched FE.
+    for (unsigned int fe=1; fe < fes.size(); fe++)
+      {
+        const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);
+        if (fe_nothing)
+          AssertThrow (fe_nothing->is_dominating(),
+                       ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched"));
+
+        AssertThrow (fes[fe]->n_components() == n_comp_base,
+                     ExcMessage("All elements must have the same number of components"));
+      }
+    return true;
+  }
+
+
+  /**
+   * Auxiliary function which determines whether the FiniteElement will be enriched.
+   */
+  template <int dim, int spacedim>
+  bool
+  check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)
+  {
+    // start from fe=1 as 0th is always non-enriched FE.
+    for (unsigned int fe=1; fe < fes.size(); fe++)
+      if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == 0)
+        // this is not FE_Nothing => there will be enrichment
+        return true;
+
+    return false;
+  }
+}
+
+
+template <int dim, int spacedim>
+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> &fe_base)
+  :
+  FE_Enriched<dim,spacedim>(fe_base,
+                            FE_Nothing<dim,spacedim>(fe_base.n_components(),true),
+                            NULL)
+{
+}
+
+
+template <int dim, int spacedim>
+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> &fe_base,
+                                        const FiniteElement<dim,spacedim> &fe_enriched,
+                                        const Function<spacedim>      *enrichment_function)
+  :
+  FE_Enriched<dim,spacedim>
+  (&fe_base, { &fe_enriched },
+{
+  {
+    [=] (const typename Triangulation<dim, spacedim>::cell_iterator &) -> const Function<spacedim> *
+    {
+      return enrichment_function;
+    }
+  }
+})
+{
+}
+
+
+template <int dim, int spacedim>
+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> *fe_base,
+                                        const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,
+                                        const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
+  :
+  FE_Enriched<dim,spacedim> (build_fes(fe_base,fe_enriched),
+                             build_multiplicities(functions),
+                             functions)
+{
+}
+
+
+template <int dim, int spacedim>
+FE_Enriched<dim,spacedim>::FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
+                                        const std::vector< unsigned int > &multiplicities,
+                                        const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
+  :
+  FiniteElement<dim,spacedim> (FETools::Compositing::multiply_dof_numbers(fes,multiplicities,false),
+                               FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities),
+                               FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)),
+  enrichments(functions),
+  is_enriched(check_if_enriched(fes)),
+  fe_system(fes,multiplicities)
+{
+  // descriptive error are thrown within the function.
+  Assert(consistency_check(fes,multiplicities,functions),
+         ExcInternalError());
+
+  initialize(fes, multiplicities);
+
+  // resize to be consistent with all FEs used to construct the FE_Enriched,
+  // even though we will never use the 0th element.
+  base_no_mult_local_enriched_dofs.resize(fes.size());
+  for (unsigned int fe=1; fe < fes.size(); fe++)
+    base_no_mult_local_enriched_dofs[fe].resize(multiplicities[fe]);
+
+  Assert (base_no_mult_local_enriched_dofs.size() == this->n_base_elements(),
+          ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),
+                               this->n_base_elements()));
+
+  // build the map: (base_no, base_m) -> vector of local element DoFs
+  for (unsigned int system_index=0; system_index<this->dofs_per_cell;
+       ++system_index)
+    {
+      const unsigned int base_no = this->system_to_base_table[system_index].first.first;
+      if (base_no == 0) // 0th is always non-enriched FE
+        continue;
+
+      const unsigned int base_m  = this->system_to_base_table[system_index].first.second;
+      const unsigned int
+      base_index = this->system_to_base_table[system_index].second;
+
+      Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(),
+              ExcMessage("Size mismatch for base_no_mult_local_enriched_dofs: "
+                         "base_index = " + std::to_string(base_index) +
+                         "; base_no = " + std::to_string(base_no) +
+                         "; base_m = " + std::to_string(base_m) +
+                         "; system_index = " + std::to_string(system_index)));
+
+      Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(),
+              ExcDimensionMismatch(base_m,
+                                   base_no_mult_local_enriched_dofs[base_no].size()));
+
+      base_no_mult_local_enriched_dofs[base_no][base_m].push_back(system_index);
+    }
+
+  // make sure that local_enriched_dofs.size() is correct, that is equals to DoFs
+  // per cell of the corresponding FE.
+  for (unsigned int base_no = 1; base_no < base_no_mult_local_enriched_dofs.size(); base_no++)
+    {
+      for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)
+        Assert ( base_no_mult_local_enriched_dofs[base_no][m].size() == fes[base_no]->dofs_per_cell,
+                 ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no][m].size(),
+                                      fes[base_no]->dofs_per_cell));
+    }
+}
+
+
+template <int dim, int spacedim>
+const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > >
+FE_Enriched<dim,spacedim>::get_enrichments() const
+{
+  return enrichments;
+}
+
+
+template <int dim, int spacedim>
+double
+FE_Enriched<dim,spacedim>::shape_value(const unsigned int   i,
+                                       const Point< dim > &p) const
+{
+  Assert(!is_enriched, ExcMessage("For enriched finite elements shape_value() can not be defined on the reference element."));
+  return fe_system.shape_value(i,p);
+}
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_Enriched<dim,spacedim>::clone() const
+{
+  std::vector< const FiniteElement< dim, spacedim > * > fes;
+  std::vector< unsigned int > multiplicities;
+
+  for (unsigned int i=0; i<this->n_base_elements(); i++)
+    {
+      fes.push_back( & base_element(i) );
+      multiplicities.push_back(this->element_multiplicity(i) );
+    }
+
+  return new FE_Enriched<dim,spacedim>(fes, multiplicities, get_enrichments());
+}
+
+
+template <int dim, int spacedim>
+UpdateFlags
+FE_Enriched<dim,spacedim>::requires_update_flags (const UpdateFlags flags) const
+{
+  UpdateFlags out = fe_system.requires_update_flags(flags);
+
+  if (is_enriched)
+    {
+      // if we ask for values or gradients, then we would need quadrature points
+      if (flags & (update_values | update_gradients))
+        out |= update_quadrature_points;
+
+      // if need gradients, add update_values due to product rule
+      if (out & update_gradients)
+        out |= update_values;
+    }
+
+  Assert (!(flags & update_3rd_derivatives),
+          ExcNotImplemented());
+
+  return out;
+}
+
+
+template <int dim, int spacedim>
+template <int dim_1>
+typename FiniteElement<dim,spacedim>::InternalDataBase *
+FE_Enriched<dim,spacedim>::setup_data (std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fes_data,
+                                       const UpdateFlags      flags,
+                                       const Quadrature<dim_1> &quadrature) const
+{
+  Assert ((dynamic_cast<typename FESystem<dim,spacedim>::InternalData *> (fes_data.get()) != NULL),
+          ExcInternalError());
+  typename FESystem<dim,spacedim>::InternalData *data_fesystem =
+    static_cast<typename FESystem<dim,spacedim>::InternalData *> (fes_data.get());
+
+  // FESystem::InternalData will be aggregated (owned) by
+  // our InternalData.
+  fes_data.release();
+  InternalData *data = new InternalData(std::unique_ptr<typename FESystem<dim,spacedim>::InternalData>(data_fesystem));
+
+  // copy update_each from FESystem data:
+  data->update_each = data_fesystem->update_each;
+
+  // resize cache array according to requested flags
+  data->enrichment.resize(this->n_base_elements());
+
+  const unsigned int n_q_points = quadrature.size();
+
+  for (unsigned int base=0; base < this->n_base_elements(); ++base)
+    {
+      data->enrichment[base].resize(this->element_multiplicity(base));
+      for (unsigned int m = 0; m < this->element_multiplicity(base); ++m)
+        {
+          if (flags & update_values)
+            data->enrichment[base][m].values.resize(n_q_points);
+
+          if (flags & update_gradients)
+            data->enrichment[base][m].gradients.resize(n_q_points);
+
+          if (flags & update_hessians)
+            data->enrichment[base][m].hessians.resize(n_q_points);
+        }
+    }
+
+  return data;
+}
+
+
+template <int dim, int spacedim>
+typename FiniteElement<dim,spacedim>::InternalDataBase *
+FE_Enriched<dim,spacedim>::get_face_data (const UpdateFlags      update_flags,
+                                          const Mapping<dim,spacedim>    &mapping,
+                                          const Quadrature<dim-1> &quadrature,
+                                          internal::FEValues::FiniteElementRelatedData< dim, spacedim >        &output_data) const
+{
+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_face_data(update_flags,mapping,quadrature,output_data)),
+                    update_flags,
+                    quadrature);
+}
+
+
+template <int dim, int spacedim>
+typename FiniteElement<dim,spacedim>::InternalDataBase *
+FE_Enriched<dim,spacedim>::get_subface_data (const UpdateFlags      update_flags,
+                                             const Mapping<dim,spacedim>    &mapping,
+                                             const Quadrature<dim-1> &quadrature,
+                                             dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const
+{
+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_subface_data(update_flags,mapping,quadrature,output_data)),
+                    update_flags,
+                    quadrature);
+}
+
+
+template <int dim, int spacedim>
+typename FiniteElement<dim,spacedim>::InternalDataBase *
+FE_Enriched<dim,spacedim>::get_data (const UpdateFlags      flags,
+                                     const Mapping<dim,spacedim>    &mapping,
+                                     const Quadrature<dim> &quadrature,
+                                     internal::FEValues::FiniteElementRelatedData< dim, spacedim >   &output_data) const
+{
+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_data(flags,mapping,quadrature,output_data)),
+                    flags,
+                    quadrature);
+}
+
+
+template <int dim, int spacedim>
+void FE_Enriched<dim,spacedim>::initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+                                            const std::vector<unsigned int> &multiplicities)
+{
+  Assert (fes.size() == multiplicities.size(),
+          ExcDimensionMismatch (fes.size(), multiplicities.size()) );
+
+  // Note that we need to skip every fe with multiplicity 0 in the following block of code
+  this->base_to_block_indices.reinit(0, 0);
+
+  for (unsigned int i=0; i<fes.size(); i++)
+    if (multiplicities[i]>0)
+      this->base_to_block_indices.push_back( multiplicities[i] );
+
+  {
+    // If the system is not primitive, these have not been initialized by
+    // FiniteElement
+    this->system_to_component_table.resize(this->dofs_per_cell);
+    this->face_system_to_component_table.resize(this->dofs_per_face);
+
+    FETools::Compositing::build_cell_tables(this->system_to_base_table,
+                                            this->system_to_component_table,
+                                            this->component_to_base_table,
+                                            *this,
+                                            false);
+
+    FETools::Compositing::build_face_tables(this->face_system_to_base_table,
+                                            this->face_system_to_component_table,
+                                            *this,
+                                            false);
+  }
+
+  // restriction and prolongation matrices are built on demand
+
+  // now set up the interface constraints for h-refinement.
+  // take them from fe_system:
+  this->interface_constraints = fe_system.interface_constraints;
+
+  // if we just wrap another FE (i.e. use FE_Nothing as a second FE)
+  // then it makes sense to have support points.
+  // However, functions like interpolate_boundary_values() need all FEs inside
+  // FECollection to be able to provide support points irrespectively whether
+  // this FE sits on the boundary or not. Thus for moment just copy support
+  // points from fe system:
+  {
+    this->unit_support_points      = fe_system.unit_support_points;
+    this->unit_face_support_points = fe_system.unit_face_support_points;
+  }
+
+  // take adjust_quad_dof_index_for_face_orientation_table from FESystem:
+  {
+    this->adjust_line_dof_index_for_line_orientation_table = fe_system.adjust_line_dof_index_for_line_orientation_table;
+  }
+}
+
+
+template <int dim, int spacedim>
+std::string
+FE_Enriched<dim,spacedim>::get_name () const
+{
+  std::ostringstream namebuf;
+
+  namebuf << "FE_Enriched<"
+          << Utilities::dim_string(dim,spacedim)
+          << ">[";
+  for (unsigned int i=0; i< this->n_base_elements(); ++i)
+    {
+      namebuf << base_element(i).get_name();
+      if (this->element_multiplicity(i) != 1)
+        namebuf << '^' << this->element_multiplicity(i);
+      if (i != this->n_base_elements()-1)
+        namebuf << '-';
+    }
+  namebuf << ']';
+
+  return namebuf.str();
+}
+
+
+template <int dim, int spacedim>
+const FiniteElement<dim,spacedim> &
+FE_Enriched<dim,spacedim>::base_element (const unsigned int index) const
+{
+  return fe_system.base_element(index);
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Enriched<dim,spacedim>::fill_fe_values (const typename Triangulation< dim, spacedim >::cell_iterator &cell,
+                                           const CellSimilarity::Similarity cell_similarity,
+                                           const Quadrature< dim > &quadrature,
+                                           const Mapping< dim, spacedim > &mapping,
+                                           const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,
+                                           const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,
+                                           const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+                                           internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data
+                                          ) const
+{
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,
+          ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  // call FESystem's method to fill everything without enrichment function
+  fe_system.fill_fe_values(cell,
+                           cell_similarity,
+                           quadrature,
+                           mapping,
+                           mapping_internal,
+                           mapping_data,
+                           *fe_data.fesystem_data,
+                           output_data);
+
+  if (is_enriched)
+    multiply_by_enrichment(quadrature,
+                           fe_data,
+                           mapping_data,
+                           cell,
+                           output_data);
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Enriched<dim,spacedim>::fill_fe_face_values
+(const typename Triangulation< dim, spacedim >::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature< dim-1 > &quadrature,
+ const Mapping< dim, spacedim > &mapping,
+ const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,
+ const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,
+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+ internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data
+) const
+{
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,
+          ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  // call FESystem's method to fill everything without enrichment function
+  fe_system.fill_fe_face_values(cell,
+                                face_no,
+                                quadrature,
+                                mapping,
+                                mapping_internal,
+                                mapping_data,
+                                *fe_data.fesystem_data,
+                                output_data);
+
+  if (is_enriched)
+    multiply_by_enrichment(quadrature,
+                           fe_data,
+                           mapping_data,
+                           cell,
+                           output_data);
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Enriched<dim,spacedim>::fill_fe_subface_values
+(const typename Triangulation< dim, spacedim >::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature< dim-1 > &quadrature,
+ const Mapping< dim, spacedim > &mapping,
+ const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,
+ const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,
+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
+ internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data
+) const
+{
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,
+          ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  // call FESystem's method to fill everything without enrichment function
+  fe_system.fill_fe_subface_values(cell,
+                                   face_no,
+                                   sub_no,
+                                   quadrature,
+                                   mapping,
+                                   mapping_internal,
+                                   mapping_data,
+                                   *fe_data.fesystem_data,
+                                   output_data);
+
+  if (is_enriched)
+    multiply_by_enrichment(quadrature,
+                           fe_data,
+                           mapping_data,
+                           cell,
+                           output_data);
+
+}
+
+
+template <int dim, int spacedim>
+template <int dim_1>
+void
+FE_Enriched<dim,spacedim>::multiply_by_enrichment
+(const Quadrature<dim_1> &quadrature,
+ const InternalData &fe_data,
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ const typename Triangulation< dim, spacedim >::cell_iterator &cell,
+ internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
+{
+  // mapping_data will contain quadrature points on the real element.
+  // fe_internal is needed to get update flags
+  // finally, output_data should store all the results we need.
+
+  // Either dim_1==dim
+  // (fill_fe_values) or dim_1==dim-1
+  // (fill_fe_(sub)face_values)
+  Assert(dim_1==dim || dim_1==dim-1, ExcInternalError());
+  const UpdateFlags flags = fe_data.update_each;
+
+  const unsigned int n_q_points  = quadrature.size();
+
+  // First, populate output_data object (that shall hold everything requested such
+  // as shape value, gradients, hessians, etc) from each base element.
+  // That is almost identical to FESystem::compute_fill_one_base(),
+  // the difference being that we do it irrespectively of cell_similarity
+  // and use base_fe_data.update_flags
+
+  // TODO: do we need it only for dim_1 == dim (i.e. fill_fe_values)?
+  if (dim_1 == dim)
+    for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
+      {
+        const FiniteElement<dim,spacedim> &
+        base_fe      = base_element(base_no);
+        typename FiniteElement<dim,spacedim>::InternalDataBase &
+        base_fe_data = fe_data.get_fe_data(base_no);
+        internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
+        base_data    = fe_data.get_fe_output_object(base_no);
+
+        const UpdateFlags base_flags = base_fe_data.update_each;
+
+        for (unsigned int system_index=0; system_index<this->dofs_per_cell;
+             ++system_index)
+          if (this->system_to_base_table[system_index].first.first == base_no)
+            {
+              const unsigned int
+              base_index = this->system_to_base_table[system_index].second;
+              Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());
+
+              // now copy. if the shape function is primitive, then there
+              // is only one value to be copied, but for non-primitive
+              // elements, there might be more values to be copied
+              //
+              // so, find out from which index to take this one value, and
+              // to which index to put
+              unsigned int out_index = 0;
+              for (unsigned int i=0; i<system_index; ++i)
+                out_index += this->n_nonzero_components(i);
+              unsigned int in_index = 0;
+              for (unsigned int i=0; i<base_index; ++i)
+                in_index += base_fe.n_nonzero_components(i);
+
+              // then loop over the number of components to be copied
+              Assert (this->n_nonzero_components(system_index) ==
+                      base_fe.n_nonzero_components(base_index),
+                      ExcInternalError());
+              for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
+                {
+                  if (base_flags & update_values)
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      output_data.shape_values[out_index+s][q] =
+                        base_data.shape_values(in_index+s,q);
+
+                  if (base_flags & update_gradients)
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      output_data.shape_gradients[out_index+s][q] =
+                        base_data.shape_gradients[in_index+s][q];
+
+                  if (base_flags & update_hessians)
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      output_data.shape_hessians[out_index+s][q] =
+                        base_data.shape_hessians[in_index+s][q];
+                }
+            }
+      }
+
+  Assert (base_no_mult_local_enriched_dofs.size() == fe_data.enrichment.size(),
+          ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),
+                               fe_data.enrichment.size()));
+  // calculate hessians, gradients and values for each function
+  for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
+    {
+      Assert (base_no_mult_local_enriched_dofs[base_no].size() == fe_data.enrichment[base_no].size(),
+              ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no].size(),
+                                   fe_data.enrichment[base_no].size()));
+      for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)
+        {
+          Assert (enrichments[base_no-1][m](cell) !=NULL,
+                  ExcMessage("The pointer to the enrichment function is NULL"));
+
+          Assert (enrichments[base_no-1][m](cell)->n_components == 1,
+                  ExcMessage("Only scalar-valued enrichment functions are allowed"));
+
+          if (flags & update_hessians)
+            {
+              Assert (fe_data.enrichment[base_no][m].hessians.size() == n_q_points,
+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].hessians.size(),
+                                           n_q_points));
+              for (unsigned int q=0; q<n_q_points; q++)
+                fe_data.enrichment[base_no][m].hessians[q]  = enrichments[base_no-1][m](cell)->hessian (mapping_data.quadrature_points[q]);
+            }
+
+          if (flags & update_gradients)
+            {
+              Assert (fe_data.enrichment[base_no][m].gradients.size() == n_q_points,
+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].gradients.size(),
+                                           n_q_points));
+              for (unsigned int q=0; q<n_q_points; q++)
+                fe_data.enrichment[base_no][m].gradients[q] = enrichments[base_no-1][m](cell)->gradient(mapping_data.quadrature_points[q]);
+            }
+
+          if (flags & update_values)
+            {
+              Assert (fe_data.enrichment[base_no][m].values.size() == n_q_points,
+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].values.size(),
+                                           n_q_points));
+              for (unsigned int q=0; q<n_q_points; q++)
+                fe_data.enrichment[base_no][m].values[q]    = enrichments[base_no-1][m](cell)->value   (mapping_data.quadrature_points[q]);
+            }
+        }
+    }
+
+  // Finally, update the standard data stored in output_data
+  // by expanding the product rule for enrichment function.
+  // note that the order if important, namely
+  // output_data.shape_XYZ contains values of standard FEM and we overwrite
+  // it with the updated one in the following order: hessians -> gradients -> values
+  if (flags & update_hessians)
+    {
+      for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
+        {
+          for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)
+            for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++)
+              {
+                const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i];
+                for (unsigned int q=0; q<n_q_points; ++q)
+                  {
+                    const Tensor<2, spacedim> grad_grad =  outer_product(output_data.shape_gradients[enriched_dof][q],
+                                                                         fe_data.enrichment[base_no][m].gradients[q]);
+                    const Tensor<2,spacedim,double> sym_grad_grad = symmetrize (grad_grad) * 2.0; // symmetrize does [s+s^T]/2
+
+                    output_data.shape_hessians   [enriched_dof][q]*= fe_data.enrichment[base_no][m].values[q];
+                    output_data.shape_hessians   [enriched_dof][q]+=
+                      sym_grad_grad +
+                      output_data.shape_values   [enriched_dof][q] * fe_data.enrichment[base_no][m].hessians[q];
+                  }
+              }
+        }
+    }
+
+  if (flags & update_gradients)
+    for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
+      {
+        for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)
+          for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++)
+            {
+              const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i];
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  output_data.shape_gradients[enriched_dof][q]*=  fe_data.enrichment[base_no][m].values[q];
+                  output_data.shape_gradients[enriched_dof][q]+=
+                    output_data.shape_values [enriched_dof][q] *  fe_data.enrichment[base_no][m].gradients[q];
+                }
+            }
+      }
+
+  if (flags & update_values)
+    for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
+      {
+        for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)
+          for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++)
+            {
+              const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i];
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  output_data.shape_values[enriched_dof][q] *= fe_data.enrichment[base_no][m].values[q];
+                }
+            }
+      }
+}
+
+
+template <int dim, int spacedim>
+const FESystem<dim,spacedim> &
+FE_Enriched<dim,spacedim>::
+get_fe_system() const
+{
+  return fe_system;
+}
+
+
+template <int dim, int spacedim>
+bool
+FE_Enriched<dim,spacedim>::
+hp_constraints_are_implemented () const
+{
+  return true;
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Enriched<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                               FullMatrix<double>                &matrix) const
+{
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&source))
+    {
+      fe_system.get_face_interpolation_matrix(fe_enr_other->get_fe_system(),
+                                              matrix);
+    }
+  else
+    {
+      typedef FiniteElement<dim,spacedim> FEL;
+      AssertThrow(false,typename FEL::ExcInterpolationNotImplemented());
+    }
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Enriched<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                                  const unsigned int        subface,
+                                  FullMatrix<double>       &matrix) const
+{
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&source))
+    {
+      fe_system.get_subface_interpolation_matrix(fe_enr_other->get_fe_system(),
+                                                 subface,
+                                                 matrix);
+    }
+  else
+    {
+      typedef FiniteElement<dim,spacedim> FEL;
+      AssertThrow(false,typename FEL::ExcInterpolationNotImplemented());
+    }
+}
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Enriched<dim,spacedim>::
+hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&fe_other))
+    {
+      return fe_system.hp_vertex_dof_identities(fe_enr_other->get_fe_system());
+    }
+  else
+    {
+      Assert (false, ExcNotImplemented());
+      return std::vector<std::pair<unsigned int, unsigned int> >();
+    }
+}
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Enriched<dim,spacedim>::
+hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&fe_other))
+    {
+      return fe_system.hp_line_dof_identities(fe_enr_other->get_fe_system());
+    }
+  else
+    {
+      Assert (false, ExcNotImplemented());
+      return std::vector<std::pair<unsigned int, unsigned int> >();
+    }
+}
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Enriched<dim,spacedim>::
+hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&fe_other))
+    {
+      return fe_system.hp_quad_dof_identities(fe_enr_other->get_fe_system());
+    }
+  else
+    {
+      Assert (false, ExcNotImplemented());
+      return std::vector<std::pair<unsigned int, unsigned int> >();
+    }
+}
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_Enriched<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+  // need to decide which element constrain another.
+  // for example Q(2) dominate Q(4) and thus some DoFs of Q(4) will be constrained.
+  // If we have Q(2) and Q(4)+POU, then it's clear that Q(2) dominates,
+  // namely our DoFs will be constrained to make field continuous.
+  // However, we need to check for situations like Q(4) vs Q(2)+POU.
+  // In that case the domination for the underlying FEs should be the other way,
+  // but this implies that we can't constrain POU dofs to make the field continuous.
+  // In that case, through an error
+
+  // if it's also enriched, do domination based on each one's FESystem
+  if (const FE_Enriched<dim,spacedim> *fe_enr_other
+      = dynamic_cast<const FE_Enriched<dim,spacedim>*>(&fe_other))
+    {
+      return fe_system.compare_for_face_domination(fe_enr_other->get_fe_system());
+    }
+  else
+    {
+      Assert (false, ExcNotImplemented());
+      return FiniteElementDomination::neither_element_dominates;
+    }
+}
+
+/* ----------------------- FESystem::InternalData ------------------- */
+
+
+template <int dim, int spacedim>
+FE_Enriched<dim,spacedim>::InternalData::InternalData(std::unique_ptr<typename FESystem<dim,spacedim>::InternalData> fesystem_data)
+  :
+  fesystem_data(std::move(fesystem_data))
+{}
+
+
+template <int dim, int spacedim>
+typename FiniteElement<dim,spacedim>::InternalDataBase &
+FE_Enriched<dim,spacedim>::
+InternalData::get_fe_data (const unsigned int base_no) const
+{
+  return fesystem_data->get_fe_data(base_no);
+}
+
+
+template <int dim, int spacedim>
+internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
+FE_Enriched<dim,spacedim>::
+InternalData::get_fe_output_object (const unsigned int base_no) const
+{
+  return fesystem_data->get_fe_output_object(base_no);
+}
+
+
+// explicit instantiations
+#include "fe_enriched.inst"
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // CXX11
diff --git a/source/fe/fe_enriched.inst.in b/source/fe/fe_enriched.inst.in
new file mode 100644 (file)
index 0000000..853f36a
--- /dev/null
@@ -0,0 +1,23 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    template class FE_Enriched<deal_II_dimension, deal_II_space_dimension>;
+#endif
+}
+
diff --git a/tests/fe/fe_enriched_01.cc b/tests/fe/fe_enriched_01.cc
new file mode 100644 (file)
index 0000000..9413d6f
--- /dev/null
@@ -0,0 +1,127 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 basic output of the class (name, n_blocks, n_compoennts, n_dofs_per_cell)
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+using namespace dealii;
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+template<int dim>
+void test_base()
+{
+  deallog << "Test basic functions:"<<std::endl;
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+  deallog << fe.get_name()<<std::endl;
+  deallog << fe.n_blocks()<<std::endl;
+  deallog << fe.n_components() << std::endl;
+  deallog << fe.n_dofs_per_cell() << std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test_base<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;
+    };
+}
diff --git a/tests/fe/fe_enriched_01.with_cxx11=on.output b/tests/fe/fe_enriched_01.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..d5e7fca
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::Test basic functions:
+DEAL::FE_Enriched<3>[FE_Q<3>(1)-FE_Q<3>(1)]
+DEAL::2
+DEAL::1
+DEAL::16
diff --git a/tests/fe/fe_enriched_02.cc b/tests/fe/fe_enriched_02.cc
new file mode 100644 (file)
index 0000000..fb6ac20
--- /dev/null
@@ -0,0 +1,163 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check FEValues::shape_value()
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+template<int dim>
+void test1()
+{
+  deallog << "FEValues.shape_value()"<<std::endl;
+  deallog << "for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)"<<std::endl;
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<dim> quadrature(1);
+  FEValues<dim> fe_values(fe, quadrature,
+                          update_values | update_gradients |
+                          update_JxW_values);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      const unsigned int n_q_points = quadrature.size();
+      const unsigned int dofs_per_cell = fe.dofs_per_cell;
+      const std::vector<dealii::Point<dim> > &q_points = fe_values.get_quadrature_points();
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          deallog   <<"dof="<<i<<" qp="<<q_points[q_point]
+                    <<" f(qp)="<<function.value(q_points[q_point])
+                    <<" N(qp)="<<fe_values.shape_value(i,q_point)<<std::endl;
+    }
+  dof_handler.clear();
+  deallog << std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test1<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;
+    };
+}
diff --git a/tests/fe/fe_enriched_02.with_cxx11=on.output b/tests/fe/fe_enriched_02.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..00053ee
--- /dev/null
@@ -0,0 +1,20 @@
+
+DEAL::FEValues.shape_value()
+DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)
+DEAL::dof=0 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=1 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=2 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=3 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=4 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=5 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=6 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=7 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=8 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=9 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=10 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=11 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=12 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=13 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::dof=14 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250
+DEAL::dof=15 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526
+DEAL::
diff --git a/tests/fe/fe_enriched_03.cc b/tests/fe/fe_enriched_03.cc
new file mode 100644 (file)
index 0000000..73ec999
--- /dev/null
@@ -0,0 +1,157 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check FEFaceValues::shape_value()
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+using namespace dealii;
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template <int dim>
+void test2()
+{
+  deallog << "FEFaceValues"<<std::endl;
+  deallog << "for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)"<<std::endl;
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<dim-1> quadrature(1);
+  FEFaceValues<dim> fe_face_values(fe, quadrature,
+                                   update_values | update_gradients |
+                                   update_JxW_values);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+      {
+        fe_face_values.reinit (cell,face);
+        const unsigned int n_q_points = quadrature.size();
+        const unsigned int dofs_per_cell = fe.dofs_per_cell;
+        const std::vector<dealii::Point<dim> > &q_points = fe_face_values.get_quadrature_points();
+
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+            deallog   <<"dof="<<i<<" qp="<<q_points[q_point]
+                      <<" f(qp)="<<function.value(q_points[q_point])
+                      <<" N(qp)="<<fe_face_values.shape_value(i,q_point)<<std::endl;
+      }
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test2<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;
+    };
+}
diff --git a/tests/fe/fe_enriched_03.with_cxx11=on.output b/tests/fe/fe_enriched_03.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..ef9efb9
--- /dev/null
@@ -0,0 +1,99 @@
+
+DEAL::FEFaceValues
+DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)
+DEAL::dof=0 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=1 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=2 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=3 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=4 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=5 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=6 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=7 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=8 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=9 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=10 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=11 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=12 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=13 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=14 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=15 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=0 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=1 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=2 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=3 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=4 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=5 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=6 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=7 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=8 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=9 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=10 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=11 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=12 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=13 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=14 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=15 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=0 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=1 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=2 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=3 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=4 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=5 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=6 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=7 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=8 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=9 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=10 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=11 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=12 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=13 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=14 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=15 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=0 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=1 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=2 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=3 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=4 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=5 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=6 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=7 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=8 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=9 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=10 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=11 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=12 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=13 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=14 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=15 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=0 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=1 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=2 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=3 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=4 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=5 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=6 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500
+DEAL::dof=7 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233
+DEAL::dof=8 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=9 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=10 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=11 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=12 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=13 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=14 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=15 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0
+DEAL::dof=0 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=1 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=2 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=3 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=4 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=5 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=6 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=7 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0
+DEAL::dof=8 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=9 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=10 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=11 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=12 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=13 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735
+DEAL::dof=14 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500
+DEAL::dof=15 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735
diff --git a/tests/fe/fe_enriched_04.cc b/tests/fe/fe_enriched_04.cc
new file mode 100644 (file)
index 0000000..2d3d921
--- /dev/null
@@ -0,0 +1,244 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 FEValues::get_function_values()
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template<int dim>
+void test3()
+{
+  deallog << "FEValues.get_function_values()"<<std::endl;
+  deallog << "for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)"<<std::endl;
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values(fe, quadrature,
+                          update_values | update_gradients |
+                          update_JxW_values);
+
+  Vector<double> solution_fe(dof_handler.n_dofs()),
+         solution_pou(dof_handler.n_dofs()),
+         solution(dof_handler.n_dofs());
+  solution_fe = 0;
+  solution_pou = 0;
+
+  // groups are one after another
+  // thus set to 1.0 the same shape function in the underlying FE
+  solution_fe[0] = 2.0;
+  solution_pou[1] = 2.0;
+  solution = solution_fe;
+  solution+= solution_pou;
+
+  const unsigned int n_q_points = quadrature.size();
+
+  std::vector<double> solution_values_fe(n_q_points),
+      solution_values_pou(n_q_points),
+      solution_values(n_q_points);
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+
+      const unsigned int dofs_per_cell = fe.dofs_per_cell;
+      const std::vector<dealii::Point<dim> > &q_points = fe_values.get_quadrature_points();
+      fe_values.get_function_values (solution_fe,
+                                     solution_values_fe);
+      fe_values.get_function_values (solution_pou,
+                                     solution_values_pou);
+      fe_values.get_function_values (solution,
+                                     solution_values);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        deallog   <<" qp="<<q_points[q_point]
+                  <<" f(qp)="<<function.value(q_points[q_point])
+                  <<" U_fe(qp)="<<solution_values_fe[q_point]
+                  <<" U_pou(qp)="<<solution_values_pou[q_point]
+                  <<" U_[fe+pou](qp)="<<solution_values[q_point]
+                  <<std::endl;
+    }
+  dof_handler.clear();
+}
+
+template <int dim>
+void plot_shape_function()
+{
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<Vector<double>> shape_functions(dof_handler.n_dofs());
+  std::vector<std::string> names;
+  for (unsigned int s=0; s < shape_functions.size(); s++)
+    {
+      names.push_back(std::string("N_") +
+                      dealii::Utilities::int_to_string(s));
+
+      shape_functions[s].reinit(dof_handler.n_dofs());
+
+      // this is ok for 1 cell only!
+      const unsigned int group = fe.system_to_base_index(s).first.first;
+
+      // zero everywhere but
+      // the DoF corresponding to this shape function
+      if (group==0)
+        shape_functions[s][s] = 1.0;
+      else
+        shape_functions[s][s] = 1.0;//can potentially put another value here
+    }
+
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+
+  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_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test3<3>();
+#ifdef DATA_OUT_FE_ENRICHED
+      plot_shape_function<1>();
+      plot_shape_function<2>();
+      plot_shape_function<3>();
+#endif
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_04.with_cxx11=on.output b/tests/fe/fe_enriched_04.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..785ca6a
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::FEValues.get_function_values()
+DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp)
+DEAL:: qp=0.2113 0.2113 0.2113 f(qp)=0.6935 U_fe(qp)=0.9811 U_pou(qp)=0.6804 U_[fe+pou](qp)=1.6615
+DEAL:: qp=0.7887 0.2113 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760
+DEAL:: qp=0.2113 0.7887 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760
+DEAL:: qp=0.7887 0.7887 0.2113 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931
+DEAL:: qp=0.2113 0.2113 0.7887 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760
+DEAL:: qp=0.7887 0.2113 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931
+DEAL:: qp=0.2113 0.7887 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931
+DEAL:: qp=0.7887 0.7887 0.7887 f(qp)=0.2551 U_fe(qp)=0.0189 U_pou(qp)=0.0048 U_[fe+pou](qp)=0.0237
diff --git a/tests/fe/fe_enriched_05.cc b/tests/fe/fe_enriched_05.cc
new file mode 100644 (file)
index 0000000..6a2a63c
--- /dev/null
@@ -0,0 +1,186 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 h-refinement with DoFHandler and print constraints.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template <int dim>
+void test4()
+{
+  deallog << "h-refinement:"<<std::endl;
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(2),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+
+  // one global refinement and 1 refinement of 1 active cell (local)
+  triangulation.refine_global();
+  triangulation.begin_active()->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+
+  dof_handler.distribute_dofs(fe);
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints.close ();
+
+  constraints.print(deallog.get_file_stream());
+
+#ifdef DATA_OUT_FE_ENRICHED
+  std::vector<Vector<double>> shape_functions;
+  std::vector<std::string> names;
+  for (unsigned int s=0; s < dof_handler.n_dofs(); s++)
+    if (!constraints.is_constrained(s))
+      {
+        names.push_back(std::string("N_") +
+                        dealii::Utilities::int_to_string(s));
+
+        Vector<double> shape_function;
+        shape_function.reinit(dof_handler.n_dofs());
+        shape_function[s] = 1.0;
+
+        // make continuous:
+        constraints.distribute(shape_function);
+        shape_functions.push_back(shape_function);
+      }
+
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+
+  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 = "h-shape_functions_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+#endif
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test4<2>();
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_05.with_cxx11=on.output b/tests/fe/fe_enriched_05.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..764434a
--- /dev/null
@@ -0,0 +1,20 @@
+
+DEAL::h-refinement:
+    42 8:  1.00000
+    43 1:  0.500000
+    43 5:  0.500000
+    44 0:  0.375000
+    44 4:  -0.125000
+    44 8:  0.750000
+    48 21:  1.00000
+    49 5:  0.500000
+    49 14:  0.500000
+    52 4:  -0.125000
+    52 13:  0.375000
+    52 21:  0.750000
+    54 0:  -0.125000
+    54 4:  0.375000
+    54 8:  0.750000
+    55 4:  0.375000
+    55 13:  -0.125000
+    55 21:  0.750000
diff --git a/tests/fe/fe_enriched_06.cc b/tests/fe/fe_enriched_06.cc
new file mode 100644 (file)
index 0000000..01b33e1
--- /dev/null
@@ -0,0 +1,223 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// no refinement, just 4 cells.
+// Make sure fe_values and shape functions are continuous
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template <int dim>
+void test5()
+{
+  deallog << "4 cells:"<<std::endl;
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe(FE_Q<dim>(1),
+                      FE_Q<dim>(1),
+                      &function);
+
+  GridGenerator::hyper_cube(triangulation);
+
+  triangulation.refine_global();
+
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<Vector<double>> shape_functions;
+  std::vector<std::string> names;
+  for (unsigned int s=0; s < dof_handler.n_dofs(); s++)
+    {
+      names.push_back(std::string("N_") +
+                      dealii::Utilities::int_to_string(s));
+
+      Vector<double> shape_function;
+      shape_function.reinit(dof_handler.n_dofs());
+      shape_function[s] = 1.0;
+
+      shape_functions.push_back(shape_function);
+    }
+
+  // output 11th:
+  {
+    const unsigned int global_dof = 11;
+    const Vector<double> &solution = shape_functions[global_dof];
+    QTrapez<dim> quadrature;
+    FEValues<dim> fe_values(fe, quadrature,
+                            update_values);
+
+    const unsigned int n_q_points = quadrature.size();
+    std::vector<double> solution_values(n_q_points);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    std::vector<dealii::types::global_dof_index> local_dof_indices;
+    local_dof_indices.resize (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active (),
+    endc = dof_handler.end ();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+
+        // find out which
+        unsigned int local_dof = 0;
+        cell->get_dof_indices (local_dof_indices);
+        for ( ; local_dof < dofs_per_cell; local_dof++)
+          if (local_dof_indices[local_dof]==global_dof)
+            break;
+
+        const std::vector<dealii::Point<dim> > &q_points = fe_values.get_quadrature_points();
+        fe_values.get_function_values (solution,
+                                       solution_values);
+
+        deallog   <<" cell="<<cell->center()<<std::endl;
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            // find non-zero shape_value
+            deallog   <<" qp="<<q_points[q_point]
+                      <<" f(qp)="<<function.value(q_points[q_point]);
+
+            // if the cell contains our global dof:
+            if (local_dof < dofs_per_cell)
+              deallog <<" N("<<local_dof<<",qp)="<<fe_values.shape_value(local_dof,q_point);
+
+            deallog <<" U(qp)="<<solution_values[q_point]
+                    <<std::endl;
+          }
+      }
+  }
+
+#ifdef DATA_OUT_FE_ENRICHED
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+
+  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 = "4cell_functions_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+#endif
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test5<2>();
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_06.with_cxx11=on.output b/tests/fe/fe_enriched_06.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..a45dcea
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL::4 cells:
+DEAL:: cell=0.2500 0.2500
+DEAL:: qp=0.0000 0.0000 f(qp)=1.0000 U(qp)=0
+DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 U(qp)=0
+DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0
+DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0
+DEAL:: cell=0.7500 0.2500
+DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 N(7,qp)=0 U(qp)=0
+DEAL:: qp=1.0000 0.0000 f(qp)=0.3679 N(7,qp)=0 U(qp)=0
+DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(7,qp)=0 U(qp)=0
+DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(7,qp)=0.3269 U(qp)=0.3269
+DEAL:: cell=0.2500 0.7500
+DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0
+DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0
+DEAL:: qp=0.0000 1.0000 f(qp)=0.3679 U(qp)=0
+DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 U(qp)=0
+DEAL:: cell=0.7500 0.7500
+DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(3,qp)=0 U(qp)=0
+DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(3,qp)=0.3269 U(qp)=0.3269
+DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 N(3,qp)=0 U(qp)=0
+DEAL:: qp=1.0000 1.0000 f(qp)=0.2431 N(3,qp)=0 U(qp)=0
diff --git a/tests/fe/fe_enriched_07.cc b/tests/fe/fe_enriched_07.cc
new file mode 100644 (file)
index 0000000..d051202
--- /dev/null
@@ -0,0 +1,237 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// hp-refinement. Output constraints.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template <int dim>
+void test6(const bool do_href,
+           const unsigned int p_feq  = 1,
+           const unsigned int p_feen = 2)
+{
+  deallog << "hp: "<<do_href<<" "<<p_feq<<" "<<p_feen<<std::endl;
+
+  Triangulation<dim> triangulation;
+  hp::DoFHandler<dim> dof_handler(triangulation);
+
+  EnrichmentFunction<dim> function;
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feq)));
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feen),
+                                           FE_Q<dim>(1),
+                                           &function));
+  // push back to be able to resolve hp constrains:
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feen)));
+
+  GridGenerator::hyper_cube(triangulation);
+
+  // one global refinement and 1 refinement of 1 active cell
+  // a sketch with material ids:
+  // ---------------
+  // |      |      |
+  // |  1   |  0   |
+  // |------|-------
+  // |1  |1 |      |
+  // |------|  0   |
+  // |1  |1 |      |
+  // |------|------|
+  triangulation.refine_global();
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active();
+    cell->set_active_fe_index(1);//POU
+    cell++;
+    cell->set_active_fe_index(1);//POU
+  }
+  dof_handler.distribute_dofs(fe_collection);
+
+  if (do_href)
+    {
+      triangulation.begin_active()->set_refine_flag ();
+      triangulation.execute_coarsening_and_refinement ();
+      dof_handler.distribute_dofs(fe_collection);
+    }
+
+  deallog << "n_cells: "<<triangulation.n_active_cells()<<std::endl;
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints.close ();
+
+  constraints.print(deallog.get_file_stream());
+
+#ifdef DATA_OUT_FE_ENRICHED
+  // 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(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_ref="+std::to_string(do_href)
+                         +"_p_feq="+std::to_string(p_feq)
+                         +"_p_feenr="+std::to_string(p_feen)
+                         +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+#endif
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test6<2>(false,1,2);// 1 vs 2+1
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_07.with_cxx11=on.output b/tests/fe/fe_enriched_07.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..e4f68ec
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::hp: 0 1 2
+DEAL::n_cells: 4
+    4 = 0
+    5 = 0
+    9 18:  0.500000
+    9 19:  0.500000
+    13 = 0
+    16 19:  0.500000
+    16 22:  0.500000
diff --git a/tests/fe/fe_enriched_08.cc b/tests/fe/fe_enriched_08.cc
new file mode 100644 (file)
index 0000000..c59363d
--- /dev/null
@@ -0,0 +1,219 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check constraints on hp dofhandler with 2 neighboring cells where one
+// cell is enriched and nother is not.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// argument for build_patches()
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+// uncomment when debugging
+// #define DATA_OUT_FE_ENRICHED
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+};
+
+
+template<int dim>
+void test2cells(const unsigned int p_feq=2,
+                const unsigned int p_feen=1)
+{
+  deallog << "2cells: "<<dim<<" "<<p_feq<<" "<<p_feen<<std::endl;
+  Triangulation<dim>   triangulation;
+  {
+    Triangulation<dim>   triangulationL;
+    Triangulation<dim>   triangulationR;
+    GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain
+    GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain
+    Point<dim> shift_vector;
+    shift_vector[0] = 1.0;
+    GridTools::shift(shift_vector,triangulationR);
+    GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation);
+  }
+
+  hp::DoFHandler<dim> dof_handler(triangulation);
+  EnrichmentFunction<dim> function;
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feq)));
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feen),
+                                           FE_Q<dim>(1),
+                                           &function));
+  // push back to be able to resolve hp constrains:
+  fe_collection.push_back(FE_Enriched<dim>(FE_Q<dim>(p_feen)));
+
+  dof_handler.begin_active()->set_active_fe_index(1);//POU
+
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints.close ();
+
+  constraints.print(deallog.get_file_stream());
+
+#ifdef DATA_OUT_FE_ENRICHED
+  // 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(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 = "2cell-shape_functions_p_feq="+std::to_string(p_feq)
+                         +"_p_feenr="+std::to_string(p_feen)
+                         +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+#endif
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test2cells<2>(1,2); // 1 vs 2+1
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_08.with_cxx11=on.output b/tests/fe/fe_enriched_08.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..3a075cc
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::2cells: 2 1 2
+    2 = 0
+    5 = 0
+    7 11:  0.500000
+    7 13:  0.500000
diff --git a/tests/fe/fe_enriched_08a.cc b/tests/fe/fe_enriched_08a.cc
new file mode 100644 (file)
index 0000000..23fd6db
--- /dev/null
@@ -0,0 +1,133 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Same as fe_enriched_08, but checks the constraints from FESystem directly.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+using namespace dealii;
+
+template<int dim>
+void test2cellsFESystem(const unsigned int p_feq=2,
+                        const unsigned int p_feen=1)
+{
+  deallog << "2cells: "<<dim<<" "<<p_feq<<" "<<p_feen<<std::endl;
+  Triangulation<dim>   triangulation;
+  {
+    Triangulation<dim>   triangulationL;
+    Triangulation<dim>   triangulationR;
+    GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain
+    GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain
+    Point<dim> shift_vector;
+    shift_vector[0] = 1.0;
+    GridTools::shift(shift_vector,triangulationR);
+    GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation);
+  }
+
+  hp::DoFHandler<dim> dof_handler(triangulation);
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p_feq),1,
+                                        FE_Nothing<dim>(),1));
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p_feen),1,
+                                        FE_Q<dim>(1),1));
+
+  // push back to be able to resolve hp constrains:
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p_feen),1,
+                                        FE_Nothing<dim>(),1));
+
+  dof_handler.begin_active()->set_active_fe_index(1);
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints.close ();
+
+  constraints.print(deallog.get_file_stream());
+
+  dof_handler.clear();
+}
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test2cellsFESystem<2>(1,2); // 1+0 vs 2+1
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_08a.with_cxx11=on.output b/tests/fe/fe_enriched_08a.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..3a075cc
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::2cells: 2 1 2
+    2 = 0
+    5 = 0
+    7 11:  0.500000
+    7 13:  0.500000
diff --git a/tests/fe/fe_enriched_compare_to_fe_system.cc b/tests/fe/fe_enriched_compare_to_fe_system.cc
new file mode 100644 (file)
index 0000000..f19cb7c
--- /dev/null
@@ -0,0 +1,320 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 FE_Enriched class by comparing with FE_System for a single enrichment
+// function.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// create a mesh without hanging nodes,
+// move nodes randomly
+// and compare FEEnriched to FESystem with explicit treatment of
+// the product rule for
+// valus, gradients, hessians on
+// elements and faces.
+// The comparison is straight forward because local dofs are enumerated
+// in the same way for FE_System and FEEnriched.
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    return std::exp(-point.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> res = point;
+    Assert (point.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    res *= -value(point)/point.norm();
+    return res;
+  }
+
+  virtual SymmetricTensor<2,dim> hessian (const Point<dim> &p,
+                                          const unsigned int component=0) const
+  {
+    Tensor<1,dim> dir = p;
+    const double r = dir.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    dir/=r;
+    SymmetricTensor<2,dim> dir_x_dir;
+    for (unsigned int i=0; i<dim; i++)
+      for (unsigned int j=i; j<dim; j++)
+        dir_x_dir[i][j] = dir[i] * dir[j];
+
+    return std::exp(-r)*( (1.0+0.5/r)*dir_x_dir - unit_symmetric_tensor<dim>()/r );
+  }
+
+};
+
+/**
+ *
+ * @param p     point
+ * @param func  function
+ * @param v_e   value of enriched FE
+ * @param g_e   gradient of enriched FE
+ * @param h_e   hessian of enriched FE
+ * @param v_s0  value of 0th component of FE_System
+ * @param g_s0  gradient of 0th component of FE_System
+ * @param h_s0  hessian of 0th component of FE_System
+ * @param v_s1  value of 1st component of FE_System
+ * @param g_s1  gradient of 1st component of FE_System
+ * @param h_s1  hessian of 1st component of FE_System
+ */
+template <int dim>
+void check_consistency(const Point<dim>    &p,
+                       const Function<dim> &func,
+                       const double        &v_e,
+                       const Tensor<1,dim> &g_e,
+                       const Tensor<2,dim> &h_e,
+                       const double        &v_s0,
+                       const Tensor<1,dim> &g_s0,
+                       const Tensor<2,dim> &h_s0,
+                       const double        &v_s1,
+                       const Tensor<1,dim> &g_s1,
+                       const Tensor<2,dim> &h_s1)
+{
+  const double                 v_f = func.value(p);
+  const Tensor<1,dim>          g_f = func.gradient(p);
+  const SymmetricTensor<2,dim> h_f = func.hessian(p);
+
+  // product rule:
+  const double                 v_s = v_s0 + v_s1 * v_f;
+  const Tensor<1,dim>          g_s = g_s0 + g_s1 * v_f + v_s1 * g_f;
+  Tensor<2,dim>                 op = outer_product(g_s1 , g_f );
+
+  const SymmetricTensor<2,dim> sp = symmetrize (op) * 2.0; // symmetrize does [s+s^T]/2
+  // Hessians should be symmetric, but due to round off errors and the fact
+  // that deal.II stores hessians in tensors, we may end up with failing check
+  // for symmetry in SymmetricTensor class.
+  // For a moment stick with usual Tensor
+  const Tensor<2,dim> h_s = h_s0 + h_s1 * v_f + sp + v_s1 * h_f;
+
+  const double        v_d = v_s - v_e;
+  AssertThrow ( std::abs(v_d) < eps,
+                ExcInternalError());
+  const Tensor<1,dim> g_d = g_s - g_e;
+  AssertThrow ( g_d.norm() < eps,
+                ExcInternalError());
+
+  // see note above.
+  const Tensor<2,dim> h_d = h_s - h_e;
+  AssertThrow ( h_d.norm() < eps,
+                ExcInternalError());
+}
+
+template<int dim>
+void test(const FiniteElement<dim> &fe1,
+          const FiniteElement<dim> &fe2,
+          const Quadrature<dim>    &volume_quad,
+          const Quadrature<dim-1>  &face_quad,
+          const bool                distort)
+{
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler_enriched(triangulation);
+  DoFHandler<dim> dof_handler_system(triangulation);
+
+  EnrichmentFunction<dim> function;
+  FE_Enriched<dim> fe_enriched(fe1,
+                               fe2,
+                               &function);
+  FESystem<dim> fe_system(fe1,1,
+                          fe2,1);
+
+  {
+    Point<dim> p1,p2;
+    std::vector<unsigned int> repetitions(dim);
+    for (unsigned int d=0; d< dim; d++)
+      {
+        p1[d]=-1.0;
+        p2[d]=2.0;
+        repetitions[d] = 3;
+      }
+    GridGenerator::subdivided_hyper_rectangle(triangulation,
+                                              repetitions,
+                                              p1,p2);
+
+    if (distort)
+      GridTools::distort_random (0.1,triangulation);
+  }
+
+  dof_handler_enriched.distribute_dofs(fe_enriched);
+  dof_handler_system.distribute_dofs(fe_system);
+
+  FEValues<dim> fe_values_enriched(fe_enriched, volume_quad,
+                                   update_values | update_gradients | update_hessians );
+
+  FEValues<dim> fe_values_system(fe_system, volume_quad,
+                                 update_values | update_gradients | update_hessians | update_quadrature_points);
+
+  FEFaceValues<dim> fe_face_values_enriched(fe_enriched, face_quad,
+                                            update_values | update_gradients | update_hessians);
+
+  FEFaceValues<dim> fe_face_values_system  (fe_system, face_quad,
+                                            update_values | update_gradients | update_hessians | update_quadrature_points);
+
+  const unsigned int dofs_per_cell = fe_enriched.dofs_per_cell;
+  Assert (fe_enriched.dofs_per_cell == fe_system.dofs_per_cell,
+          ExcInternalError());
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell_enriched = dof_handler_enriched.begin_active (),
+  endc_enriched = dof_handler_enriched.end (),
+  cell_system   = dof_handler_system.begin_active(),
+  endc_system   = dof_handler_system.end();
+  for (; cell_enriched!=endc_enriched; ++cell_enriched, ++cell_system)
+    {
+      fe_values_enriched.reinit (cell_enriched);
+      fe_values_system.reinit   (cell_system);
+      const unsigned int n_q_points = volume_quad.size();
+      const std::vector<dealii::Point<dim> > &q_points = fe_values_system.get_quadrature_points();
+
+      // check shape functions
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          check_consistency(q_points[q_point], function,
+                            fe_values_enriched.shape_value  (i,q_point),
+                            fe_values_enriched.shape_grad   (i,q_point),
+                            fe_values_enriched.shape_hessian(i,q_point),
+                            fe_values_system.shape_value_component   (i,q_point,0),
+                            fe_values_system.shape_grad_component    (i,q_point,0),
+                            fe_values_system.shape_hessian_component (i,q_point,0),
+                            fe_values_system.shape_value_component   (i,q_point,1),
+                            fe_values_system.shape_grad_component    (i,q_point,1),
+                            fe_values_system.shape_hessian_component (i,q_point,1));
+
+      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+        {
+          fe_face_values_enriched.reinit (cell_enriched,face);
+          fe_face_values_system.reinit (cell_system,face);
+          const unsigned int n_q_points_face = face_quad.size();
+          const std::vector<dealii::Point<dim> > &q_points = fe_face_values_system.get_quadrature_points();
+
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int q_point=0; q_point<n_q_points_face; ++q_point)
+              check_consistency(q_points[q_point], function,
+                                fe_face_values_enriched.shape_value  (i,q_point),
+                                fe_face_values_enriched.shape_grad   (i,q_point),
+                                fe_face_values_enriched.shape_hessian(i,q_point),
+                                fe_face_values_system.shape_value_component   (i,q_point,0),
+                                fe_face_values_system.shape_grad_component    (i,q_point,0),
+                                fe_face_values_system.shape_hessian_component (i,q_point,0),
+                                fe_face_values_system.shape_value_component   (i,q_point,1),
+                                fe_face_values_system.shape_grad_component    (i,q_point,1),
+                                fe_face_values_system.shape_hessian_component (i,q_point,1));
+        }
+    }
+
+
+  deallog << "Ok" << std::endl;
+  dof_handler_enriched.clear();
+  dof_handler_system.clear();
+}
+
+
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      {
+        const unsigned int dim = 2;
+        test(FE_Q<dim>(3),FE_Q<dim>(2),QGauss<dim>(2),QGauss<dim-1>(2),false);
+        test(FE_Q<dim>(3),FE_Q<dim>(2),QGauss<dim>(2),QGauss<dim-1>(2),true);
+      }
+
+      {
+        const unsigned int dim = 3;
+        test(FE_Q<dim>(3),FE_Q<dim>(2),QGauss<dim>(2),QGauss<dim-1>(2),false);
+        test(FE_Q<dim>(3),FE_Q<dim>(2),QGauss<dim>(2),QGauss<dim-1>(2),true);
+      }
+
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..f2b98e8
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.cc b/tests/fe/fe_enriched_compare_to_fe_system_2.cc
new file mode 100644 (file)
index 0000000..ca0e3c8
--- /dev/null
@@ -0,0 +1,403 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 FE_Enriched class by comparing with FE_System for multiple enrichment
+// functions.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+// create a mesh without hanging nodes,
+// move nodes randomly
+// and compare FEEnriched to FESystem with explicit treatment of
+// the product rule for
+// valus, gradients, hessians on
+// elements and faces.
+// The comparison is straight forward because local dofs are enumerated
+// in the same way for FE_System and FEEnriched.
+// Same as fe_enriched_compare_to_fe_system.cc but now multiple enrichment functions.
+const unsigned int patches = 10;
+
+using namespace dealii;
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction(const Point<dim> &origin)
+    : Function<dim>(1),
+      origin(origin)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0 ) const
+  {
+    Tensor<1,dim> d = point-origin;
+    return std::exp(-d.norm());
+  }
+
+  virtual Tensor<1,dim> gradient(const Point<dim> &point,
+                                 const unsigned int component = 0) const
+  {
+    Tensor<1,dim> d = point -origin;
+    Assert (d.norm() > 0,
+            dealii::ExcMessage("gradient is not defined at zero"));
+    d *= -std::exp(-d.norm())/d.norm();
+    return d;
+  }
+
+  virtual SymmetricTensor<2,dim> hessian (const Point<dim> &p,
+                                          const unsigned int component=0) const
+  {
+    Tensor<1,dim> dir = p - origin;
+    const double r = dir.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    dir/=r;
+    SymmetricTensor<2,dim> dir_x_dir;
+    for (unsigned int i=0; i<dim; i++)
+      for (unsigned int j=i; j<dim; j++)
+        dir_x_dir[i][j] = dir[i] * dir[j];
+
+    return std::exp(-r)*( (1.0+0.5/r)*dir_x_dir - unit_symmetric_tensor<dim>()/r );
+  }
+
+private:
+  Point<dim> origin;
+
+};
+
+/**
+ *
+ * @param p     point
+ * @param func  function
+ * @param v_e   value of enriched FE
+ * @param g_e   gradient of enriched FE
+ * @param h_e   hessian of enriched FE
+ * @param v_s0  value of 0th component of FE_System
+ * @param g_s0  gradient of 0th component of FE_System
+ * @param h_s0  hessian of 0th component of FE_System
+ * @param v_s1  value of 1st component of FE_System
+ * @param g_s1  gradient of 1st component of FE_System
+ * @param h_s1  hessian of 1st component of FE_System
+ */
+template <int dim>
+void check_consistency(const Point<dim>    &p,
+                       const Function<dim> &func1,
+                       const Function<dim> &func2,
+                       const Function<dim> &func3,
+                       const double        &v_e,
+                       const Tensor<1,dim> &g_e,
+                       const Tensor<2,dim> &h_e,
+                       const double        &v_s0,
+                       const Tensor<1,dim> &g_s0,
+                       const Tensor<2,dim> &h_s0,
+                       const double        &v_s1,
+                       const Tensor<1,dim> &g_s1,
+                       const Tensor<2,dim> &h_s1,
+                       const double        &v_s2,
+                       const Tensor<1,dim> &g_s2,
+                       const Tensor<2,dim> &h_s2,
+                       const double        &v_s3,
+                       const Tensor<1,dim> &g_s3,
+                       const Tensor<2,dim> &h_s3)
+{
+  const double                 v_f1 = func1.value(p);
+  const Tensor<1,dim>          g_f1 = func1.gradient(p);
+  const SymmetricTensor<2,dim> h_f1 = func1.hessian(p);
+
+  const double                 v_f2 = func2.value(p);
+  const Tensor<1,dim>          g_f2 = func2.gradient(p);
+  const SymmetricTensor<2,dim> h_f2 = func2.hessian(p);
+
+  const double                 v_f3 = func3.value(p);
+  const Tensor<1,dim>          g_f3 = func3.gradient(p);
+  const SymmetricTensor<2,dim> h_f3 = func3.hessian(p);
+
+
+
+  // product rule:
+  const double                 v_s = v_s0 + v_s1 * v_f1
+                                     + v_s2 * v_f2
+                                     + v_s3 * v_f3;
+  const Tensor<1,dim>          g_s = g_s0 + g_s1 * v_f1 + v_s1 * g_f1
+                                     + g_s2 * v_f2 + v_s2 * g_f2
+                                     + g_s3 * v_f3 + v_s3 * g_f3;
+  Tensor<2,dim> op1 = outer_product(g_s1 , g_f1 ),
+                op2 = outer_product(g_s2 , g_f2 ),
+                op3 = outer_product(g_s3 , g_f3 );
+
+  const SymmetricTensor<2,dim> sp1 = symmetrize (op1) * 2.0; // symmetrize does [s+s^T]/2
+  const SymmetricTensor<2,dim> sp2 = symmetrize (op2) * 2.0; // symmetrize does [s+s^T]/2
+  const SymmetricTensor<2,dim> sp3 = symmetrize (op3) * 2.0; // symmetrize does [s+s^T]/2
+  // Hessians should be symmetric, but due to round off errors and the fact
+  // that deal.II stores hessians in tensors, we may end up with failing check
+  // for symmetry in SymmetricTensor class.
+  // For a moment stick with usual Tensor
+  const Tensor<2,dim> h_s = h_s0 + h_s1 * v_f1 + sp1 + v_s1 * h_f1
+                            + h_s2 * v_f2 + sp2 + v_s2 * h_f2
+                            + h_s3 * v_f3 + sp3 + v_s3 * h_f3;
+
+  const double        v_d = v_s - v_e;
+  AssertThrow ( std::abs(v_d) < eps,
+                ExcInternalError());
+  const Tensor<1,dim> g_d = g_s - g_e;
+  AssertThrow ( g_d.norm() < eps,
+                ExcInternalError());
+
+  // see note above.
+  const Tensor<2,dim> h_d = h_s - h_e;
+  AssertThrow ( h_d.norm() < eps,
+                ExcInternalError());
+}
+
+/**
+ * {fe_en1,1}, {fe_en2,2} => 3 functions
+ */
+template<int dim>
+void test(const FiniteElement<dim> &fe_base,
+          const FiniteElement<dim> &fe_en1,
+          const FiniteElement<dim> &fe_en2,
+          const Quadrature<dim>    &volume_quad,
+          const Quadrature<dim-1>  &face_quad,
+          const bool                distort)
+{
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler_enriched(triangulation);
+  DoFHandler<dim> dof_handler_system(triangulation);
+
+  Point<dim> p1, p2, p3;
+  for (unsigned int d=0; d< dim; d++)
+    {
+      p1[d] = 0.0;
+      p2[d] = 0.0;
+      p3[d] = 0.0;
+    }
+  p2[0] = 10.0;
+  p3[0] = 5.0;
+  EnrichmentFunction<dim> fun1(p1), fun2(p2),fun3(p3);
+  std::vector<const FiniteElement<dim> *> fe_enrichements(2);
+  fe_enrichements[0] = &fe_en1;
+  fe_enrichements[1] = &fe_en2;
+  std::vector<std::vector<std::function<const Function<dim> *(const typename Triangulation<dim, dim>::cell_iterator &) > > > functions(2);
+  functions[0].resize(1);
+  functions[0][0] = [&](const typename Triangulation<dim, dim>::cell_iterator &)
+  {
+    return &fun1;
+  };
+  functions[1].resize(2);
+  functions[1][0] = [&](const typename Triangulation<dim, dim>::cell_iterator &)
+  {
+    return &fun2;
+  };
+  functions[1][1] = [&](const typename Triangulation<dim, dim>::cell_iterator &)
+  {
+    return &fun3;
+  };
+
+  FE_Enriched<dim> fe_enriched(&fe_base,
+                               fe_enrichements,
+                               functions);
+  FESystem<dim> fe_system(fe_base,1,
+                          fe_en1, 1,
+                          fe_en2, 2);
+
+  {
+    Point<dim> p1,p2;
+    std::vector<unsigned int> repetitions(dim);
+    for (unsigned int d=0; d< dim; d++)
+      {
+        p1[d]=-1.0;
+        p2[d]=2.0;
+        repetitions[d] = 3;
+      }
+    GridGenerator::subdivided_hyper_rectangle(triangulation,
+                                              repetitions,
+                                              p1,p2);
+
+    if (distort)
+      GridTools::distort_random (0.1,triangulation);
+  }
+
+  dof_handler_enriched.distribute_dofs(fe_enriched);
+  dof_handler_system.distribute_dofs(fe_system);
+
+  FEValues<dim> fe_values_enriched(fe_enriched, volume_quad,
+                                   update_values | update_gradients | update_hessians );
+
+  FEValues<dim> fe_values_system(fe_system, volume_quad,
+                                 update_values | update_gradients | update_hessians | update_quadrature_points);
+
+  FEFaceValues<dim> fe_face_values_enriched(fe_enriched, face_quad,
+                                            update_values | update_gradients | update_hessians);
+
+  FEFaceValues<dim> fe_face_values_system  (fe_system, face_quad,
+                                            update_values | update_gradients | update_hessians | update_quadrature_points);
+
+  AssertThrow(dof_handler_enriched.n_dofs() == dof_handler_system.n_dofs(),
+              ExcInternalError());
+  Vector<double> solution(dof_handler_enriched.n_dofs());
+
+  const unsigned int dofs_per_cell = fe_enriched.dofs_per_cell;
+  Assert (fe_enriched.dofs_per_cell == fe_system.dofs_per_cell,
+          ExcInternalError());
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell_enriched = dof_handler_enriched.begin_active (),
+  endc_enriched = dof_handler_enriched.end (),
+  cell_system   = dof_handler_system.begin_active(),
+  endc_system   = dof_handler_system.end();
+  for (; cell_enriched!=endc_enriched; ++cell_enriched, ++cell_system)
+    {
+      fe_values_enriched.reinit (cell_enriched);
+      fe_values_system.reinit   (cell_system);
+      const unsigned int n_q_points = volume_quad.size();
+      const std::vector<dealii::Point<dim> > &q_points = fe_values_system.get_quadrature_points();
+
+      // check shape functions
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          check_consistency(q_points[q_point], fun1,fun2,fun3,
+                            fe_values_enriched.shape_value  (i,q_point),
+                            fe_values_enriched.shape_grad   (i,q_point),
+                            fe_values_enriched.shape_hessian(i,q_point),
+                            fe_values_system.shape_value_component   (i,q_point,0),
+                            fe_values_system.shape_grad_component    (i,q_point,0),
+                            fe_values_system.shape_hessian_component (i,q_point,0),
+                            fe_values_system.shape_value_component   (i,q_point,1),
+                            fe_values_system.shape_grad_component    (i,q_point,1),
+                            fe_values_system.shape_hessian_component (i,q_point,1),
+                            fe_values_system.shape_value_component   (i,q_point,2),
+                            fe_values_system.shape_grad_component    (i,q_point,2),
+                            fe_values_system.shape_hessian_component (i,q_point,2),
+                            fe_values_system.shape_value_component   (i,q_point,3),
+                            fe_values_system.shape_grad_component    (i,q_point,3),
+                            fe_values_system.shape_hessian_component (i,q_point,3));
+
+      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+        {
+          fe_face_values_enriched.reinit (cell_enriched,face);
+          fe_face_values_system.reinit (cell_system,face);
+          const unsigned int n_q_points_face = face_quad.size();
+          const std::vector<dealii::Point<dim> > &q_points = fe_face_values_system.get_quadrature_points();
+
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int q_point=0; q_point<n_q_points_face; ++q_point)
+              check_consistency(q_points[q_point], fun1,fun2,fun3,
+                                fe_face_values_enriched.shape_value  (i,q_point),
+                                fe_face_values_enriched.shape_grad   (i,q_point),
+                                fe_face_values_enriched.shape_hessian(i,q_point),
+                                fe_face_values_system.shape_value_component   (i,q_point,0),
+                                fe_face_values_system.shape_grad_component    (i,q_point,0),
+                                fe_face_values_system.shape_hessian_component (i,q_point,0),
+                                fe_face_values_system.shape_value_component   (i,q_point,1),
+                                fe_face_values_system.shape_grad_component    (i,q_point,1),
+                                fe_face_values_system.shape_hessian_component (i,q_point,1),
+                                fe_face_values_system.shape_value_component   (i,q_point,2),
+                                fe_face_values_system.shape_grad_component    (i,q_point,2),
+                                fe_face_values_system.shape_hessian_component (i,q_point,2),
+                                fe_face_values_system.shape_value_component   (i,q_point,3),
+                                fe_face_values_system.shape_grad_component    (i,q_point,3),
+                                fe_face_values_system.shape_hessian_component (i,q_point,3));
+        }
+    }
+
+
+  deallog << "Ok" << std::endl;
+  dof_handler_enriched.clear();
+  dof_handler_system.clear();
+}
+
+
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      {
+        const unsigned int dim = 2;
+        test(FE_Q<dim>(3),FE_Q<dim>(2),FE_Q<dim>(1),QGauss<dim>(2),QGauss<dim-1>(2),false);
+        test(FE_Q<dim>(3),FE_Q<dim>(2),FE_Q<dim>(1),QGauss<dim>(2),QGauss<dim-1>(2),true);
+      }
+
+      {
+        const unsigned int dim = 3;
+        test(FE_Q<dim>(3),FE_Q<dim>(2),FE_Q<dim>(1),QGauss<dim>(2),QGauss<dim-1>(2),false);
+        test(FE_Q<dim>(3),FE_Q<dim>(2),FE_Q<dim>(1),QGauss<dim>(2),QGauss<dim-1>(2),true);
+      }
+
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..f2b98e8
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
diff --git a/tests/fe/fe_enriched_step-36.cc b/tests/fe/fe_enriched_step-36.cc
new file mode 100644 (file)
index 0000000..a22183a
--- /dev/null
@@ -0,0 +1,779 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 FE_Enriched in real-life application on eigenvalue problem similar
+// to Step-36. That involves assembly (shape values and gradients) and
+// error estimator (Kelly - > face gradients) and MPI run.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/slepc_solver.h>
+
+const unsigned int dim = 3;
+
+using namespace dealii;
+
+/**
+ * Coulomb potential
+ */
+template<int dim>
+class PotentialFunction : public dealii::Function<dim>
+{
+public:
+  PotentialFunction()
+    : dealii::Function<dim>(1)
+  {}
+
+  virtual double value(const dealii::Point<dim> &point,
+                       const unsigned int component = 0 ) const;
+};
+
+template<int dim>
+double PotentialFunction<dim>::value(const dealii::Point<dim> &p,
+                                     const unsigned int ) const
+{
+  return -1.0 / std::sqrt(p.square());
+}
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction(const Point<dim> &origin,
+                     const double     &Z,
+                     const double     &radius)
+    : Function<dim>(1),
+      origin(origin),
+      Z(Z),
+      radius(radius)
+  {}
+
+  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 std::exp(-Z*r);
+  }
+
+  bool is_enriched(const Point<dim> &point) const
+  {
+    if (origin.distance(point) < radius)
+      return true;
+    else
+      return false;
+  }
+
+  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();
+    dist/=r;
+    return -Z*std::exp(-Z*r)*dist;
+  }
+
+private:
+  /**
+   * origin
+   */
+  const Point<dim> origin;
+
+  /**
+   * charge
+   */
+  const double Z;
+
+  /**
+   * enrichment radius
+   */
+  const double radius;
+};
+
+namespace Step36
+{
+
+  /**
+   * Main class
+   */
+  template <int dim>
+  class EigenvalueProblem
+  {
+  public:
+    EigenvalueProblem ();
+    virtual ~EigenvalueProblem();
+    void run ();
+
+  private:
+    bool cell_is_pou(const typename hp::DoFHandler<dim>::cell_iterator &cell) const;
+
+    std::pair<unsigned int, unsigned int> setup_system ();
+    void assemble_system ();
+    std::pair<unsigned int, double> solve ();
+    void constrain_pou_dofs();
+    void estimate_error ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>  triangulation;
+    hp::DoFHandler<dim> dof_handler;
+    hp::FECollection<dim> fe_collection;
+    hp::QCollection<dim> q_collection;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
+    std::vector<PETScWrappers::MPI::Vector> eigenfunctions_locally_relevant;
+    std::vector<PetscScalar>                eigenvalues;
+    PETScWrappers::MPI::SparseMatrix        stiffness_matrix, mass_matrix;
+
+    ConstraintMatrix constraints;
+
+    MPI_Comm mpi_communicator;
+    const unsigned int n_mpi_processes;
+    const unsigned int this_mpi_process;
+
+    ConditionalOStream pcout;
+
+    const unsigned int number_of_eigenvalues;
+
+    PotentialFunction<dim> potential;
+
+    EnrichmentFunction<dim> enrichment;
+
+    const FEValuesExtractors::Scalar fe_extractor;
+    const unsigned int               fe_group;
+    const unsigned int               fe_fe_index;
+    const unsigned int               fe_material_id;
+    const FEValuesExtractors::Scalar pou_extractor;
+    const unsigned int               pou_group;
+    const unsigned int               pou_fe_index;
+    const unsigned int               pou_material_id;
+
+    std::vector<Vector<float> > vec_estimated_error_per_cell;
+    Vector<float> estimated_error_per_cell;
+  };
+
+  template <int dim>
+  EigenvalueProblem<dim>::EigenvalueProblem ()
+    :
+    dof_handler (triangulation),
+    mpi_communicator(MPI_COMM_WORLD),
+    n_mpi_processes(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator)),
+    this_mpi_process(dealii::Utilities::MPI::this_mpi_process(mpi_communicator)),
+    pcout (std::cout,
+           (this_mpi_process == 0)),
+    number_of_eigenvalues(1),
+    enrichment(Point<dim>(),
+           /*Z*/1.0,
+           /*radius*/2.5), // radius is set such that 8 cells are marked as enriched
+    fe_extractor(/*dofs start at...*/0),
+    fe_group(/*in FE*/0),
+    fe_fe_index(0),
+    fe_material_id(0),
+    pou_extractor(/*dofs start at (scalar fields!)*/1),
+    pou_group(/*FE*/1),
+    pou_fe_index(1),
+    pou_material_id(1)
+  {
+    dealii::GridGenerator::hyper_cube (triangulation, -10, 10);
+    triangulation.refine_global (2); // 64 cells
+
+    for (auto cell= triangulation.begin_active(); cell != triangulation.end(); ++cell)
+      if (std::sqrt(cell->center().square()) < 5.0)
+        cell->set_refine_flag();
+
+    triangulation.execute_coarsening_and_refinement(); // 120 cells
+
+    q_collection.push_back(QGauss<dim>(4));
+    q_collection.push_back(QGauss<dim>(10));
+
+    // usual elements (active_fe_index ==0):
+    fe_collection.push_back (FE_Enriched<dim> (FE_Q<dim>(2)));
+
+
+    // enriched elements (active_fe_index==1):
+    fe_collection.push_back (FE_Enriched<dim> (FE_Q<dim>(2),
+                                               FE_Q<dim>(1),
+                                               &enrichment));
+
+    // assign fe index in the constructor so that FE/FE+POU is determined
+    // on the coarsest mesh to avoid issues like
+    // +---------+----+----+
+    // |         | fe |    |
+    // |    fe   +----+----+
+    // |         | pou|    |
+    // +---------+----+----+
+    // see discussion in Step46.
+    for (typename hp::DoFHandler<dim>::cell_iterator cell= dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (enrichment.is_enriched(cell->center()))
+        cell->set_material_id(pou_material_id);
+      else
+        cell->set_material_id(fe_material_id);
+  }
+
+  template <int dim>
+  std::pair<unsigned int, unsigned int>
+  EigenvalueProblem<dim>::setup_system ()
+  {
+    for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      {
+        if (cell->material_id() == fe_material_id)
+          cell->set_active_fe_index (fe_fe_index);
+        else if (cell->material_id() == pou_material_id)
+          cell->set_active_fe_index (pou_fe_index);
+        else
+          Assert (false, ExcNotImplemented());
+      }
+
+    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_processor
+      = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
+    locally_owned_dofs = locally_owned_dofs_per_processor[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);
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                              0,
+                                              ZeroFunction<dim> (1),
+                                              constraints);
+    constraints.close ();
+
+    // Initialise the stiffness and mass matrices
+    DynamicSparsityPattern csp (locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints,
+                                     /* keep constrained dofs */ 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_processor[i].n_elements();
+
+    SparsityTools::distribute_sparsity_pattern
+    (csp,
+     n_locally_owned_dofs,
+     mpi_communicator,
+     locally_relevant_dofs);
+
+    stiffness_matrix.reinit (locally_owned_dofs,
+                             locally_owned_dofs,
+                             csp,
+                             mpi_communicator);
+
+    mass_matrix.reinit (locally_owned_dofs,
+                        locally_owned_dofs,
+                        csp,
+                        mpi_communicator);
+
+    // reinit vectors
+    eigenfunctions.resize (number_of_eigenvalues);
+    eigenfunctions_locally_relevant.resize(number_of_eigenvalues);
+    vec_estimated_error_per_cell.resize(number_of_eigenvalues);
+    for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+      {
+        eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+        eigenfunctions_locally_relevant[i].reinit (locally_owned_dofs,
+                                                   locally_relevant_dofs,
+                                                   mpi_communicator);
+
+        vec_estimated_error_per_cell[i].reinit(triangulation.n_active_cells());
+      }
+
+    eigenvalues.resize (eigenfunctions.size ());
+
+    estimated_error_per_cell.reinit (triangulation.n_active_cells());
+
+    unsigned int n_pou_cells = 0,
+                 n_fem_cells = 0;
+
+    for (typename hp::DoFHandler<dim>::cell_iterator cell= dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell_is_pou(cell))
+        n_pou_cells++;
+      else
+        n_fem_cells++;
+
+    return std::make_pair (n_fem_cells,
+                           n_pou_cells);
+  }
+
+  template <int dim>
+  bool
+  EigenvalueProblem<dim>::cell_is_pou(const typename hp::DoFHandler<dim>::cell_iterator &cell) const
+  {
+    return cell->material_id() == pou_material_id;
+  }
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::constrain_pou_dofs()
+  {
+    std::vector<unsigned int> local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face);
+    for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell_is_pou(cell))
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (!cell->at_boundary(f))
+            {
+              bool face_is_on_interface = false;
+              if ((cell->neighbor(f)->has_children() == false /* => it is active */)
+                  &&
+                  (! cell_is_pou(cell->neighbor(f))))
+                face_is_on_interface = true;
+              else if (cell->neighbor(f)->has_children() == true)
+                {
+                  // The neighbor does
+                  // have
+                  // children. See if
+                  // any of the cells
+                  // on the other
+                  // side are vanilla FEM
+                  for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
+                    if (!cell_is_pou (cell->neighbor_child_on_subface(f, sf)))
+                      {
+                        face_is_on_interface = true;
+                        break;
+                      }
+                }
+              // add constraints
+              if (face_is_on_interface)
+                {
+                  cell->face(f)->get_dof_indices (local_face_dof_indices, pou_fe_index);
+                  for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
+                    if (fe_collection[pou_fe_index].face_system_to_base_index(i).first.first == pou_group)
+                      //if (fe_collection[1].face_system_to_component_index(i).first /*component*/ > 0)
+                      constraints.add_line (local_face_dof_indices[i]);
+                }
+            }
+  }
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::assemble_system()
+  {
+    stiffness_matrix = 0;
+    mass_matrix = 0;
+
+    dealii::FullMatrix<double> cell_stiffness_matrix;
+    dealii::FullMatrix<double> cell_mass_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<double> potential_values;
+    std::vector<double> enrichment_values;
+    std::vector<Tensor<1,dim>> enrichment_gradients;
+
+    hp::FEValues<dim> fe_values_hp(fe_collection, q_collection,
+                                   update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values);
+
+
+    typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active (),
+    endc = dof_handler.end ();
+    for (; cell!=endc; ++cell)
+      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;
+
+          potential_values.resize(n_q_points);
+          enrichment_values.resize(n_q_points);
+          enrichment_gradients.resize(n_q_points);
+
+          local_dof_indices.resize     (dofs_per_cell);
+          cell_stiffness_matrix.reinit (dofs_per_cell,dofs_per_cell);
+          cell_mass_matrix.reinit      (dofs_per_cell,dofs_per_cell);
+
+          cell_stiffness_matrix = 0;
+          cell_mass_matrix      = 0;
+
+          potential.value_list (fe_values.get_quadrature_points(),
+                                potential_values);
+
+          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_stiffness_matrix (i, j)
+                  += (fe_values[fe_extractor].gradient (i, q_point) *
+                      fe_values[fe_extractor].gradient (j, q_point) *
+                      0.5
+                      +
+                      potential_values[q_point] *
+                      fe_values[fe_extractor].value (i, q_point) *
+                      fe_values[fe_extractor].value (j, q_point)
+                     ) * fe_values.JxW (q_point);
+
+                  cell_mass_matrix (i, j)
+                  += (fe_values[fe_extractor].value (i, q_point) *
+                      fe_values[fe_extractor].value (j, 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_stiffness_matrix (j, i) = cell_stiffness_matrix (i, j);
+                cell_mass_matrix (j, i)      = cell_mass_matrix (i, j);
+              }
+
+          cell->get_dof_indices (local_dof_indices);
+
+          constraints
+          .distribute_local_to_global (cell_stiffness_matrix,
+                                       local_dof_indices,
+                                       stiffness_matrix);
+          constraints
+          .distribute_local_to_global (cell_mass_matrix,
+                                       local_dof_indices,
+                                       mass_matrix);
+        }
+
+    stiffness_matrix.compress (dealii::VectorOperation::add);
+    mass_matrix.compress (dealii::VectorOperation::add);
+  }
+
+  template<int dim>
+  std::pair<unsigned int, double>
+  EigenvalueProblem<dim>::solve()
+  {
+    dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false);
+    dealii::SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control,
+                                                          mpi_communicator);
+
+    eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);
+    eigensolver.set_problem_type (EPS_GHEP);
+
+    eigensolver.solve (stiffness_matrix, mass_matrix,
+                       eigenvalues, eigenfunctions,
+                       eigenfunctions.size());
+
+    for (unsigned int i = 0; i < eigenfunctions.size(); i++)
+      {
+        constraints.distribute(eigenfunctions[i]);
+        eigenfunctions_locally_relevant[i] = eigenfunctions[i];
+      }
+
+    return std::make_pair (solver_control.last_step (),
+                           solver_control.last_value ());
+
+  }
+
+  template <int dim>
+  EigenvalueProblem<dim>::~EigenvalueProblem()
+  {
+    dof_handler.clear ();
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::estimate_error ()
+  {
+    {
+      std::vector<const PETScWrappers::MPI::Vector *> sol (number_of_eigenvalues);
+      std::vector<dealii::Vector<float> *>   error (number_of_eigenvalues);
+
+      for (unsigned int i = 0; i < number_of_eigenvalues; i++)
+        {
+          sol[i] = &eigenfunctions_locally_relevant[i];
+          error[i] = &vec_estimated_error_per_cell[i];
+        }
+
+      hp::QCollection<dim-1> face_quadrature_formula;
+      face_quadrature_formula.push_back (dealii::QGauss<dim-1>(3));
+      face_quadrature_formula.push_back (dealii::QGauss<dim-1>(3));
+
+      KellyErrorEstimator<dim>::estimate (dof_handler,
+                                          face_quadrature_formula,
+                                          typename FunctionMap<dim>::type(),
+                                          sol,
+                                          error);
+    }
+
+    // sum up for a global:
+    for (unsigned int c=0; c < estimated_error_per_cell.size(); c++)
+      {
+        double er = 0.0;
+        for (unsigned int i = 0; i < number_of_eigenvalues; i++)
+          er += vec_estimated_error_per_cell[i][c] * vec_estimated_error_per_cell[i][c];
+
+        estimated_error_per_cell[c] = sqrt(er);
+      }
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::refine_grid ()
+  {
+    const double threshold = 0.9 * estimated_error_per_cell.linfty_norm();
+    GridRefinement::refine (triangulation,
+                            estimated_error_per_cell,
+                            threshold);
+
+    triangulation.prepare_coarsening_and_refinement ();
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+  template <int dim>
+  class Postprocessor : public DataPostprocessorScalar<dim>
+  {
+  public:
+    Postprocessor(const EnrichmentFunction<dim> &enrichment);
+
+    virtual
+    void
+    compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                       const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                       const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                       const std::vector<Point<dim> >                  &normals,
+                                       const std::vector<Point<dim> >                  &evaluation_points,
+                                       std::vector<Vector<double> >                    &computed_quantities) const;
+
+  private:
+    const EnrichmentFunction<dim> &enrichment;
+  };
+
+  template <int dim>
+  Postprocessor<dim>::Postprocessor(const EnrichmentFunction<dim> &enrichment)
+    :
+    DataPostprocessorScalar<dim> ("total_solution",
+                                  update_values | update_q_points),
+    enrichment(enrichment)
+  {}
+
+  template <int dim>
+  void
+  Postprocessor<dim>::
+  compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                     const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
+                                     const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+                                     const std::vector<Point<dim> >                  &/*normals*/,
+                                     const std::vector<Point<dim> >                  &evaluation_points,
+                                     std::vector<Vector<double> >                    &computed_quantities) const
+  {
+    const unsigned int n_quadrature_points = uh.size();
+    Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
+    for (unsigned int q=0; q<n_quadrature_points; ++q)
+      {
+        Assert(uh[q].size() == 2, ExcDimensionMismatch (uh[q].size(), 2)); // FESystem with 2 components
+
+        computed_quantities[q](0)
+          = (uh[q](0)
+             +
+             uh[q](1) * enrichment.value(evaluation_points[q])); // for FE_Nothing uh[q](1) will be zero
+      }
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    dealii::Vector<float> fe_index(triangulation.n_active_cells());
+    {
+      typename dealii::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();
+    }
+
+    Assert (cycle < 10, ExcNotImplemented());
+    if (this_mpi_process==0)
+      {
+        std::string filename = "solution-";
+        filename += ('0' + cycle);
+        filename += ".vtk";
+        std::ofstream output (filename.c_str());
+
+        Postprocessor<dim> postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration
+        DataOut<dim,hp::DoFHandler<dim> > data_out;
+        data_out.attach_dof_handler (dof_handler);
+        data_out.add_data_vector (eigenfunctions_locally_relevant[0], "solution");
+        data_out.build_patches (6);
+        data_out.write_vtk (output);
+      }
+
+    // second output without making mesh finer
+    if (this_mpi_process==0)
+      {
+        std::string filename = "mesh-";
+        filename += ('0' + 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 (fe_index, "fe_index");
+        data_out.add_data_vector (estimated_error_per_cell, "estimated_error");
+        data_out.build_patches ();
+        data_out.write_vtk (output);
+      }
+
+    // scalar data for plotting
+    //output scalar data (eigenvalues, energies, ndofs, etc)
+    if (this_mpi_process == 0)
+      {
+        const std::string scalar_fname = "scalar-data.txt";
+
+        std::ofstream output (scalar_fname.c_str(),
+                              std::ios::out |
+                              (cycle==0
+                               ?
+                               std::ios::trunc : std::ios::app));
+
+        std::streamsize max_precision  = std::numeric_limits<double>::digits10;
+
+        std::string sep("  ");
+        output << cycle << sep
+               << std::setprecision(max_precision)
+               << triangulation.n_active_cells() << sep
+               << dof_handler.n_dofs () << sep
+               << std::scientific;
+
+        for (unsigned int i = 0; i < eigenvalues.size(); i++)
+          output << eigenvalues[i] << sep;
+
+        output<<std::endl;
+        output.close();
+      } //end scope
+
+  }
+
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::run()
+  {
+    for (unsigned int cycle = 0; cycle < 5; cycle++)
+      {
+        pcout << "Cycle "<<cycle <<std::endl;
+        const std::pair<unsigned int, unsigned int> n_cells = setup_system ();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_active_cells ()
+              << std::endl
+              << "     FE / POUFE :                "
+              << n_cells.first << " / " << n_cells.second
+              << std::endl
+              << "   Number of degrees of freedom: "
+              << dof_handler.n_dofs ()
+              << std::endl;
+
+        assemble_system ();
+
+        const std::pair<unsigned int, double> res = solve ();
+        AssertThrow(res.second < 5e-8,
+                    ExcInternalError());
+
+        estimate_error ();
+        // output_results(cycle);
+        refine_grid ();
+
+        pcout << std::endl;
+        for (unsigned int i=0; i<eigenvalues.size(); ++i)
+          pcout << "      Eigenvalue " << i
+                << " : " << eigenvalues[i]
+                << std::endl;
+      }
+
+  }
+}
+
+int main (int argc,char **argv)
+{
+
+  try
+    {
+      dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        Step36::EigenvalueProblem<dim> step36;
+        dealii::PETScWrappers::set_option_value("-eps_target","-1.0");
+        dealii::PETScWrappers::set_option_value("-st_type","sinvert");
+        dealii::PETScWrappers::set_option_value("-st_ksp_type","cg");
+        dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg");
+        dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11");
+        step36.run();
+      }
+
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..fa33870
--- /dev/null
@@ -0,0 +1,30 @@
+Cycle 0
+   Number of active cells:       120
+     FE / POUFE :                120 / 8
+   Number of degrees of freedom: 1432
+
+      Eigenvalue 0 : -0.483398
+Cycle 1
+   Number of active cells:       176
+     FE / POUFE :                120 / 72
+   Number of degrees of freedom: 2206
+
+      Eigenvalue 0 : -0.498786
+Cycle 2
+   Number of active cells:       512
+     FE / POUFE :                504 / 72
+   Number of degrees of freedom: 5734
+
+      Eigenvalue 0 : -0.499788
+Cycle 3
+   Number of active cells:       568
+     FE / POUFE :                504 / 136
+   Number of degrees of freedom: 6508
+
+      Eigenvalue 0 : -0.499818
+Cycle 4
+   Number of active cells:       736
+     FE / POUFE :                504 / 328
+   Number of degrees of freedom: 8542
+
+      Eigenvalue 0 : -0.499947
diff --git a/tests/fe/fe_enriched_step-36b.cc b/tests/fe/fe_enriched_step-36b.cc
new file mode 100644 (file)
index 0000000..4a92b81
--- /dev/null
@@ -0,0 +1,903 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 FE_Enriched in real-life application on eigenvalue problem similar
+// to Step-36. That involves assembly (shape values and gradients) and
+// error estimator (Kelly - > face gradients) and MPI run.
+//
+// same as fe_enriched_step-36, but using FESystem only. Note that the results
+// of KellyErrorEstimator won't be the same, thus this tests does only a single
+// refinement cycle where the refinement path is the same.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/slepc_solver.h>
+
+const unsigned int dim = 3;
+
+using namespace dealii;
+
+/**
+ * Coulomb potential
+ */
+template<int dim>
+class PotentialFunction : public dealii::Function<dim>
+{
+public:
+  PotentialFunction()
+    : dealii::Function<dim>(1)
+  {}
+
+  virtual double value(const dealii::Point<dim> &point,
+                       const unsigned int component = 0 ) const;
+};
+
+template<int dim>
+double PotentialFunction<dim>::value(const dealii::Point<dim> &p,
+                                     const unsigned int ) const
+{
+  return -1.0 / std::sqrt(p.square());
+}
+
+template<int dim>
+class EnrichmentFunction : public Function<dim>
+{
+public:
+  EnrichmentFunction(const Point<dim> &origin,
+                     const double     &Z,
+                     const double     &radius)
+    : Function<dim>(1),
+      origin(origin),
+      Z(Z),
+      radius(radius)
+  {}
+
+  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 std::exp(-Z*r);
+  }
+
+  bool is_enriched(const Point<dim> &point) const
+  {
+    if (origin.distance(point) < radius)
+      return true;
+    else
+      return false;
+  }
+
+  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();
+    dist/=r;
+    return -Z*std::exp(-Z*r)*dist;
+  }
+
+private:
+  /**
+   * origin
+   */
+  const Point<dim> origin;
+
+  /**
+   * charge
+   */
+  const double Z;
+
+  /**
+   * enrichment radius
+   */
+  const double radius;
+};
+
+namespace Step36
+{
+
+  /**
+   * Main class
+   */
+  template <int dim>
+  class EigenvalueProblem
+  {
+  public:
+    EigenvalueProblem ();
+    virtual ~EigenvalueProblem();
+    void run ();
+
+  private:
+    bool cell_is_pou(const typename hp::DoFHandler<dim>::cell_iterator &cell) const;
+
+    std::pair<unsigned int, unsigned int> setup_system ();
+    void assemble_system ();
+    std::pair<unsigned int, double> solve ();
+    void constrain_pou_dofs();
+    void estimate_error ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>  triangulation;
+    hp::DoFHandler<dim> dof_handler;
+    hp::FECollection<dim> fe_collection;
+    hp::QCollection<dim> q_collection;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
+    std::vector<PETScWrappers::MPI::Vector> eigenfunctions_locally_relevant;
+    std::vector<PetscScalar>                eigenvalues;
+    PETScWrappers::MPI::SparseMatrix        stiffness_matrix, mass_matrix;
+
+    ConstraintMatrix constraints;
+
+    MPI_Comm mpi_communicator;
+    const unsigned int n_mpi_processes;
+    const unsigned int this_mpi_process;
+
+    ConditionalOStream pcout;
+
+    const unsigned int number_of_eigenvalues;
+
+    PotentialFunction<dim> potential;
+
+    EnrichmentFunction<dim> enrichment;
+
+    const FEValuesExtractors::Scalar fe_extractor;
+    const unsigned int               fe_group;
+    const unsigned int               fe_fe_index;
+    const unsigned int               fe_material_id;
+    const FEValuesExtractors::Scalar pou_extractor;
+    const unsigned int               pou_group;
+    const unsigned int               pou_fe_index;
+    const unsigned int               pou_material_id;
+
+    std::vector<Vector<float> > vec_estimated_error_per_cell;
+    Vector<float> estimated_error_per_cell;
+  };
+
+  template <int dim>
+  EigenvalueProblem<dim>::EigenvalueProblem ()
+    :
+    dof_handler (triangulation),
+    mpi_communicator(MPI_COMM_WORLD),
+    n_mpi_processes(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator)),
+    this_mpi_process(dealii::Utilities::MPI::this_mpi_process(mpi_communicator)),
+    pcout (std::cout,
+           (this_mpi_process == 0)),
+    number_of_eigenvalues(1),
+    enrichment(Point<dim>(),
+           /*Z*/1.0,
+           /*radius*/2.5), // radius is set such that 8 cells are marked as enriched
+    fe_extractor(/*dofs start at...*/0),
+    fe_group(/*in FE*/0),
+    fe_fe_index(0),
+    fe_material_id(0),
+    pou_extractor(/*dofs start at (scalar fields!)*/1),
+    pou_group(/*FE*/1),
+    pou_fe_index(1),
+    pou_material_id(1)
+  {
+    dealii::GridGenerator::hyper_cube (triangulation, -10, 10);
+    triangulation.refine_global (2); // 64 cells
+
+    for (auto cell= triangulation.begin_active(); cell != triangulation.end(); ++cell)
+      if (std::sqrt(cell->center().square()) < 5.0)
+        cell->set_refine_flag();
+
+    triangulation.execute_coarsening_and_refinement(); // 120 cells
+
+    q_collection.push_back(QGauss<dim>(4));
+    q_collection.push_back(QGauss<dim>(10));
+
+    // usual elements (active_fe_index ==0):
+    fe_collection.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                            FE_Nothing<dim>(), 1));
+
+    // enriched elements (active_fe_index==1):
+    fe_collection.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                            FE_Q<dim>(1), 1));
+
+    // assign fe index in the constructor so that FE/FE+POU is determined
+    // on the coarsest mesh to avoid issues like
+    // +---------+----+----+
+    // |         | fe |    |
+    // |    fe   +----+----+
+    // |         | pou|    |
+    // +---------+----+----+
+    // see discussion in Step46.
+    for (typename hp::DoFHandler<dim>::cell_iterator cell= dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (enrichment.is_enriched(cell->center()))
+        cell->set_material_id(pou_material_id);
+      else
+        cell->set_material_id(fe_material_id);
+  }
+
+  template <int dim>
+  std::pair<unsigned int, unsigned int>
+  EigenvalueProblem<dim>::setup_system ()
+  {
+    for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      {
+        if (cell->material_id() == fe_material_id)
+          cell->set_active_fe_index (fe_fe_index);
+        else if (cell->material_id() == pou_material_id)
+          cell->set_active_fe_index (pou_fe_index);
+        else
+          Assert (false, ExcNotImplemented());
+      }
+
+    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_processor
+      = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
+    locally_owned_dofs = locally_owned_dofs_per_processor[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);
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                              0,
+                                              ZeroFunction<dim> (2),
+                                              constraints);
+    constrain_pou_dofs();
+    constraints.close ();
+
+    // Initialise the stiffness and mass matrices
+    DynamicSparsityPattern csp (locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints,
+                                     /* keep constrained dofs */ 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_processor[i].n_elements();
+
+    SparsityTools::distribute_sparsity_pattern
+    (csp,
+     n_locally_owned_dofs,
+     mpi_communicator,
+     locally_relevant_dofs);
+
+    stiffness_matrix.reinit (locally_owned_dofs,
+                             locally_owned_dofs,
+                             csp,
+                             mpi_communicator);
+
+    mass_matrix.reinit (locally_owned_dofs,
+                        locally_owned_dofs,
+                        csp,
+                        mpi_communicator);
+
+    // reinit vectors
+    eigenfunctions.resize (number_of_eigenvalues);
+    eigenfunctions_locally_relevant.resize(number_of_eigenvalues);
+    vec_estimated_error_per_cell.resize(number_of_eigenvalues);
+    for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+      {
+        eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+        eigenfunctions_locally_relevant[i].reinit (locally_owned_dofs,
+                                                   locally_relevant_dofs,
+                                                   mpi_communicator);
+
+        vec_estimated_error_per_cell[i].reinit(triangulation.n_active_cells());
+      }
+
+    eigenvalues.resize (eigenfunctions.size ());
+
+    estimated_error_per_cell.reinit (triangulation.n_active_cells());
+
+    unsigned int n_pou_cells = 0,
+                 n_fem_cells = 0;
+
+    for (typename hp::DoFHandler<dim>::cell_iterator cell= dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell_is_pou(cell))
+        n_pou_cells++;
+      else
+        n_fem_cells++;
+
+    return std::make_pair (n_fem_cells,
+                           n_pou_cells);
+  }
+
+  template <int dim>
+  bool
+  EigenvalueProblem<dim>::cell_is_pou(const typename hp::DoFHandler<dim>::cell_iterator &cell) const
+  {
+    return cell->material_id() == pou_material_id;
+  }
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::constrain_pou_dofs()
+  {
+    std::vector<unsigned int> local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face);
+    for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell_is_pou(cell))
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (!cell->at_boundary(f))
+            {
+              bool face_is_on_interface = false;
+              if ((cell->neighbor(f)->has_children() == false /* => it is active */)
+                  &&
+                  (! cell_is_pou(cell->neighbor(f))))
+                face_is_on_interface = true;
+              else if (cell->neighbor(f)->has_children() == true)
+                {
+                  // The neighbor does
+                  // have
+                  // children. See if
+                  // any of the cells
+                  // on the other
+                  // side are vanilla FEM
+                  for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
+                    if (!cell_is_pou (cell->neighbor_child_on_subface(f, sf)))
+                      {
+                        face_is_on_interface = true;
+                        break;
+                      }
+                }
+              // add constraints
+              if (face_is_on_interface)
+                {
+                  cell->face(f)->get_dof_indices (local_face_dof_indices, pou_fe_index);
+                  for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
+                    if (fe_collection[pou_fe_index].face_system_to_base_index(i).first.first == pou_group)
+                      //if (fe_collection[1].face_system_to_component_index(i).first /*component*/ > 0)
+                      constraints.add_line (local_face_dof_indices[i]);
+                }
+            }
+  }
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::assemble_system()
+  {
+    stiffness_matrix = 0;
+    mass_matrix = 0;
+
+    dealii::FullMatrix<double> cell_stiffness_matrix;
+    dealii::FullMatrix<double> cell_mass_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<double> potential_values;
+    std::vector<double> enrichment_values;
+    std::vector<Tensor<1,dim>> enrichment_gradients;
+
+    hp::FEValues<dim> fe_values_hp(fe_collection, q_collection,
+                                   update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values);
+
+
+    typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active (),
+    endc = dof_handler.end ();
+    for (; cell!=endc; ++cell)
+      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;
+
+          potential_values.resize(n_q_points);
+          enrichment_values.resize(n_q_points);
+          enrichment_gradients.resize(n_q_points);
+
+          local_dof_indices.resize     (dofs_per_cell);
+          cell_stiffness_matrix.reinit (dofs_per_cell,dofs_per_cell);
+          cell_mass_matrix.reinit      (dofs_per_cell,dofs_per_cell);
+
+          cell_stiffness_matrix = 0;
+          cell_mass_matrix      = 0;
+
+          potential.value_list (fe_values.get_quadrature_points(),
+                                potential_values);
+
+          if (cell->active_fe_index() == 0) // plain FE
+            {
+              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_stiffness_matrix (i, j)
+                      += (fe_values[fe_extractor].gradient (i, q_point) *
+                          fe_values[fe_extractor].gradient (j, q_point) *
+                          0.5
+                          +
+                          potential_values[q_point] *
+                          fe_values[fe_extractor].value (i, q_point) *
+                          fe_values[fe_extractor].value (j, q_point)
+                         ) * fe_values.JxW (q_point);
+
+                      cell_mass_matrix (i, j)
+                      += (fe_values[fe_extractor].value (i, q_point) *
+                          fe_values[fe_extractor].value (j, q_point)
+                         ) * fe_values.JxW (q_point);
+                    }
+            }
+          else // POUFEM
+            {
+              Assert (cell->active_fe_index() == 1, ExcInternalError());
+
+              enrichment.value_list (fe_values.get_quadrature_points(),
+                                     enrichment_values);
+              enrichment.gradient_list (fe_values.get_quadrature_points(),
+                                        enrichment_gradients);
+
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                  const unsigned int i_group = cell->get_fe().system_to_base_index(i).first.first;
+                  for (unsigned int j=i; j<dofs_per_cell; ++j)
+                    {
+                      const unsigned int j_group     = cell->get_fe().system_to_base_index(j).first.first;
+
+                      if ( i_group == fe_group  && j_group == fe_group ) // fe - fe
+                        {
+                          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                            {
+                              cell_stiffness_matrix (i, j)
+                              += (fe_values[fe_extractor].gradient (i, q_point) *
+                                  fe_values[fe_extractor].gradient (j, q_point) *
+                                  0.5
+                                  +
+                                  potential_values[q_point] *
+                                  fe_values[fe_extractor].value (i, q_point) *
+                                  fe_values[fe_extractor].value (j, q_point)
+                                 ) * fe_values.JxW (q_point);
+
+                              cell_mass_matrix (i, j)
+                              += (fe_values[fe_extractor].value (i, q_point) *
+                                  fe_values[fe_extractor].value (j, q_point)
+                                 ) * fe_values.JxW (q_point);
+                            }
+                        }
+                      else if ( i_group == fe_group  && j_group == pou_group ) // fe - pou
+                        {
+                          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                            {
+                              cell_stiffness_matrix (i, j)
+                              += (fe_values [fe_extractor].gradient  (i, q_point) *
+                                  (fe_values[pou_extractor].gradient (j, q_point) * enrichment_values   [q_point] +
+                                   fe_values[pou_extractor].value    (j, q_point) * enrichment_gradients[q_point]) *
+                                  0.5
+                                  +
+                                  potential_values[q_point] *
+                                  fe_values[fe_extractor].value  (i, q_point) *
+                                  fe_values[pou_extractor].value (j, q_point) *
+                                  enrichment_values[q_point]
+                                 ) * fe_values.JxW (q_point);
+
+                              cell_mass_matrix (i, j)
+                              += (fe_values[fe_extractor].value (i, q_point) *
+                                  fe_values[pou_extractor].value (j, q_point) *
+                                  enrichment_values[q_point]
+                                 ) * fe_values.JxW (q_point);
+                            }
+                        }
+                      else if ( i_group == pou_group  && j_group == fe_group ) // pou - fe
+                        {
+                          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                            {
+                              cell_stiffness_matrix (i, j)
+                              += ((fe_values[pou_extractor].gradient (i, q_point) * enrichment_values   [q_point] +
+                                   fe_values[pou_extractor].value    (i, q_point) * enrichment_gradients[q_point]) *
+                                  fe_values [fe_extractor].gradient  (j, q_point) *
+                                  0.5
+                                  +
+                                  potential_values[q_point] *
+                                  fe_values[pou_extractor].value (i, q_point) *
+                                  enrichment_values[q_point] *
+                                  fe_values[fe_extractor].value (j, q_point)
+                                 ) * fe_values.JxW (q_point);
+
+                              cell_mass_matrix (i, j)
+                              += (fe_values[pou_extractor].value (i, q_point) *
+                                  enrichment_values[q_point] *
+                                  fe_values[fe_extractor].value (j, q_point)
+                                 ) * fe_values.JxW (q_point);
+                            }
+                        }
+                      else // pou - pou
+                        {
+                          Assert (i_group == pou_group  && j_group == pou_group,
+                                  ExcInternalError());
+
+                          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                            {
+                              cell_stiffness_matrix (i, j)
+                              += ((fe_values[pou_extractor].gradient (i, q_point) * enrichment_values   [q_point] +
+                                   fe_values[pou_extractor].value    (i, q_point) * enrichment_gradients[q_point]) *
+                                  (fe_values[pou_extractor].gradient (j, q_point) * enrichment_values   [q_point] +
+                                   fe_values[pou_extractor].value    (j, q_point) * enrichment_gradients[q_point]) *
+                                  0.5
+                                  +
+                                  potential_values[q_point] *
+                                  fe_values[pou_extractor].value (i, q_point) *
+                                  enrichment_values[q_point] *
+                                  fe_values[pou_extractor].value (j, q_point) *
+                                  enrichment_values[q_point]
+                                 ) * fe_values.JxW (q_point);
+
+                              cell_mass_matrix (i, j)
+                              += (fe_values[pou_extractor].value (i, q_point) *
+                                  enrichment_values[q_point] *
+                                  fe_values[pou_extractor].value (j, q_point) *
+                                  enrichment_values[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_stiffness_matrix (j, i) = cell_stiffness_matrix (i, j);
+                cell_mass_matrix (j, i)      = cell_mass_matrix (i, j);
+              }
+
+          cell->get_dof_indices (local_dof_indices);
+
+          constraints
+          .distribute_local_to_global (cell_stiffness_matrix,
+                                       local_dof_indices,
+                                       stiffness_matrix);
+          constraints
+          .distribute_local_to_global (cell_mass_matrix,
+                                       local_dof_indices,
+                                       mass_matrix);
+        }
+
+    stiffness_matrix.compress (dealii::VectorOperation::add);
+    mass_matrix.compress (dealii::VectorOperation::add);
+  }
+
+  template<int dim>
+  std::pair<unsigned int, double>
+  EigenvalueProblem<dim>::solve()
+  {
+    dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false);
+    dealii::SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control,
+                                                          mpi_communicator);
+
+    eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);
+    eigensolver.set_problem_type (EPS_GHEP);
+
+    eigensolver.solve (stiffness_matrix, mass_matrix,
+                       eigenvalues, eigenfunctions,
+                       eigenfunctions.size());
+
+    for (unsigned int i = 0; i < eigenfunctions.size(); i++)
+      {
+        constraints.distribute(eigenfunctions[i]);
+        eigenfunctions_locally_relevant[i] = eigenfunctions[i];
+      }
+
+    return std::make_pair (solver_control.last_step (),
+                           solver_control.last_value ());
+
+  }
+
+  template <int dim>
+  EigenvalueProblem<dim>::~EigenvalueProblem()
+  {
+    dof_handler.clear ();
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::estimate_error ()
+  {
+    {
+      std::vector<const PETScWrappers::MPI::Vector *> sol (number_of_eigenvalues);
+      std::vector<dealii::Vector<float> *>   error (number_of_eigenvalues);
+
+      for (unsigned int i = 0; i < number_of_eigenvalues; i++)
+        {
+          sol[i] = &eigenfunctions_locally_relevant[i];
+          error[i] = &vec_estimated_error_per_cell[i];
+        }
+
+      hp::QCollection<dim-1> face_quadrature_formula;
+      face_quadrature_formula.push_back (dealii::QGauss<dim-1>(3));
+      face_quadrature_formula.push_back (dealii::QGauss<dim-1>(3));
+
+      KellyErrorEstimator<dim>::estimate (dof_handler,
+                                          face_quadrature_formula,
+                                          typename FunctionMap<dim>::type(),
+                                          sol,
+                                          error);
+    }
+
+    // sum up for a global:
+    for (unsigned int c=0; c < estimated_error_per_cell.size(); c++)
+      {
+        double er = 0.0;
+        for (unsigned int i = 0; i < number_of_eigenvalues; i++)
+          er += vec_estimated_error_per_cell[i][c] * vec_estimated_error_per_cell[i][c];
+
+        estimated_error_per_cell[c] = sqrt(er);
+      }
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::refine_grid ()
+  {
+    const double threshold = 0.9 * estimated_error_per_cell.linfty_norm();
+    GridRefinement::refine (triangulation,
+                            estimated_error_per_cell,
+                            threshold);
+
+    triangulation.prepare_coarsening_and_refinement ();
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+  template <int dim>
+  class Postprocessor : public DataPostprocessorScalar<dim>
+  {
+  public:
+    Postprocessor(const EnrichmentFunction<dim> &enrichment);
+
+    virtual
+    void
+    compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                       const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                       const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                       const std::vector<Point<dim> >                  &normals,
+                                       const std::vector<Point<dim> >                  &evaluation_points,
+                                       std::vector<Vector<double> >                    &computed_quantities) const;
+
+  private:
+    const EnrichmentFunction<dim> &enrichment;
+  };
+
+  template <int dim>
+  Postprocessor<dim>::Postprocessor(const EnrichmentFunction<dim> &enrichment)
+    :
+    DataPostprocessorScalar<dim> ("total_solution",
+                                  update_values | update_q_points),
+    enrichment(enrichment)
+  {}
+
+  template <int dim>
+  void
+  Postprocessor<dim>::
+  compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                     const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
+                                     const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+                                     const std::vector<Point<dim> >                  &/*normals*/,
+                                     const std::vector<Point<dim> >                  &evaluation_points,
+                                     std::vector<Vector<double> >                    &computed_quantities) const
+  {
+    const unsigned int n_quadrature_points = uh.size();
+    Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
+    for (unsigned int q=0; q<n_quadrature_points; ++q)
+      {
+        Assert(uh[q].size() == 2, ExcDimensionMismatch (uh[q].size(), 2)); // FESystem with 2 components
+
+        computed_quantities[q](0)
+          = (uh[q](0)
+             +
+             uh[q](1) * enrichment.value(evaluation_points[q])); // for FE_Nothing uh[q](1) will be zero
+      }
+  }
+
+  template <int dim>
+  void EigenvalueProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    dealii::Vector<float> fe_index(triangulation.n_active_cells());
+    {
+      typename dealii::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();
+    }
+
+    Assert (cycle < 10, ExcNotImplemented());
+    if (this_mpi_process==0)
+      {
+        std::string filename = "solution-";
+        filename += ('0' + cycle);
+        filename += ".vtk";
+        std::ofstream output (filename.c_str());
+
+        Postprocessor<dim> postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration
+        DataOut<dim,hp::DoFHandler<dim> > data_out;
+        data_out.attach_dof_handler (dof_handler);
+        data_out.add_data_vector (eigenfunctions_locally_relevant[0], "solution");
+        data_out.add_data_vector (eigenfunctions_locally_relevant[0], postprocessor);
+        data_out.build_patches (6);
+        data_out.write_vtk (output);
+      }
+
+    // second output without making mesh finer
+    if (this_mpi_process==0)
+      {
+        std::string filename = "mesh-";
+        filename += ('0' + 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 (fe_index, "fe_index");
+        data_out.add_data_vector (estimated_error_per_cell, "estimated_error");
+        data_out.build_patches ();
+        data_out.write_vtk (output);
+      }
+
+    // scalar data for plotting
+    //output scalar data (eigenvalues, energies, ndofs, etc)
+    if (this_mpi_process == 0)
+      {
+        const std::string scalar_fname = "scalar-data.txt";
+
+        std::ofstream output (scalar_fname.c_str(),
+                              std::ios::out |
+                              (cycle==0
+                               ?
+                               std::ios::trunc : std::ios::app));
+
+        std::streamsize max_precision  = std::numeric_limits<double>::digits10;
+
+        std::string sep("  ");
+        output << cycle << sep
+               << std::setprecision(max_precision)
+               << triangulation.n_active_cells() << sep
+               << dof_handler.n_dofs () << sep
+               << std::scientific;
+
+        for (unsigned int i = 0; i < eigenvalues.size(); i++)
+          output << eigenvalues[i] << sep;
+
+        output<<std::endl;
+        output.close();
+      } //end scope
+
+  }
+
+
+  template <int dim>
+  void
+  EigenvalueProblem<dim>::run()
+  {
+    for (unsigned int cycle = 0; cycle < 2; cycle++)
+      {
+        pcout << "Cycle "<<cycle <<std::endl;
+        const std::pair<unsigned int, unsigned int> n_cells = setup_system ();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_active_cells ()
+              << std::endl
+              << "     FE / POUFE :                "
+              << n_cells.first << " / " << n_cells.second
+              << std::endl
+              << "   Number of degrees of freedom: "
+              << dof_handler.n_dofs ()
+              << std::endl;
+
+        assemble_system ();
+
+        const std::pair<unsigned int, double> res = solve ();
+        AssertThrow(res.second < 5e-8,
+                    ExcInternalError());
+
+        estimate_error ();
+        // output_results(cycle);
+        refine_grid ();
+
+        pcout << std::endl;
+        for (unsigned int i=0; i<eigenvalues.size(); ++i)
+          pcout << "      Eigenvalue " << i
+                << " : " << eigenvalues[i]
+                << std::endl;
+      }
+
+  }
+}
+
+int main (int argc,char **argv)
+{
+
+  try
+    {
+      dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        Step36::EigenvalueProblem<dim> step36;
+        dealii::PETScWrappers::set_option_value("-eps_target","-1.0");
+        dealii::PETScWrappers::set_option_value("-st_type","sinvert");
+        dealii::PETScWrappers::set_option_value("-st_ksp_type","cg");
+        dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg");
+        dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11");
+        step36.run();
+      }
+
+    }
+  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;
+    };
+}
diff --git a/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..516893b
--- /dev/null
@@ -0,0 +1,12 @@
+Cycle 0
+   Number of active cells:       120
+     FE / POUFE :                120 / 8
+   Number of degrees of freedom: 1432
+
+      Eigenvalue 0 : -0.483398
+Cycle 1
+   Number of active cells:       176
+     FE / POUFE :                120 / 72
+   Number of degrees of freedom: 2206
+
+      Eigenvalue 0 : -0.498786

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.