From: Denis Davydov Date: Wed, 28 Sep 2016 09:52:07 +0000 (+0200) Subject: add FE_Enriched finite element X-Git-Tag: v8.5.0-rc1~625^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3fc80271ca27f808368c276bff5e892526c469be;p=dealii.git add FE_Enriched finite element --- diff --git a/doc/doxygen/images/fe_enriched_1d.png b/doc/doxygen/images/fe_enriched_1d.png new file mode 100644 index 0000000000..6d1af3172d Binary files /dev/null and b/doc/doxygen/images/fe_enriched_1d.png differ diff --git a/doc/doxygen/images/fe_enriched_h-refinement.png b/doc/doxygen/images/fe_enriched_h-refinement.png new file mode 100644 index 0000000000..909cee39fe Binary files /dev/null and b/doc/doxygen/images/fe_enriched_h-refinement.png differ diff --git a/doc/doxygen/stylesheet.css b/doc/doxygen/stylesheet.css index 1562747f8c..36c1fa4fed 100644 --- a/doc/doxygen/stylesheet.css +++ b/doc/doxygen/stylesheet.css @@ -26,3 +26,11 @@ div.doxygen-generated-exception-message { padding: 5px; font-size: 95%; } + +div.image img[src="fe_enriched_1d.png"] { + width:500px; +} + +div.image img[src="fe_enriched_h-refinement.png"] { + width:500px; +} diff --git a/doc/news/changes.h b/doc/news/changes.h index de2cc39093..9859b2b5d2 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -206,6 +206,12 @@ inconvenience this causes.

General

    +
  1. 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. +
    + (Denis Davydov, 2016/09/28) +
  2. Improved: The doxygen documentation now contains nicely formatted boxes containing the text message of each exception. Several messages diff --git a/include/deal.II/fe/fe_enriched.h b/include/deal.II/fe/fe_enriched.h new file mode 100644 index 0000000000..b3131e46da --- /dev/null +++ b/include/deal.II/fe/fe_enriched.h @@ -0,0 +1,613 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__fe_enriched_h +#define dealii__fe_enriched_h + +#include + +#ifdef DEAL_II_WITH_CXX11 +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +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 fe(FE_Q(2), + * FE_Q(1), + * function) + * @endcode + * + * In this case, standard DoFs are distributed by FE_Q(2), + * whereas enriched DoFs are coming from a single finite element + * FE_Q(1) used with a single enrichment function + * function. In this case, the total number of DoFs on the + * enriched element is the sum of DoFs from FE_Q(2) and + * FE_Q(1). + * + * As an example of an enrichment function, consider $ \exp(-x) $, which + * leads to the following shape functions on the unit element: + * + * + * + * + * + * + * + * + * + *
    + * @image html fe_enriched_1d.png + * + * @image html fe_enriched_h-refinement.png + *
    + * 1d element, base and enriched shape functions. + * + * enriched shape function corresponding to the central vertex. + *
    + * + * 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} + * + *

    Using enriched and non-enriched FEs together

    + * + * 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 without actually enriching it. + * This can be done as follows: + * @code + * FE_Enriched fe_non_enriched(FE_Q(1)); + * @endcode + * This constructor is equivalent to calling + * @code + * FE_Enriched fe_non_enriched(FE_Q(1), + * FE_Nothing(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. + * + *

    References

    + * + * 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 + * + *

    Implementation

    + * + * The implementation of the class is based on FESystem which is aggregated as + * a private member. The simplest constructor FE_Enriched fe(FE_Q(2), FE_Q(1),function) + * will internally initialize FESystem as + * + * @code + * FESystem fe_system(FE_Q(2),1, + * FE_Q(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 +class FE_Enriched : public FiniteElement +{ +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 &fe_base, + const FiniteElement &fe_enriched, + const Function *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 &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 *fe_base, + const std::vector * > &fe_enriched, + const std::vector *(const typename Triangulation::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 + * multiplicities[0]=1, 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 *(const typename Triangulation::cell_iterator &) > > > &functions); +public: + + /** + * @p clone function instead of a copy constructor. + * + * This function is needed by the constructors of @p FESystem. + */ + virtual FiniteElement *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 & + 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 + * source.dofs_per_face times this->dofs_per_face. + * + * 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::ExcInterpolationNotImplemented, which + * will get propagated out from this element. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &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 + * source.dofs_per_face times this->dofs_per_face. + * + * 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::ExcInterpolationNotImplemented, which + * will get propagated out from this element. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &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 > + hp_vertex_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as hp_vertex_dof_indices(), except that the function treats degrees + * of freedom on lines. + */ + virtual + std::vector > + hp_line_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as hp_vertex_dof_indices(), except that the function treats degrees + * of freedom on quads. + */ + virtual + std::vector > + hp_quad_dof_identities (const FiniteElement &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 &fe_other) const; + //@} + + + /** + * Return enrichment functions + */ + const std::vector *(const typename Triangulation::cell_iterator &) > > > + get_enrichments() const; + + /** + * Return the underlying FESystem object. + */ + const FESystem & + get_fe_system() const; + +protected: + + /** + * A class to hold internal data needed for evaluation of this FE at quadrature points. + */ + class InternalData : public FiniteElement::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 values; + std::vector > gradients; + std::vector > 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::InternalData> fesystem_data); + + /** + * Gives read-access to the pointer to a @p InternalData of the @p + * base_noth base element of FESystem's data. + */ + typename FiniteElement::InternalDataBase & + get_fe_data (const unsigned int base_no) const; + + /** + * Gives read-access to the pointer to an object into which the + * base_noth base element will write its output when calling + * FiniteElement::fill_fe_values() and similar functions. + */ + internal::FEValues::FiniteElementRelatedData & + 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::InternalData> fesystem_data; + + /** + * For each FE used in enrichment (base number i) and each enrichment function + * (base multiplicity j), enrichment_values[i][j] will be used to store + * possibly requested values, gradients and hessians of enrichment function j. + * + * 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 > 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 > > 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 *(const typename Triangulation::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 false, otherwise it is + * true. + */ + 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 + typename FiniteElement::InternalDataBase * + setup_data (std::unique_ptr::InternalDataBase> fes_data, + const UpdateFlags flags, + const Quadrature &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::InternalDataBase * + get_data (const UpdateFlags flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const; + + virtual typename FiniteElement::InternalDataBase * + get_face_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const; + + virtual typename FiniteElement::InternalDataBase * + get_subface_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData &output_data) const; + + virtual + void fill_fe_values (const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &output_data + ) const; + + virtual + void + fill_fe_face_values ( const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &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 &quadrature, + const Mapping &mapping, + const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &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*> &fes, + const std::vector &multiplicities); + + /** + * The underlying FESystem object. + */ + FESystem 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 + void + multiply_by_enrichment (const Quadrature &quadrature, + const InternalData &fe_data, + const internal::FEValues::MappingRelatedData &mapping_data, + const typename Triangulation< dim, spacedim >::cell_iterator &cell, + internal::FEValues::FiniteElementRelatedData &output_data) const; +}; + +//} +DEAL_II_NAMESPACE_CLOSE + +#endif // CXX11 + +#endif // dealii__fe_enriched_h + diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 1f87be95e3..89b017bd1e 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -28,6 +28,8 @@ DEAL_II_NAMESPACE_OPEN +template class FE_Enriched; + /** * This class provides an interface to group several elements together into @@ -959,7 +961,6 @@ private: unsigned int> > base_elements; - /** * Initialize the @p unit_support_points field of the FiniteElement class. * Called from the constructor. @@ -1071,6 +1072,8 @@ private: * Mutex for protecting initialization of restriction and embedding matrix. */ mutable Threads::Mutex mutex; + + friend class FE_Enriched; }; diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 0e8f549797..f0a674c08c 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -44,6 +44,7 @@ SET(_src fe_raviart_thomas_nodal.cc fe_series.cc fe_system.cc + fe_enriched.cc fe_tools.cc fe_tools_interpolate.cc fe_trace.cc @@ -86,6 +87,7 @@ SET(_inst fe_raviart_thomas.inst.in fe_raviart_thomas_nodal.inst.in fe_system.inst.in + fe_enriched.inst.in fe_tools.inst.in fe_tools_interpolate.inst.in fe_trace.inst.in diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc new file mode 100644 index 0000000000..7486d6bedd --- /dev/null +++ b/source/fe/fe_enriched.cc @@ -0,0 +1,926 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +#include +#ifdef DEAL_II_WITH_CXX11 + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace +{ + /** + * Auxiliary function to create multiplicity vector from input enrichment functions. + */ + template + std::vector + build_multiplicities(const std::vector > &functions ) + { + std::vector 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 + std::vector< const FiniteElement< dim, spacedim > * > + build_fes(const FiniteElement *fe_base, + const std::vector * > &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 + bool + consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename Triangulation::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 *fe_nothing = dynamic_cast*>(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 + 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*>(fes[fe]) == 0) + // this is not FE_Nothing => there will be enrichment + return true; + + return false; + } +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement &fe_base) + : + FE_Enriched(fe_base, + FE_Nothing(fe_base.n_components(),true), + NULL) +{ +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement &fe_base, + const FiniteElement &fe_enriched, + const Function *enrichment_function) + : + FE_Enriched + (&fe_base, { &fe_enriched }, +{ + { + [=] (const typename Triangulation::cell_iterator &) -> const Function * + { + return enrichment_function; + } + } +}) +{ +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement *fe_base, + const std::vector * > &fe_enriched, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) + : + FE_Enriched (build_fes(fe_base,fe_enriched), + build_multiplicities(functions), + functions) +{ +} + + +template +FE_Enriched::FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) + : + FiniteElement (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_indexdofs_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 +const std::vector *(const typename Triangulation::cell_iterator &) > > > +FE_Enriched::get_enrichments() const +{ + return enrichments; +} + + +template +double +FE_Enriched::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 +FiniteElement * +FE_Enriched::clone() const +{ + std::vector< const FiniteElement< dim, spacedim > * > fes; + std::vector< unsigned int > multiplicities; + + for (unsigned int i=0; in_base_elements(); i++) + { + fes.push_back( & base_element(i) ); + multiplicities.push_back(this->element_multiplicity(i) ); + } + + return new FE_Enriched(fes, multiplicities, get_enrichments()); +} + + +template +UpdateFlags +FE_Enriched::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 +template +typename FiniteElement::InternalDataBase * +FE_Enriched::setup_data (std::unique_ptr::InternalDataBase> fes_data, + const UpdateFlags flags, + const Quadrature &quadrature) const +{ + Assert ((dynamic_cast::InternalData *> (fes_data.get()) != NULL), + ExcInternalError()); + typename FESystem::InternalData *data_fesystem = + static_cast::InternalData *> (fes_data.get()); + + // FESystem::InternalData will be aggregated (owned) by + // our InternalData. + fes_data.release(); + InternalData *data = new InternalData(std::unique_ptr::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 +typename FiniteElement::InternalDataBase * +FE_Enriched::get_face_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_face_data(update_flags,mapping,quadrature,output_data)), + update_flags, + quadrature); +} + + +template +typename FiniteElement::InternalDataBase * +FE_Enriched::get_subface_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_subface_data(update_flags,mapping,quadrature,output_data)), + update_flags, + quadrature); +} + + +template +typename FiniteElement::InternalDataBase * +FE_Enriched::get_data (const UpdateFlags flags, + const Mapping &mapping, + const Quadrature &quadrature, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_data(flags,mapping,quadrature,output_data)), + flags, + quadrature); +} + + +template +void FE_Enriched::initialize (const std::vector*> &fes, + const std::vector &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; i0) + 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 +std::string +FE_Enriched::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 +const FiniteElement & +FE_Enriched::base_element (const unsigned int index) const +{ + return fe_system.base_element(index); +} + + +template +void +FE_Enriched::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::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data + ) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (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 +void +FE_Enriched::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::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data +) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (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 +void +FE_Enriched::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::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data +) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (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 +template +void +FE_Enriched::multiply_by_enrichment +(const Quadrature &quadrature, + const InternalData &fe_data, + const internal::FEValues::MappingRelatedData &mapping_data, + const typename Triangulation< dim, spacedim >::cell_iterator &cell, + internal::FEValues::FiniteElementRelatedData &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 & + base_fe = base_element(base_no); + typename FiniteElement::InternalDataBase & + base_fe_data = fe_data.get_fe_data(base_no); + internal::FEValues::FiniteElementRelatedData & + 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_indexdofs_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_indexn_nonzero_components(i); + unsigned int in_index = 0; + for (unsigned int i=0; in_nonzero_components(system_index) == + base_fe.n_nonzero_components(base_index), + ExcInternalError()); + for (unsigned int s=0; sn_nonzero_components(system_index); ++s) + { + if (base_flags & update_values) + for (unsigned int q=0; qn_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; qhessian (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; qgradient(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; qvalue (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 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; qn_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 +const FESystem & +FE_Enriched:: +get_fe_system() const +{ + return fe_system; +} + + +template +bool +FE_Enriched:: +hp_constraints_are_implemented () const +{ + return true; +} + + +template +void +FE_Enriched:: +get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&source)) + { + fe_system.get_face_interpolation_matrix(fe_enr_other->get_fe_system(), + matrix); + } + else + { + typedef FiniteElement FEL; + AssertThrow(false,typename FEL::ExcInterpolationNotImplemented()); + } +} + + +template +void +FE_Enriched:: +get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&source)) + { + fe_system.get_subface_interpolation_matrix(fe_enr_other->get_fe_system(), + subface, + matrix); + } + else + { + typedef FiniteElement FEL; + AssertThrow(false,typename FEL::ExcInterpolationNotImplemented()); + } +} + + +template +std::vector > +FE_Enriched:: +hp_vertex_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_vertex_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +std::vector > +FE_Enriched:: +hp_line_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_line_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +std::vector > +FE_Enriched:: +hp_quad_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_quad_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +FiniteElementDomination::Domination +FE_Enriched:: +compare_for_face_domination (const FiniteElement &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 *fe_enr_other + = dynamic_cast*>(&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 +FE_Enriched::InternalData::InternalData(std::unique_ptr::InternalData> fesystem_data) + : + fesystem_data(std::move(fesystem_data)) +{} + + +template +typename FiniteElement::InternalDataBase & +FE_Enriched:: +InternalData::get_fe_data (const unsigned int base_no) const +{ + return fesystem_data->get_fe_data(base_no); +} + + +template +internal::FEValues::FiniteElementRelatedData & +FE_Enriched:: +InternalData::get_fe_output_object (const unsigned int base_no) const +{ + return fesystem_data->get_fe_output_object(base_no); +} + + +// explicit instantiations +#include "fe_enriched.inst" + +DEAL_II_NAMESPACE_CLOSE + +#endif // CXX11 diff --git a/source/fe/fe_enriched.inst.in b/source/fe/fe_enriched.inst.in new file mode 100644 index 0000000000..853f36a622 --- /dev/null +++ b/source/fe/fe_enriched.inst.in @@ -0,0 +1,23 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension + template class FE_Enriched; +#endif +} + diff --git a/tests/fe/fe_enriched_01.cc b/tests/fe/fe_enriched_01.cc new file mode 100644 index 0000000000..9413d6fe4a --- /dev/null +++ b/tests/fe/fe_enriched_01.cc @@ -0,0 +1,127 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test basic output of the class (name, n_blocks, n_compoennts, n_dofs_per_cell) + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test_base() +{ + deallog << "Test basic functions:"< function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + deallog << fe.get_name()<(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_01.with_cxx11=on.output b/tests/fe/fe_enriched_01.with_cxx11=on.output new file mode 100644 index 0000000000..d5e7fca629 --- /dev/null +++ b/tests/fe/fe_enriched_01.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::Test basic functions: +DEAL::FE_Enriched<3>[FE_Q<3>(1)-FE_Q<3>(1)] +DEAL::2 +DEAL::1 +DEAL::16 diff --git a/tests/fe/fe_enriched_02.cc b/tests/fe/fe_enriched_02.cc new file mode 100644 index 0000000000..fb6ac2069e --- /dev/null +++ b/tests/fe/fe_enriched_02.cc @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check FEValues::shape_value() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test1() +{ + deallog << "FEValues.shape_value()"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(1); + FEValues fe_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + typename DoFHandler::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 > &q_points = fe_values.get_quadrature_points(); + + for (unsigned int i=0; i(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_02.with_cxx11=on.output b/tests/fe/fe_enriched_02.with_cxx11=on.output new file mode 100644 index 0000000000..00053ee909 --- /dev/null +++ b/tests/fe/fe_enriched_02.with_cxx11=on.output @@ -0,0 +1,20 @@ + +DEAL::FEValues.shape_value() +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL::dof=0 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=1 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=2 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=3 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=4 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=5 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=6 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=7 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=8 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=9 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=10 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=11 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=12 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=13 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=14 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=15 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL:: diff --git a/tests/fe/fe_enriched_03.cc b/tests/fe/fe_enriched_03.cc new file mode 100644 index 0000000000..73ec999994 --- /dev/null +++ b/tests/fe/fe_enriched_03.cc @@ -0,0 +1,157 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check FEFaceValues::shape_value() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test2() +{ + deallog << "FEFaceValues"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(1); + FEFaceValues fe_face_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + for (unsigned int face = 0; face < GeometryInfo::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 > &q_points = fe_face_values.get_quadrature_points(); + + for (unsigned int i=0; i(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_03.with_cxx11=on.output b/tests/fe/fe_enriched_03.with_cxx11=on.output new file mode 100644 index 0000000000..ef9efb93fa --- /dev/null +++ b/tests/fe/fe_enriched_03.with_cxx11=on.output @@ -0,0 +1,99 @@ + +DEAL::FEFaceValues +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL::dof=0 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=3 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=4 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=5 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=6 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=7 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=8 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=9 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=10 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=11 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=12 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=13 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=14 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=3 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=4 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=5 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=6 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=7 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=8 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=9 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=10 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=11 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=12 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=13 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=14 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=0 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=3 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=4 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=5 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=6 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=7 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=8 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=9 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=10 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=11 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=12 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=13 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=14 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=3 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=4 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=5 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=6 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=7 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=8 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=9 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=10 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=11 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=12 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=13 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=14 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=0 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=3 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=4 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=5 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=6 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=7 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=8 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=9 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=10 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=11 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=12 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=13 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=14 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=3 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=4 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=5 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=6 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=7 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=8 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=9 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=10 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=11 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=12 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=13 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=14 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 diff --git a/tests/fe/fe_enriched_04.cc b/tests/fe/fe_enriched_04.cc new file mode 100644 index 0000000000..2d3d921658 --- /dev/null +++ b/tests/fe/fe_enriched_04.cc @@ -0,0 +1,244 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test FEValues::get_function_values() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test3() +{ + deallog << "FEValues.get_function_values()"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(2); + FEValues fe_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + Vector 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 solution_values_fe(n_q_points), + solution_values_pou(n_q_points), + solution_values(n_q_points); + typename DoFHandler::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 > &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 data_out; + data_out.attach_dof_handler (dof_handler); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + { + data_out.add_data_vector (shape_functions[i], names[i]); + } + + data_out.build_patches(patches); + + std::string filename = "shape_functions_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test3<3>(); +#ifdef DATA_OUT_FE_ENRICHED + plot_shape_function<1>(); + plot_shape_function<2>(); + plot_shape_function<3>(); +#endif + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_04.with_cxx11=on.output b/tests/fe/fe_enriched_04.with_cxx11=on.output new file mode 100644 index 0000000000..785ca6ae19 --- /dev/null +++ b/tests/fe/fe_enriched_04.with_cxx11=on.output @@ -0,0 +1,11 @@ + +DEAL::FEValues.get_function_values() +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL:: qp=0.2113 0.2113 0.2113 f(qp)=0.6935 U_fe(qp)=0.9811 U_pou(qp)=0.6804 U_[fe+pou](qp)=1.6615 +DEAL:: qp=0.7887 0.2113 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.2113 0.7887 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.7887 0.7887 0.2113 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.2113 0.2113 0.7887 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.7887 0.2113 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.2113 0.7887 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.7887 0.7887 0.7887 f(qp)=0.2551 U_fe(qp)=0.0189 U_pou(qp)=0.0048 U_[fe+pou](qp)=0.0237 diff --git a/tests/fe/fe_enriched_05.cc b/tests/fe/fe_enriched_05.cc new file mode 100644 index 0000000000..6a2a63ce30 --- /dev/null +++ b/tests/fe/fe_enriched_05.cc @@ -0,0 +1,186 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test h-refinement with DoFHandler and print constraints. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test4() +{ + deallog << "h-refinement:"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(2), + FE_Q(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> shape_functions; + std::vector 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 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 data_out; + data_out.attach_dof_handler (dof_handler); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector (shape_functions[i], names[i]); + + data_out.build_patches(patches); + + std::string filename = "h-shape_functions_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test4<2>(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_05.with_cxx11=on.output b/tests/fe/fe_enriched_05.with_cxx11=on.output new file mode 100644 index 0000000000..764434aac4 --- /dev/null +++ b/tests/fe/fe_enriched_05.with_cxx11=on.output @@ -0,0 +1,20 @@ + +DEAL::h-refinement: + 42 8: 1.00000 + 43 1: 0.500000 + 43 5: 0.500000 + 44 0: 0.375000 + 44 4: -0.125000 + 44 8: 0.750000 + 48 21: 1.00000 + 49 5: 0.500000 + 49 14: 0.500000 + 52 4: -0.125000 + 52 13: 0.375000 + 52 21: 0.750000 + 54 0: -0.125000 + 54 4: 0.375000 + 54 8: 0.750000 + 55 4: 0.375000 + 55 13: -0.125000 + 55 21: 0.750000 diff --git a/tests/fe/fe_enriched_06.cc b/tests/fe/fe_enriched_06.cc new file mode 100644 index 0000000000..01b33e104e --- /dev/null +++ b/tests/fe/fe_enriched_06.cc @@ -0,0 +1,223 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// no refinement, just 4 cells. +// Make sure fe_values and shape functions are continuous + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test5() +{ + deallog << "4 cells:"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + + triangulation.refine_global(); + + dof_handler.distribute_dofs(fe); + + std::vector> shape_functions; + std::vector 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 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 &solution = shape_functions[global_dof]; + QTrapez quadrature; + FEValues fe_values(fe, quadrature, + update_values); + + const unsigned int n_q_points = quadrature.size(); + std::vector solution_values(n_q_points); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + std::vector local_dof_indices; + local_dof_indices.resize (dofs_per_cell); + + typename DoFHandler::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 > &q_points = fe_values.get_quadrature_points(); + fe_values.get_function_values (solution, + solution_values); + + deallog <<" cell="<center()<(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_06.with_cxx11=on.output b/tests/fe/fe_enriched_06.with_cxx11=on.output new file mode 100644 index 0000000000..a45dcea52b --- /dev/null +++ b/tests/fe/fe_enriched_06.with_cxx11=on.output @@ -0,0 +1,22 @@ + +DEAL::4 cells: +DEAL:: cell=0.2500 0.2500 +DEAL:: qp=0.0000 0.0000 f(qp)=1.0000 U(qp)=0 +DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0 +DEAL:: cell=0.7500 0.2500 +DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 N(7,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.0000 f(qp)=0.3679 N(7,qp)=0 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(7,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(7,qp)=0.3269 U(qp)=0.3269 +DEAL:: cell=0.2500 0.7500 +DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0 +DEAL:: qp=0.0000 1.0000 f(qp)=0.3679 U(qp)=0 +DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 U(qp)=0 +DEAL:: cell=0.7500 0.7500 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(3,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(3,qp)=0.3269 U(qp)=0.3269 +DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 N(3,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 1.0000 f(qp)=0.2431 N(3,qp)=0 U(qp)=0 diff --git a/tests/fe/fe_enriched_07.cc b/tests/fe/fe_enriched_07.cc new file mode 100644 index 0000000000..d0512022fe --- /dev/null +++ b/tests/fe/fe_enriched_07.cc @@ -0,0 +1,237 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// hp-refinement. Output constraints. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test6(const bool do_href, + const unsigned int p_feq = 1, + const unsigned int p_feen = 2) +{ + deallog << "hp: "< triangulation; + hp::DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + + hp::FECollection fe_collection; + fe_collection.push_back(FE_Enriched(FE_Q(p_feq))); + fe_collection.push_back(FE_Enriched(FE_Q(p_feen), + FE_Q(1), + &function)); + // push back to be able to resolve hp constrains: + fe_collection.push_back(FE_Enriched(FE_Q(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::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: "<> shape_functions; + std::vector names; + for (unsigned int s=0; s < dof_handler.n_dofs(); s++) + { + Vector 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> data_out; + data_out.attach_dof_handler (dof_handler); + + // get material ids: + Vector fe_index(triangulation.n_active_cells()); + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (unsigned int index=0; cell!=endc; ++cell,++index) + { + fe_index[index] = cell->active_fe_index(); + } + data_out.add_data_vector(fe_index, "fe_index"); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector (shape_functions[i], names[i]); + + data_out.build_patches(patches); + + std::string filename = "hp-shape_functions_ref="+std::to_string(do_href) + +"_p_feq="+std::to_string(p_feq) + +"_p_feenr="+std::to_string(p_feen) + +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test6<2>(false,1,2);// 1 vs 2+1 + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_07.with_cxx11=on.output b/tests/fe/fe_enriched_07.with_cxx11=on.output new file mode 100644 index 0000000000..e4f68ec21f --- /dev/null +++ b/tests/fe/fe_enriched_07.with_cxx11=on.output @@ -0,0 +1,10 @@ + +DEAL::hp: 0 1 2 +DEAL::n_cells: 4 + 4 = 0 + 5 = 0 + 9 18: 0.500000 + 9 19: 0.500000 + 13 = 0 + 16 19: 0.500000 + 16 22: 0.500000 diff --git a/tests/fe/fe_enriched_08.cc b/tests/fe/fe_enriched_08.cc new file mode 100644 index 0000000000..c59363d6bc --- /dev/null +++ b/tests/fe/fe_enriched_08.cc @@ -0,0 +1,219 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check constraints on hp dofhandler with 2 neighboring cells where one +// cell is enriched and nother is not. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 +void test2cells(const unsigned int p_feq=2, + const unsigned int p_feen=1) +{ + deallog << "2cells: "< triangulation; + { + Triangulation triangulationL; + Triangulation 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 shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + } + + hp::DoFHandler dof_handler(triangulation); + EnrichmentFunction function; + + hp::FECollection fe_collection; + fe_collection.push_back(FE_Enriched(FE_Q(p_feq))); + fe_collection.push_back(FE_Enriched(FE_Q(p_feen), + FE_Q(1), + &function)); + // push back to be able to resolve hp constrains: + fe_collection.push_back(FE_Enriched(FE_Q(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> shape_functions; + std::vector names; + for (unsigned int s=0; s < dof_handler.n_dofs(); s++) + { + Vector 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> data_out; + data_out.attach_dof_handler (dof_handler); + + // get material ids: + Vector fe_index(triangulation.n_active_cells()); + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (unsigned int index=0; cell!=endc; ++cell,++index) + { + fe_index[index] = cell->active_fe_index(); + } + data_out.add_data_vector(fe_index, "fe_index"); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector (shape_functions[i], names[i]); + + data_out.build_patches(patches); + + std::string filename = "2cell-shape_functions_p_feq="+std::to_string(p_feq) + +"_p_feenr="+std::to_string(p_feen) + +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test2cells<2>(1,2); // 1 vs 2+1 + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_08.with_cxx11=on.output b/tests/fe/fe_enriched_08.with_cxx11=on.output new file mode 100644 index 0000000000..3a075cc441 --- /dev/null +++ b/tests/fe/fe_enriched_08.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::2cells: 2 1 2 + 2 = 0 + 5 = 0 + 7 11: 0.500000 + 7 13: 0.500000 diff --git a/tests/fe/fe_enriched_08a.cc b/tests/fe/fe_enriched_08a.cc new file mode 100644 index 0000000000..23fd6dbdae --- /dev/null +++ b/tests/fe/fe_enriched_08a.cc @@ -0,0 +1,133 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Same as fe_enriched_08, but checks the constraints from FESystem directly. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +using namespace dealii; + +template +void test2cellsFESystem(const unsigned int p_feq=2, + const unsigned int p_feen=1) +{ + deallog << "2cells: "< triangulation; + { + Triangulation triangulationL; + Triangulation 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 shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + } + + hp::DoFHandler dof_handler(triangulation); + + hp::FECollection fe_collection; + fe_collection.push_back(FESystem(FE_Q(p_feq),1, + FE_Nothing(),1)); + fe_collection.push_back(FESystem(FE_Q(p_feen),1, + FE_Q(1),1)); + + // push back to be able to resolve hp constrains: + fe_collection.push_back(FESystem(FE_Q(p_feen),1, + FE_Nothing(),1)); + + dof_handler.begin_active()->set_active_fe_index(1); + + dof_handler.distribute_dofs(fe_collection); + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + constraints.print(deallog.get_file_stream()); + + dof_handler.clear(); +} + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test2cellsFESystem<2>(1,2); // 1+0 vs 2+1 + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_08a.with_cxx11=on.output b/tests/fe/fe_enriched_08a.with_cxx11=on.output new file mode 100644 index 0000000000..3a075cc441 --- /dev/null +++ b/tests/fe/fe_enriched_08a.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::2cells: 2 1 2 + 2 = 0 + 5 = 0 + 7 11: 0.500000 + 7 13: 0.500000 diff --git a/tests/fe/fe_enriched_compare_to_fe_system.cc b/tests/fe/fe_enriched_compare_to_fe_system.cc new file mode 100644 index 0000000000..f19cb7c120 --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system.cc @@ -0,0 +1,320 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test FE_Enriched class by comparing with FE_System for a single enrichment +// function. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &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 &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()/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 +void check_consistency(const Point &p, + const Function &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 +void test(const FiniteElement &fe1, + const FiniteElement &fe2, + const Quadrature &volume_quad, + const Quadrature &face_quad, + const bool distort) +{ + Triangulation triangulation; + DoFHandler dof_handler_enriched(triangulation); + DoFHandler dof_handler_system(triangulation); + + EnrichmentFunction function; + FE_Enriched fe_enriched(fe1, + fe2, + &function); + FESystem fe_system(fe1,1, + fe2,1); + + { + Point p1,p2; + std::vector 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 fe_values_enriched(fe_enriched, volume_quad, + update_values | update_gradients | update_hessians ); + + FEValues fe_values_system(fe_system, volume_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + FEFaceValues fe_face_values_enriched(fe_enriched, face_quad, + update_values | update_gradients | update_hessians); + + FEFaceValues 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::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 > &q_points = fe_values_system.get_quadrature_points(); + + // check shape functions + for (unsigned int i=0; i::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 > &q_points = fe_face_values_system.get_quadrature_points(); + + for (unsigned int i=0; i(3),FE_Q(2),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),true); + } + + { + const unsigned int dim = 3; + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),true); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output new file mode 100644 index 0000000000..f2b98e8f5b --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output @@ -0,0 +1,5 @@ + +DEAL::Ok +DEAL::Ok +DEAL::Ok +DEAL::Ok diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.cc b/tests/fe/fe_enriched_compare_to_fe_system_2.cc new file mode 100644 index 0000000000..ca0e3c8468 --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system_2.cc @@ -0,0 +1,403 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test FE_Enriched class by comparing with FE_System for multiple enrichment +// functions. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +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 +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin) + : Function(1), + origin(origin) + {} + + virtual double value(const Point &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 &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 &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()/r ); + } + +private: + Point 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 +void check_consistency(const Point &p, + const Function &func1, + const Function &func2, + const Function &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 +void test(const FiniteElement &fe_base, + const FiniteElement &fe_en1, + const FiniteElement &fe_en2, + const Quadrature &volume_quad, + const Quadrature &face_quad, + const bool distort) +{ + Triangulation triangulation; + DoFHandler dof_handler_enriched(triangulation); + DoFHandler dof_handler_system(triangulation); + + Point 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 fun1(p1), fun2(p2),fun3(p3); + std::vector *> fe_enrichements(2); + fe_enrichements[0] = &fe_en1; + fe_enrichements[1] = &fe_en2; + std::vector *(const typename Triangulation::cell_iterator &) > > > functions(2); + functions[0].resize(1); + functions[0][0] = [&](const typename Triangulation::cell_iterator &) + { + return &fun1; + }; + functions[1].resize(2); + functions[1][0] = [&](const typename Triangulation::cell_iterator &) + { + return &fun2; + }; + functions[1][1] = [&](const typename Triangulation::cell_iterator &) + { + return &fun3; + }; + + FE_Enriched fe_enriched(&fe_base, + fe_enrichements, + functions); + FESystem fe_system(fe_base,1, + fe_en1, 1, + fe_en2, 2); + + { + Point p1,p2; + std::vector 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 fe_values_enriched(fe_enriched, volume_quad, + update_values | update_gradients | update_hessians ); + + FEValues fe_values_system(fe_system, volume_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + FEFaceValues fe_face_values_enriched(fe_enriched, face_quad, + update_values | update_gradients | update_hessians); + + FEFaceValues 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 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::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 > &q_points = fe_values_system.get_quadrature_points(); + + // check shape functions + for (unsigned int i=0; i::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 > &q_points = fe_face_values_system.get_quadrature_points(); + + for (unsigned int i=0; i(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),true); + } + + { + const unsigned int dim = 3; + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),true); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output new file mode 100644 index 0000000000..f2b98e8f5b --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output @@ -0,0 +1,5 @@ + +DEAL::Ok +DEAL::Ok +DEAL::Ok +DEAL::Ok diff --git a/tests/fe/fe_enriched_step-36.cc b/tests/fe/fe_enriched_step-36.cc new file mode 100644 index 0000000000..a22183a967 --- /dev/null +++ b/tests/fe/fe_enriched_step-36.cc @@ -0,0 +1,779 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test FE_Enriched in real-life application on eigenvalue problem similar +// to Step-36. That involves assembly (shape values and gradients) and +// error estimator (Kelly - > face gradients) and MPI run. + +#include "../tests.h" + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +const unsigned int dim = 3; + +using namespace dealii; + +/** + * Coulomb potential + */ +template +class PotentialFunction : public dealii::Function +{ +public: + PotentialFunction() + : dealii::Function(1) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +double PotentialFunction::value(const dealii::Point &p, + const unsigned int ) const +{ + return -1.0 / std::sqrt(p.square()); +} + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin, + const double &Z, + const double &radius) + : Function(1), + origin(origin), + Z(Z), + radius(radius) + {} + + virtual double value(const Point &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 &point) const + { + if (origin.distance(point) < radius) + return true; + else + return false; + } + + virtual Tensor< 1, dim> gradient (const Point &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 origin; + + /** + * charge + */ + const double Z; + + /** + * enrichment radius + */ + const double radius; +}; + +namespace Step36 +{ + + /** + * Main class + */ + template + class EigenvalueProblem + { + public: + EigenvalueProblem (); + virtual ~EigenvalueProblem(); + void run (); + + private: + bool cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const; + + std::pair setup_system (); + void assemble_system (); + std::pair solve (); + void constrain_pou_dofs(); + void estimate_error (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection q_collection; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector eigenfunctions_locally_relevant; + std::vector 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 potential; + + EnrichmentFunction 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 > vec_estimated_error_per_cell; + Vector estimated_error_per_cell; + }; + + template + EigenvalueProblem::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(), + /*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(4)); + q_collection.push_back(QGauss(10)); + + // usual elements (active_fe_index ==0): + fe_collection.push_back (FE_Enriched (FE_Q(2))); + + + // enriched elements (active_fe_index==1): + fe_collection.push_back (FE_Enriched (FE_Q(2), + FE_Q(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::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 + std::pair + EigenvalueProblem::setup_system () + { + for (typename hp::DoFHandler::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 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 (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 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::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 + bool + EigenvalueProblem::cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const + { + return cell->material_id() == pou_material_id; + } + + template + void + EigenvalueProblem::constrain_pou_dofs() + { + std::vector local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face); + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + for (unsigned int f=0; f::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; sfface(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 0) + constraints.add_line (local_face_dof_indices[i]); + } + } + } + + template + void + EigenvalueProblem::assemble_system() + { + stiffness_matrix = 0; + mass_matrix = 0; + + dealii::FullMatrix cell_stiffness_matrix; + dealii::FullMatrix cell_mass_matrix; + std::vector local_dof_indices; + std::vector potential_values; + std::vector enrichment_values; + std::vector> enrichment_gradients; + + hp::FEValues fe_values_hp(fe_collection, q_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + + typename hp::DoFHandler::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 &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_pointget_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 + std::pair + EigenvalueProblem::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 + EigenvalueProblem::~EigenvalueProblem() + { + dof_handler.clear (); + } + + template + void EigenvalueProblem::estimate_error () + { + { + std::vector sol (number_of_eigenvalues); + std::vector *> 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 face_quadrature_formula; + face_quadrature_formula.push_back (dealii::QGauss(3)); + face_quadrature_formula.push_back (dealii::QGauss(3)); + + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_formula, + typename FunctionMap::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 + void EigenvalueProblem::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 + class Postprocessor : public DataPostprocessorScalar + { + public: + Postprocessor(const EnrichmentFunction &enrichment); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + private: + const EnrichmentFunction &enrichment; + }; + + template + Postprocessor::Postprocessor(const EnrichmentFunction &enrichment) + : + DataPostprocessorScalar ("total_solution", + update_values | update_q_points), + enrichment(enrichment) + {} + + template + void + Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &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 + void EigenvalueProblem::output_results (const unsigned int cycle) const + { + dealii::Vector fe_index(triangulation.n_active_cells()); + { + typename dealii::hp::DoFHandler::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 postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration + DataOut > 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 > 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::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< + void + EigenvalueProblem::run() + { + for (unsigned int cycle = 0; cycle < 5; cycle++) + { + pcout << "Cycle "< 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 res = solve (); + AssertThrow(res.second < 5e-8, + ExcInternalError()); + + estimate_error (); + // output_results(cycle); + refine_grid (); + + pcout << std::endl; + for (unsigned int i=0; i step36; + dealii::PETScWrappers::set_option_value("-eps_target","-1.0"); + dealii::PETScWrappers::set_option_value("-st_type","sinvert"); + dealii::PETScWrappers::set_option_value("-st_ksp_type","cg"); + dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg"); + dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); + step36.run(); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output new file mode 100644 index 0000000000..fa33870cc1 --- /dev/null +++ b/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output @@ -0,0 +1,30 @@ +Cycle 0 + Number of active cells: 120 + FE / POUFE : 120 / 8 + Number of degrees of freedom: 1432 + + Eigenvalue 0 : -0.483398 +Cycle 1 + Number of active cells: 176 + FE / POUFE : 120 / 72 + Number of degrees of freedom: 2206 + + Eigenvalue 0 : -0.498786 +Cycle 2 + Number of active cells: 512 + FE / POUFE : 504 / 72 + Number of degrees of freedom: 5734 + + Eigenvalue 0 : -0.499788 +Cycle 3 + Number of active cells: 568 + FE / POUFE : 504 / 136 + Number of degrees of freedom: 6508 + + Eigenvalue 0 : -0.499818 +Cycle 4 + Number of active cells: 736 + FE / POUFE : 504 / 328 + Number of degrees of freedom: 8542 + + Eigenvalue 0 : -0.499947 diff --git a/tests/fe/fe_enriched_step-36b.cc b/tests/fe/fe_enriched_step-36b.cc new file mode 100644 index 0000000000..4a92b81d87 --- /dev/null +++ b/tests/fe/fe_enriched_step-36b.cc @@ -0,0 +1,903 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test FE_Enriched in real-life application on eigenvalue problem similar +// to Step-36. That involves assembly (shape values and gradients) and +// error estimator (Kelly - > face gradients) and MPI run. +// +// same as fe_enriched_step-36, but using FESystem only. Note that the results +// of KellyErrorEstimator won't be the same, thus this tests does only a single +// refinement cycle where the refinement path is the same. + +#include "../tests.h" + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +const unsigned int dim = 3; + +using namespace dealii; + +/** + * Coulomb potential + */ +template +class PotentialFunction : public dealii::Function +{ +public: + PotentialFunction() + : dealii::Function(1) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +double PotentialFunction::value(const dealii::Point &p, + const unsigned int ) const +{ + return -1.0 / std::sqrt(p.square()); +} + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin, + const double &Z, + const double &radius) + : Function(1), + origin(origin), + Z(Z), + radius(radius) + {} + + virtual double value(const Point &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 &point) const + { + if (origin.distance(point) < radius) + return true; + else + return false; + } + + virtual Tensor< 1, dim> gradient (const Point &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 origin; + + /** + * charge + */ + const double Z; + + /** + * enrichment radius + */ + const double radius; +}; + +namespace Step36 +{ + + /** + * Main class + */ + template + class EigenvalueProblem + { + public: + EigenvalueProblem (); + virtual ~EigenvalueProblem(); + void run (); + + private: + bool cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const; + + std::pair setup_system (); + void assemble_system (); + std::pair solve (); + void constrain_pou_dofs(); + void estimate_error (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection q_collection; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector eigenfunctions_locally_relevant; + std::vector 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 potential; + + EnrichmentFunction 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 > vec_estimated_error_per_cell; + Vector estimated_error_per_cell; + }; + + template + EigenvalueProblem::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(), + /*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(4)); + q_collection.push_back(QGauss(10)); + + // usual elements (active_fe_index ==0): + fe_collection.push_back (FESystem (FE_Q(2), 1, + FE_Nothing(), 1)); + + // enriched elements (active_fe_index==1): + fe_collection.push_back (FESystem (FE_Q(2), 1, + FE_Q(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::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 + std::pair + EigenvalueProblem::setup_system () + { + for (typename hp::DoFHandler::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 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 (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 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::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 + bool + EigenvalueProblem::cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const + { + return cell->material_id() == pou_material_id; + } + + template + void + EigenvalueProblem::constrain_pou_dofs() + { + std::vector local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face); + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + for (unsigned int f=0; f::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; sfface(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 0) + constraints.add_line (local_face_dof_indices[i]); + } + } + } + + template + void + EigenvalueProblem::assemble_system() + { + stiffness_matrix = 0; + mass_matrix = 0; + + dealii::FullMatrix cell_stiffness_matrix; + dealii::FullMatrix cell_mass_matrix; + std::vector local_dof_indices; + std::vector potential_values; + std::vector enrichment_values; + std::vector> enrichment_gradients; + + hp::FEValues fe_values_hp(fe_collection, q_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + + typename hp::DoFHandler::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 &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_pointactive_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; iget_fe().system_to_base_index(i).first.first; + for (unsigned int j=i; jget_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_pointget_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 + std::pair + EigenvalueProblem::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 + EigenvalueProblem::~EigenvalueProblem() + { + dof_handler.clear (); + } + + template + void EigenvalueProblem::estimate_error () + { + { + std::vector sol (number_of_eigenvalues); + std::vector *> 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 face_quadrature_formula; + face_quadrature_formula.push_back (dealii::QGauss(3)); + face_quadrature_formula.push_back (dealii::QGauss(3)); + + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_formula, + typename FunctionMap::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 + void EigenvalueProblem::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 + class Postprocessor : public DataPostprocessorScalar + { + public: + Postprocessor(const EnrichmentFunction &enrichment); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + private: + const EnrichmentFunction &enrichment; + }; + + template + Postprocessor::Postprocessor(const EnrichmentFunction &enrichment) + : + DataPostprocessorScalar ("total_solution", + update_values | update_q_points), + enrichment(enrichment) + {} + + template + void + Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &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 + void EigenvalueProblem::output_results (const unsigned int cycle) const + { + dealii::Vector fe_index(triangulation.n_active_cells()); + { + typename dealii::hp::DoFHandler::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 postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration + DataOut > 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 > 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::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< + void + EigenvalueProblem::run() + { + for (unsigned int cycle = 0; cycle < 2; cycle++) + { + pcout << "Cycle "< 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 res = solve (); + AssertThrow(res.second < 5e-8, + ExcInternalError()); + + estimate_error (); + // output_results(cycle); + refine_grid (); + + pcout << std::endl; + for (unsigned int i=0; i step36; + dealii::PETScWrappers::set_option_value("-eps_target","-1.0"); + dealii::PETScWrappers::set_option_value("-st_type","sinvert"); + dealii::PETScWrappers::set_option_value("-st_ksp_type","cg"); + dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg"); + dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); + step36.run(); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output new file mode 100644 index 0000000000..516893b783 --- /dev/null +++ b/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output @@ -0,0 +1,12 @@ +Cycle 0 + Number of active cells: 120 + FE / POUFE : 120 / 8 + Number of degrees of freedom: 1432 + + Eigenvalue 0 : -0.483398 +Cycle 1 + Number of active cells: 176 + FE / POUFE : 120 / 72 + Number of degrees of freedom: 2206 + + Eigenvalue 0 : -0.498786