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;
+}
<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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
+
DEAL_II_NAMESPACE_OPEN
+template <int dim, int spacedim> class FE_Enriched;
+
/**
* This class provides an interface to group several elements together into
unsigned int> >
base_elements;
-
/**
* Initialize the @p unit_support_points field of the FiniteElement class.
* Called from the constructor.
* Mutex for protecting initialization of restriction and embedding matrix.
*/
mutable Threads::Mutex mutex;
+
+ friend class FE_Enriched<dim,spacedim>;
};
fe_raviart_thomas_nodal.cc
fe_series.cc
fe_system.cc
+ fe_enriched.cc
fe_tools.cc
fe_tools_interpolate.cc
fe_trace.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::Test basic functions:
+DEAL::FE_Enriched<3>[FE_Q<3>(1)-FE_Q<3>(1)]
+DEAL::2
+DEAL::1
+DEAL::16
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::2cells: 2 1 2
+ 2 = 0
+ 5 = 0
+ 7 11: 0.500000
+ 7 13: 0.500000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::2cells: 2 1 2
+ 2 = 0
+ 5 = 0
+ 7 11: 0.500000
+ 7 13: 0.500000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
+DEAL::Ok
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+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