From: kronbichler Date: Tue, 11 Feb 2014 08:33:28 +0000 (+0000) Subject: Also create 1D version of FE_FaceQ/P to allow step-51 to run in 1D (in principle... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e4ae99d6aca9faabfff6f0adeb3d8581302b0f02;p=dealii-svn.git Also create 1D version of FE_FaceQ/P to allow step-51 to run in 1D (in principle, output of faces does not work there and is not that easy to fix). Add a test for step-51. Test data of FE_FaceQ/P<1>. Fix output of test due to fix in r32451. git-svn-id: https://svn.dealii.org/trunk@32454 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 6a14027983..e62fcd2e49 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,16 @@ inconvenience this causes.

Specific improvements

    +
  1. New: FE_FaceQ and FE_FaceP now also work in 1D (with a single dof + on each vertex). +
    + (Martin Kronbichler, 2014/02/11) + +
  2. Fixed: FE_DGQ::has_support_on_face returned a wrong number for element + degree larger than 1 in 1D. This is now fixed. +
    + (Martin Kronbichler, 2014/02/10) +
  3. Changed: DerivativeApproximation used to be a class that only had static members. It is now a namespace.
    diff --git a/deal.II/include/deal.II/fe/fe_face.h b/deal.II/include/deal.II/fe/fe_face.h index 6e80608548..0831942687 100644 --- a/deal.II/include/deal.II/fe/fe_face.h +++ b/deal.II/include/deal.II/fe/fe_face.h @@ -124,7 +124,6 @@ public: FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement &fe_other) const; - /** * Returns a list of constant modes of the element. For this element, it * simply returns one row with all entries set to true. @@ -140,6 +139,197 @@ private: +/** + * Specialization of FE_FaceQ for 1D. In that case, the finite element only + * consists of one degree of freedom in each of the two faces (= vertices) of + * a cell, irrespective of the degree. However, this element still accepts a + * degree in its constructor and also returns that degree. This way, + * dimension-independent programming with trace elements is also possible in + * 1D (even though there is no computational benefit at all from it in 1D). + * + * @ingroup fe + * @author Guido Kanschat, Martin Kronbichler + * @date 2014 + */ +template +class FE_FaceQ<1,spacedim> : public FiniteElement<1,spacedim> +{ +public: + /** + * Constructor. + */ + FE_FaceQ (const unsigned int p); + + /** + * Clone method. + */ + virtual FiniteElement<1,spacedim> *clone() const; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_FaceQ(degree), with dim and + * degree replaced by appropriate values. + */ + virtual std::string get_name () 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. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source, + FullMatrix &matrix) 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. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &source, + const unsigned int subface, + FullMatrix &matrix) const; + + /** + * Check for non-zero values on a face. + * + * This function returns @p true, if the shape function @p shape_index has + * non-zero values on the face @p face_index. + * + * Implementation of the interface in FiniteElement + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + /** + * Return whether this element implements its hanging node constraints in + * the new way, which has to be used to make elements "hp compatible". + */ + virtual bool hp_constraints_are_implemented () 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<1,spacedim> &fe_other) const; + + /** + * Returns a list of constant modes of the element. For this element, it + * simply returns one row with all entries set to true. + */ + virtual Table<2,bool> get_constant_modes () const; + +protected: + virtual + typename Mapping<1,spacedim>::InternalDataBase * + get_data (const UpdateFlags, + const Mapping<1,spacedim> &mapping, + const Quadrature<1> &quadrature) const ; + + typename Mapping<1,spacedim>::InternalDataBase * + get_face_data (const UpdateFlags, + const Mapping<1,spacedim> &mapping, + const Quadrature<0>& quadrature) const ; + + typename Mapping<1,spacedim>::InternalDataBase * + get_subface_data (const UpdateFlags, + const Mapping<1,spacedim> &mapping, + const Quadrature<0>& quadrature) const ; + + virtual void + fill_fe_values (const Mapping<1,spacedim> &mapping, + const typename Triangulation<1,spacedim>::cell_iterator &cell, + const Quadrature<1> &quadrature, + typename Mapping<1,spacedim>::InternalDataBase &mapping_internal, + typename Mapping<1,spacedim>::InternalDataBase &fe_internal, + FEValuesData<1,spacedim> &data, + CellSimilarity::Similarity &cell_similarity) const; + + virtual void + fill_fe_face_values (const Mapping<1,spacedim> &mapping, + const typename Triangulation<1,spacedim>::cell_iterator &cell, + const unsigned int face_no, + const Quadrature<0> &quadrature, + typename Mapping<1,spacedim>::InternalDataBase &mapping_internal, + typename Mapping<1,spacedim>::InternalDataBase &fe_internal, + FEValuesData<1,spacedim> &data) const ; + + virtual void + fill_fe_subface_values (const Mapping<1,spacedim> &mapping, + const typename Triangulation<1,spacedim>::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature<0> &quadrature, + typename Mapping<1,spacedim>::InternalDataBase &mapping_internal, + typename Mapping<1,spacedim>::InternalDataBase &fe_internal, + FEValuesData<1,spacedim> &data) const ; + + + /** + * Determine the values that need to be computed on the unit cell to be able + * to compute all values required by flags. + * + * For the purpuse of this function, refer to the documentation in + * FiniteElement. + * + * This class assumes that shape functions of this FiniteElement do + * not depend on the actual shape of the cells in real + * space. Therefore, the effect in this element is as follows: if + * update_values is set in flags, copy it to the + * result. All other flags of the result are cleared, since everything else + * must be computed for each cell. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + + /** + * Determine the values that need to be computed on every cell to be able to + * compute all values required by flags. + * + * For the purpuse of this function, refer to the documentation in + * FiniteElement. + * + * This class assumes that shape functions of this FiniteElement do + * not depend on the actual shape of the cells in real space. + * + * The effect in this element is as follows: + *
      + *
    • if update_gradients is set, the result will contain + * update_gradients and update_covariant_transformation. + * The latter is required to transform the gradient on the unit cell to the + * real cell. Remark, that the action required by + * update_covariant_transformation is actually performed by the + * Mapping object used in conjunction with this finite element.
    • if + * update_hessians is set, the result will contain + * update_hessians and update_covariant_transformation. + * The rationale is the same as above and no higher derivatives of the + * transformation are required, since we use difference quotients for the + * actual computation. + * + *
    + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; + +private: + /** + * Return vector with dofs per vertex, line, quad, hex. + */ + static std::vector get_dpo_vector (const unsigned int deg); +}; + + /** * A finite element, which is a Legendre element of complete polynomials on @@ -257,6 +447,26 @@ private: }; + +/** + * FE_FaceP in 1D, i.e., with degrees of freedom on the element vertices. + */ +template +class FE_FaceP<1,spacedim> : public FE_FaceQ<1,spacedim> +{ +public: + /** + * Constructor. + */ + FE_FaceP (const unsigned int p); + + /** + * Returns the name of the element + */ + std::string get_name() const; +}; + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/include/deal.II/fe/fe_poly_face.templates.h b/deal.II/include/deal.II/fe/fe_poly_face.templates.h index 39651c0170..c1067d7521 100644 --- a/deal.II/include/deal.II/fe/fe_poly_face.templates.h +++ b/deal.II/include/deal.II/fe/fe_poly_face.templates.h @@ -56,11 +56,8 @@ template UpdateFlags FE_PolyFace::update_once (const UpdateFlags) const { - // for this kind of elements, only - // the values can be precomputed - // once and for all. set this flag - // if the values are requested at - // all + // for this kind of elements, only the values can be precomputed once and + // for all. set this flag if the values are requested at all return update_default; } diff --git a/deal.II/source/fe/fe_face.cc b/deal.II/source/fe/fe_face.cc index 728a6d2498..9d62ff19d5 100644 --- a/deal.II/source/fe/fe_face.cc +++ b/deal.II/source/fe/fe_face.cc @@ -300,6 +300,270 @@ FE_FaceQ::get_constant_modes () const +// ----------------------------- FE_FaceQ<1,spacedim> ------------------------ + +template +FE_FaceQ<1,spacedim>::FE_FaceQ (const unsigned int degree) + : + FiniteElement<1,spacedim> (FiniteElementData<1>(get_dpo_vector(degree), 1, degree, FiniteElementData<1>::L2), + std::vector(1,true), + std::vector (1, ComponentMask(1,true))) +{ + this->unit_face_support_points.resize(1); + + // initialize unit support points (this makes it possible to assign initial + // values to FE_FaceQ) + this->unit_support_points.resize(GeometryInfo<1>::faces_per_cell); + this->unit_support_points[1] = Point<1>(1.); +} + + + +template +FiniteElement<1,spacedim> * +FE_FaceQ<1,spacedim>::clone() const +{ + return new FE_FaceQ<1,spacedim>(this->degree); +} + + + +template +std::string +FE_FaceQ<1,spacedim>::get_name () const +{ + // note that the FETools::get_fe_from_name function depends on the + // particular format of the string this function returns, so they have to be + // kept in synch + std::ostringstream namebuf; + namebuf << "FE_FaceQ<1>(" << this->degree << ")"; + + return namebuf.str(); +} + + + +template +void +FE_FaceQ<1,spacedim>:: +get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source_fe, + FullMatrix &interpolation_matrix) const +{ + get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int, + interpolation_matrix); +} + + + +template +void +FE_FaceQ<1,spacedim>:: +get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &x_source_fe, + const unsigned int subface, + FullMatrix &interpolation_matrix) const +{ + Assert (interpolation_matrix.n() == this->dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.n(), + this->dofs_per_face)); + Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_face)); + interpolation_matrix(0,0) = 1.; +} + + + +template +bool +FE_FaceQ<1,spacedim>::has_support_on_face ( + const unsigned int shape_index, + const unsigned int face_index) const +{ + AssertIndexRange(shape_index, 2); + return (face_index == shape_index); +} + + + +template +std::vector +FE_FaceQ<1,spacedim>::get_dpo_vector (const unsigned int) +{ + std::vector dpo(2, 0U); + dpo[0] = 1; + return dpo; +} + + + +template +bool +FE_FaceQ<1,spacedim>::hp_constraints_are_implemented () const +{ + return true; +} + + + +template +FiniteElementDomination::Domination +FE_FaceQ<1,spacedim>:: +compare_for_face_domination (const FiniteElement<1,spacedim> &fe_other) const +{ + return FiniteElementDomination::no_requirements; +} + + + +template +Table<2,bool> +FE_FaceQ<1,spacedim>::get_constant_modes () const +{ + Table<2,bool> constant_modes(1, this->dofs_per_cell); + for (unsigned int i=0; idofs_per_cell; ++i) + constant_modes(0,i) = true; + return constant_modes; +} + + + +template +UpdateFlags +FE_FaceQ<1,spacedim>::update_once (const UpdateFlags) const +{ + return update_default; +} + + + +template +UpdateFlags +FE_FaceQ<1,spacedim>::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = flags & update_values; + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_hessians) + out |= update_hessians | update_covariant_transformation; + if (flags & update_cell_normal_vectors) + out |= update_cell_normal_vectors | update_JxW_values; + + return out; +} + + + +template +typename Mapping<1,spacedim>::InternalDataBase * +FE_FaceQ<1,spacedim>::get_data ( + const UpdateFlags, + const Mapping<1,spacedim> &, + const Quadrature<1> &) const +{ + return new typename Mapping<1,spacedim>::InternalDataBase; +} + + +template +typename Mapping<1,spacedim>::InternalDataBase * +FE_FaceQ<1,spacedim>::get_face_data ( + const UpdateFlags update_flags, + const Mapping<1,spacedim> &, + const Quadrature<0>& quadrature) const +{ + // generate a new data object and initialize some fields + typename Mapping<1,spacedim>::InternalDataBase *data = + new typename Mapping<1,spacedim>::InternalDataBase; + + // check what needs to be initialized only once and what on every + // cell/face/subface we visit + data->update_once = update_once(update_flags); + data->update_each = update_each(update_flags); + data->update_flags = data->update_once | data->update_each; + + const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.size(); + AssertDimension(n_q_points, 1); + + // No derivatives of this element are implemented. + if (flags & update_gradients || flags & update_hessians) + { + Assert(false, ExcNotImplemented()); + } + + return data; +} + + + +template +typename Mapping<1,spacedim>::InternalDataBase * +FE_FaceQ<1,spacedim>::get_subface_data ( + const UpdateFlags flags, + const Mapping<1,spacedim> &mapping, + const Quadrature<0>& quadrature) const +{ + return get_face_data (flags, mapping, quadrature); +} + + + +template +void +FE_FaceQ<1,spacedim>::fill_fe_values +(const Mapping<1,spacedim> &, + const typename Triangulation<1,spacedim>::cell_iterator &, + const Quadrature<1> &, + typename Mapping<1,spacedim>::InternalDataBase &, + typename Mapping<1,spacedim>::InternalDataBase &, + FEValuesData<1,spacedim> &, + CellSimilarity::Similarity &) const +{ + // Do nothing, since we do not have values in the interior +} + + + +template +void +FE_FaceQ<1,spacedim>::fill_fe_face_values ( + const Mapping<1,spacedim> &, + const typename Triangulation<1,spacedim>::cell_iterator &, + const unsigned int face, + const Quadrature<0>& quadrature, + typename Mapping<1,spacedim>::InternalDataBase &, + typename Mapping<1,spacedim>::InternalDataBase &fedata, + FEValuesData<1,spacedim> &data) const +{ + const UpdateFlags flags(fedata.update_once | fedata.update_each); + + const unsigned int foffset = face; + if (flags & update_values) + { + for (unsigned int k=0; kdofs_per_cell; ++k) + data.shape_values(k,0) = 0.; + data.shape_values(foffset,0) = 1; + } +} + + +template +void +FE_FaceQ<1,spacedim>::fill_fe_subface_values ( + const Mapping<1,spacedim> &, + const typename Triangulation<1,spacedim>::cell_iterator &, + const unsigned int , + const unsigned int , + const Quadrature<0>& , + typename Mapping<1,spacedim>::InternalDataBase &, + typename Mapping<1,spacedim>::InternalDataBase &, + FEValuesData<1,spacedim> &) const +{ + Assert(false, ExcMessage("Should not fille subface values in 1D")); +} + + + // --------------------------------------- FE_FaceP -------------------------- template @@ -312,6 +576,7 @@ FE_FaceP::FE_FaceP (const unsigned int degree) {} + template FiniteElement * FE_FaceP::clone() const @@ -320,6 +585,7 @@ FE_FaceP::clone() const } + template std::string FE_FaceP::get_name () const @@ -525,6 +791,29 @@ FE_FaceP::get_constant_modes () const +template +FE_FaceP<1,spacedim>::FE_FaceP (const unsigned int degree) + : + FE_FaceQ<1,spacedim> (degree) +{} + + + +template +std::string +FE_FaceP<1,spacedim>::get_name () const +{ + // note that the FETools::get_fe_from_name function depends on the + // particular format of the string this function returns, so they have to be + // kept in synch + std::ostringstream namebuf; + namebuf << "FE_FaceP<1>(" << this->degree << ")"; + + return namebuf.str(); +} + + + // explicit instantiations #include "fe_face.inst" diff --git a/deal.II/source/fe/fe_face.inst.in b/deal.II/source/fe/fe_face.inst.in index 513eb0030b..c1af401441 100644 --- a/deal.II/source/fe/fe_face.inst.in +++ b/deal.II/source/fe/fe_face.inst.in @@ -22,8 +22,8 @@ for (deal_II_dimension : DIMENSIONS) template class FE_PolyFace >; template class FE_PolyFace, deal_II_dimension>; //template class FE_PolyFace, deal_II_dimension>; - template class FE_FaceQ; - template class FE_FaceP; #endif + template class FE_FaceQ; + template class FE_FaceP; } diff --git a/tests/bits/step-51.cc b/tests/bits/step-51.cc new file mode 100644 index 0000000000..ed203d6d76 --- /dev/null +++ b/tests/bits/step-51.cc @@ -0,0 +1,1005 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2005 - 2013 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. +// +// --------------------------------------------------------------------- + + + +// tests step-51 (without output, without ChunkSparseMatrix, without +// adaptivity) in 1d and 2d + + +#include "../tests.h" +#include +#include +std::ofstream logfile("output"); + + +#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 + + +namespace Step51 +{ + + using namespace dealii; + + template + class SolutionBase + { + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; + }; + + + template <> + const Point<1> + SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] + = { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) + }; + + + template <> + const Point<2> + SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] + = { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) + }; + + template <> + const Point<3> + SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers] + = { Point<3>(-0.5, +0.5, 0.25), + Point<3>(-0.6, -0.5, -0.125), + Point<3>(+0.5, -0.5, 0.5) + }; + + template + const double SolutionBase::width = 1./5.; + + + template + class Solution : public Function, + protected SolutionBase + { + public: + Solution () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; + + + + template + double Solution::value (const Point &p, + const unsigned int) const + { + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + return_value += std::exp(-x_minus_xi.square() / + (this->width * this->width)); + } + + return return_value / + Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); + } + + + + template + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const + { + Tensor<1,dim> return_value; + + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += (-2 / (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width)) * + x_minus_xi); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * + this->width); + } + + + + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} + + virtual void vector_value (const Point &p, + Vector &v) const; + }; + + template + void SolutionAndGradient::vector_value (const Point &p, + Vector &v) const + { + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} + + virtual Tensor<1,dim> value (const Point &p) const; + }; + + + + template + Tensor<1,dim> + ConvectionVelocity::value(const Point &p) const + { + Tensor<1,dim> convection; + switch (dim) + { + case 1: + convection[0] = 1; + break; + case 2: + convection[0] = p[1]; + convection[1] = -p[0]; + break; + case 3: + convection[0] = p[1]; + convection[1] = -p[0]; + convection[2] = 1; + break; + default: + Assert(false, ExcNotImplemented()); + } + return convection; + } + + + + template + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + private: + const ConvectionVelocity convection_velocity; + }; + + + template + double RightHandSide::value (const Point &p, + const unsigned int) const + { + Tensor<1,dim> convection = convection_velocity.value(p); + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += + ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/ + (this->width * this->width)) / + (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width))); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } + + + template + class HDG + { + public: + HDG (const unsigned int degree); + void run (); + + private: + + void setup_system (); + void assemble_system (const bool reconstruct_trace = false); + void solve (); + void postprocess (); + void refine_grid (const unsigned int cylce); + void output_results (const unsigned int cycle); + + struct PerTaskData; + struct ScratchData; + + struct PostProcessScratchData; + + void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data); + + void copy_local_to_global(const PerTaskData &data); + + void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &empty_data); + + + Triangulation triangulation; + + FESystem fe_local; + DoFHandler dof_handler_local; + Vector solution_local; + + FE_FaceQ fe; + DoFHandler dof_handler; + Vector solution; + Vector system_rhs; + + FE_DGQ fe_u_post; + DoFHandler dof_handler_u_post; + Vector solution_u_post; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConvergenceTable convergence_table; + }; + + + template + HDG::HDG (const unsigned int degree) : + fe_local (FE_DGQ(degree), dim+1), + dof_handler_local (triangulation), + fe (degree), + dof_handler (triangulation), + fe_u_post (degree+1), + dof_handler_u_post (triangulation) + {} + + + + template + void + HDG::setup_system () + { + dof_handler_local.distribute_dofs(fe_local); + dof_handler.distribute_dofs(fe); + dof_handler_u_post.distribute_dofs(fe_u_post); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + solution_local.reinit (dof_handler_local.n_dofs()); + solution_u_post.reinit (dof_handler_u_post.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + typename FunctionMap::type boundary_functions; + Solution solution_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, + boundary_functions, + QGauss(fe.degree+1), + constraints); + constraints.close (); + + { + CompressedSimpleSparsityPattern csp (dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + sparsity_pattern.copy_from(csp); + } + system_matrix.reinit (sparsity_pattern); + } + + + + template + struct HDG::PerTaskData + { + FullMatrix cell_matrix; + Vector cell_vector; + std::vector dof_indices; + + bool trace_reconstruct; + + PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) + : cell_matrix(n_dofs, n_dofs), + cell_vector(n_dofs), + dof_indices(n_dofs), + trace_reconstruct(trace_reconstruct) + {} + }; + + + + template + struct HDG::ScratchData + { + FEValues fe_values_local; + FEFaceValues fe_face_values_local; + FEFaceValues fe_face_values; + + FullMatrix ll_matrix; + FullMatrix lf_matrix; + FullMatrix fl_matrix; + FullMatrix tmp_matrix; + Vector l_rhs; + Vector tmp_rhs; + + std::vector > q_phi; + std::vector q_phi_div; + std::vector u_phi; + std::vector > u_phi_grad; + std::vector tr_phi; + std::vector trace_values; + + std::vector > fe_local_support_on_face; + std::vector > fe_support_on_face; + + ConvectionVelocity convection_velocity; + RightHandSide right_hand_side; + const Solution exact_solution; + + ScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const QGauss &face_quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags local_face_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags), + fe_face_values (fe, face_quadrature_formula, flags), + ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell), + lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell), + fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + l_rhs (fe_local.dofs_per_cell), + tmp_rhs (fe_local.dofs_per_cell), + q_phi (fe_local.dofs_per_cell), + q_phi_div (fe_local.dofs_per_cell), + u_phi (fe_local.dofs_per_cell), + u_phi_grad (fe_local.dofs_per_cell), + tr_phi (fe.dofs_per_cell), + trace_values(face_quadrature_formula.size()), + fe_local_support_on_face(GeometryInfo::faces_per_cell), + fe_support_on_face(GeometryInfo::faces_per_cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::faces_per_cell; ++face) + for (unsigned int i=0; i + struct HDG::PostProcessScratchData + { + FEValues fe_values_local; + FEValues fe_values; + + std::vector u_values; + std::vector > u_gradients; + FullMatrix cell_matrix; + + Vector cell_rhs; + Vector cell_sol; + + PostProcessScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_values (fe, quadrature_formula, flags), + u_values (quadrature_formula.size()), + u_gradients (quadrature_formula.size()), + cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell), + cell_rhs (fe.dofs_per_cell), + cell_sol (fe.dofs_per_cell) + {} + + PostProcessScratchData(const PostProcessScratchData &sd) + : + fe_values_local (sd.fe_values_local.get_fe(), + sd.fe_values_local.get_quadrature(), + sd.fe_values_local.get_update_flags()), + fe_values (sd.fe_values.get_fe(), + sd.fe_values.get_quadrature(), + sd.fe_values.get_update_flags()), + u_values (sd.u_values), + u_gradients (sd.u_gradients), + cell_matrix (sd.cell_matrix), + cell_rhs (sd.cell_rhs), + cell_sol (sd.cell_sol) + {} + }; + + + + template + void + HDG::assemble_system (const bool trace_reconstruct) + { + const QGauss quadrature_formula(fe.degree+1); + const QGauss face_quadrature_formula(fe.degree+1); + + const UpdateFlags local_flags (update_values | update_gradients | + update_JxW_values | update_quadrature_points); + + const UpdateFlags local_face_flags (update_values); + + const UpdateFlags flags ( update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + PerTaskData task_data (fe.dofs_per_cell, + trace_reconstruct); + ScratchData scratch (fe, fe_local, + quadrature_formula, + face_quadrature_formula, + local_flags, + local_face_flags, + flags); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &HDG::assemble_system_one_cell, + &HDG::copy_local_to_global, + scratch, + task_data); + } + + + + template + void + HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data) + { + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size(); + const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size(); + + const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; + + const FEValuesExtractors::Vector fluxes (0); + const FEValuesExtractors::Scalar scalar (dim); + + scratch.ll_matrix = 0; + scratch.l_rhs = 0; + if (!task_data.trace_reconstruct) + { + scratch.lf_matrix = 0; + scratch.fl_matrix = 0; + task_data.cell_matrix = 0; + task_data.cell_vector = 0; + } + scratch.fe_values_local.reinit (loc_cell); + + for (unsigned int q=0; q convection + = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); + const double JxW = scratch.fe_values_local.JxW(q); + for (unsigned int k=0; k::faces_per_cell; ++face) + { + scratch.fe_face_values_local.reinit(loc_cell, face); + scratch.fe_face_values.reinit(cell, face); + + if (task_data.trace_reconstruct) + scratch.fe_face_values.get_function_values (solution, scratch.trace_values); + + for (unsigned int q=0; q quadrature_point = + scratch.fe_face_values.quadrature_point(q); + const Point normal = scratch.fe_face_values.normal_vector(q); + const Tensor<1,dim> convection + = scratch.convection_velocity.value(quadrature_point); + + const double tau_stab = (5. + + std::abs(convection * normal)); + + for (unsigned int k=0; kface(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) + { + const double neumann_value = + - scratch.exact_solution.gradient (quadrature_point) * normal + + convection * normal * scratch.exact_solution.value(quadrature_point); + for (unsigned int i=0; iget_dof_indices(task_data.dof_indices); + } + else + { + scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); + loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); + } + } + + + + template + void HDG::copy_local_to_global(const PerTaskData &data) + { + if (data.trace_reconstruct == false) + constraints.distribute_local_to_global (data.cell_matrix, + data.cell_vector, + data.dof_indices, + system_matrix, system_rhs); + } + + + + template + void HDG::solve () + { + SolverControl solver_control (system_matrix.m()*10, + 1e-11*system_rhs.l2_norm()); + SolverBicgstab<> solver (solver_control, false); + solver.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); + + system_matrix.clear(); + + constraints.distribute(solution); + + assemble_system(true); + } + + + + + template + void + HDG::postprocess() + { + { + const QGauss quadrature_formula(fe_u_post.degree+1); + const UpdateFlags local_flags (update_values); + const UpdateFlags flags ( update_values | update_gradients | + update_JxW_values); + + PostProcessScratchData scratch (fe_u_post, fe_local, + quadrature_formula, + local_flags, + flags); + + WorkStream::run(dof_handler_u_post.begin_active(), + dof_handler_u_post.end(), + std_cxx1x::bind (&HDG::postprocess_one_cell, + std_cxx1x::ref(*this), + std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), + std_cxx1x::function(), + scratch, + 0U); + } + + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &value_select); + const double L2_error = difference_per_cell.l2_norm(); + + ComponentSelectFunction gradient_select (std::pair(0, dim), + dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &gradient_select); + const double grad_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler_u_post, + solution_u_post, + Solution(), + difference_per_cell, + QGauss(fe.degree+3), + VectorTools::L2_norm); + const double post_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", triangulation.n_active_cells()); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("val L2", L2_error); + convergence_table.add_value("grad L2", grad_error); + convergence_table.add_value("val L2-post", post_error); + } + + + + template + void + HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &) + { + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + scratch.fe_values_local.reinit (loc_cell); + scratch.fe_values.reinit(cell); + + FEValuesExtractors::Vector fluxes(0); + FEValuesExtractors::Scalar scalar(dim); + + const unsigned int n_q_points = scratch.fe_values.get_quadrature().size(); + const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell; + + scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values); + scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients); + + double sum = 0; + for (unsigned int i=1; idistribute_local_to_global(scratch.cell_sol, solution_u_post); + } + + + + template + void HDG::output_results (const unsigned int cycle) + { + // not included in test + } + + + + template + void HDG::refine_grid (const unsigned int cycle) + { + // only global refinement + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); + + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + for (unsigned int d=0; dface(face)->center()(d) - (1)) < 1e-12)) + cell->face(face)->set_boundary_indicator (1); + } + + + + template + void HDG::run () + { + for (unsigned int cycle=0; cycle<14-4*dim; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + refine_grid (cycle); + setup_system (); + assemble_system (false); + solve (); + postprocess(); + output_results (cycle); + } + + convergence_table.set_precision("val L2", 3); + convergence_table.set_scientific("val L2", true); + convergence_table.set_precision("grad L2", 3); + convergence_table.set_scientific("grad L2", true); + convergence_table.set_precision("val L2-post", 3); + convergence_table.set_scientific("val L2-post", true); + + convergence_table + .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table.write_text(deallog.get_file_stream()); + } + +} // end of namespace Step51 + + + + +int main () +{ + deallog << std::setprecision(6); + logfile << std::setprecision(6); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + { + Step51::HDG<1> hdg_problem (1); + hdg_problem.run (); + } + + { + Step51::HDG<1> hdg_problem (3); + hdg_problem.run (); + } + + { + Step51::HDG<2> hdg_problem (1); + hdg_problem.run (); + } + + { + Step51::HDG<2> hdg_problem (2); + hdg_problem.run (); + } + + return 0; +} diff --git a/tests/bits/step-51.output b/tests/bits/step-51.output new file mode 100644 index 0000000000..768c5da026 --- /dev/null +++ b/tests/bits/step-51.output @@ -0,0 +1,165 @@ + +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 9 +DEAL:Bicgstab::Starting value 8.19443 +DEAL:Bicgstab::Convergence step 8 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 13 +DEAL:Bicgstab::Starting value 13.9539 +DEAL:Bicgstab::Convergence step 12 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 17 +DEAL:Bicgstab::Starting value 13.0731 +DEAL:Bicgstab::Convergence step 16 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 25 +DEAL:Bicgstab::Starting value 12.1492 +DEAL:Bicgstab::Convergence step 24 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 33 +DEAL:Bicgstab::Starting value 11.0492 +DEAL:Bicgstab::Convergence step 34 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 49 +DEAL:Bicgstab::Starting value 9.35310 +DEAL:Bicgstab::Convergence step 51 value 0 +DEAL::Cycle 6: +DEAL:: Number of degrees of freedom: 65 +DEAL:Bicgstab::Starting value 8.20475 +DEAL:Bicgstab::Convergence step 74 value 0 +DEAL::Cycle 7: +DEAL:: Number of degrees of freedom: 97 +DEAL:Bicgstab::Starting value 6.76136 +DEAL:Bicgstab::Convergence step 104 value 0 +DEAL::Cycle 8: +DEAL:: Number of degrees of freedom: 129 +DEAL:Bicgstab::Starting value 5.87456 +DEAL:Bicgstab::Convergence step 137 value 0 +DEAL::Cycle 9: +DEAL:: Number of degrees of freedom: 193 +DEAL:Bicgstab::Starting value 4.80772 +DEAL:Bicgstab::Convergence step 213 value 0 +cells dofs val L2 grad L2 val L2-post +8 9 4.080e-01 - 1.722e+00 - 3.362e-01 - +12 13 1.659e-01 2.22 8.766e-01 1.67 2.151e-02 6.78 +16 17 7.996e-02 2.54 4.531e-01 2.29 1.125e-02 2.25 +24 25 3.705e-02 1.90 2.028e-01 1.98 3.372e-03 2.97 +32 33 2.125e-02 1.93 1.140e-01 2.00 1.430e-03 2.98 +48 49 9.610e-03 1.96 5.048e-02 2.01 4.242e-04 3.00 +64 65 5.449e-03 1.97 2.831e-02 2.01 1.787e-04 3.00 +96 97 2.440e-03 1.98 1.253e-02 2.01 5.281e-05 3.01 +128 129 1.377e-03 1.99 7.035e-03 2.01 2.224e-05 3.01 +192 193 6.141e-04 1.99 3.119e-03 2.01 6.573e-06 3.01 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 9 +DEAL:Bicgstab::Starting value 7.96477 +DEAL:Bicgstab::Convergence step 8 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 13 +DEAL:Bicgstab::Starting value 13.4076 +DEAL:Bicgstab::Convergence step 12 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 17 +DEAL:Bicgstab::Starting value 12.9549 +DEAL:Bicgstab::Convergence step 16 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 25 +DEAL:Bicgstab::Starting value 12.1266 +DEAL:Bicgstab::Convergence step 24 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 33 +DEAL:Bicgstab::Starting value 11.0426 +DEAL:Bicgstab::Convergence step 34 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 49 +DEAL:Bicgstab::Starting value 9.35200 +DEAL:Bicgstab::Convergence step 50 value 0 +DEAL::Cycle 6: +DEAL:: Number of degrees of freedom: 65 +DEAL:Bicgstab::Starting value 8.20444 +DEAL:Bicgstab::Convergence step 73 value 0 +DEAL::Cycle 7: +DEAL:: Number of degrees of freedom: 97 +DEAL:Bicgstab::Starting value 6.76131 +DEAL:Bicgstab::Convergence step 105 value 0 +DEAL::Cycle 8: +DEAL:: Number of degrees of freedom: 129 +DEAL:Bicgstab::Starting value 5.87455 +DEAL:Bicgstab::Convergence step 143 value 0 +DEAL::Cycle 9: +DEAL:: Number of degrees of freedom: 193 +DEAL:Bicgstab::Starting value 4.80772 +DEAL:Bicgstab::Convergence step 225 value 0 +cells dofs val L2 grad L2 val L2-post +8 9 2.141e-02 - 1.398e-01 - 3.487e-03 - +12 13 6.701e-03 2.87 3.811e-02 3.21 4.976e-04 4.80 +16 17 1.693e-03 4.78 9.999e-03 4.65 1.085e-04 5.29 +24 25 3.501e-04 3.89 2.009e-03 3.96 1.440e-05 4.98 +32 33 1.126e-04 3.94 6.373e-04 3.99 3.412e-06 5.01 +48 49 2.256e-05 3.97 1.259e-04 4.00 4.472e-07 5.01 +64 65 7.179e-06 3.98 3.980e-05 4.00 1.057e-07 5.01 +96 97 1.425e-06 3.99 7.848e-06 4.00 1.388e-08 5.01 +128 129 4.519e-07 3.99 2.481e-06 4.00 3.302e-09 4.99 +192 193 8.946e-08 3.99 4.894e-07 4.00 7.398e-10 3.69 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 80 +DEAL:Bicgstab::Starting value 12.2150 +DEAL:Bicgstab::Convergence step 37 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 168 +DEAL:Bicgstab::Starting value 8.50278 +DEAL:Bicgstab::Convergence step 48 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 288 +DEAL:Bicgstab::Starting value 10.5333 +DEAL:Bicgstab::Convergence step 62 value 1.01190e-10 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 624 +DEAL:Bicgstab::Starting value 8.10478 +DEAL:Bicgstab::Convergence step 98 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 1088 +DEAL:Bicgstab::Starting value 6.61947 +DEAL:Bicgstab::Convergence step 120 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 2400 +DEAL:Bicgstab::Starting value 4.75284 +DEAL:Bicgstab::Convergence step 181 value 0 +cells dofs val L2 grad L2 val L2-post +16 80 1.066e+01 - 1.571e+01 - 1.055e+01 - +36 168 3.727e+00 2.59 7.998e+00 1.67 3.660e+00 2.61 +64 288 8.086e-01 5.31 4.228e+00 2.22 4.492e-01 7.29 +144 624 2.730e-01 2.68 1.867e+00 2.02 6.110e-02 4.92 +256 1088 1.493e-01 2.10 1.046e+00 2.01 2.878e-02 2.62 +576 2400 6.965e-02 1.88 4.847e-01 1.90 9.202e-03 2.81 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 120 +DEAL:Bicgstab::Starting value 10.6959 +DEAL:Bicgstab::Convergence step 43 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 252 +DEAL:Bicgstab::Starting value 13.2620 +DEAL:Bicgstab::Convergence step 58 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 432 +DEAL:Bicgstab::Starting value 11.2567 +DEAL:Bicgstab::Convergence step 74 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 936 +DEAL:Bicgstab::Starting value 9.46542 +DEAL:Bicgstab::Convergence step 105 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 1632 +DEAL:Bicgstab::Starting value 7.48176 +DEAL:Bicgstab::Convergence step 138 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 3600 +DEAL:Bicgstab::Starting value 5.13012 +DEAL:Bicgstab::Convergence step 189 value 0 +cells dofs val L2 grad L2 val L2-post +16 120 7.645e-01 - 5.413e+00 - 2.667e-01 - +36 252 4.745e-01 1.18 2.375e+00 2.03 3.501e-01 -0.67 +64 432 9.635e-02 5.54 7.772e-01 3.88 5.032e-02 6.74 +144 936 3.487e-02 2.51 2.718e-01 2.59 5.339e-03 5.53 +256 1632 1.920e-02 2.08 1.390e-01 2.33 1.731e-03 3.92 +576 3600 6.021e-03 2.86 4.291e-02 2.90 3.431e-04 3.99 diff --git a/tests/bits/step-51p.cc b/tests/bits/step-51p.cc new file mode 100644 index 0000000000..fae43af7b6 --- /dev/null +++ b/tests/bits/step-51p.cc @@ -0,0 +1,1005 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2005 - 2013 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. +// +// --------------------------------------------------------------------- + + + +// tests step-51 (without output, without ChunkSparseMatrix, without +// adaptivity) in 1d and 2d with FE_DGP/FE_FaceP + + +#include "../tests.h" +#include +#include +std::ofstream logfile("output"); + + +#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 + + +namespace Step51 +{ + + using namespace dealii; + + template + class SolutionBase + { + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; + }; + + + template <> + const Point<1> + SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] + = { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) + }; + + + template <> + const Point<2> + SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] + = { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) + }; + + template <> + const Point<3> + SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers] + = { Point<3>(-0.5, +0.5, 0.25), + Point<3>(-0.6, -0.5, -0.125), + Point<3>(+0.5, -0.5, 0.5) + }; + + template + const double SolutionBase::width = 1./5.; + + + template + class Solution : public Function, + protected SolutionBase + { + public: + Solution () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; + + + + template + double Solution::value (const Point &p, + const unsigned int) const + { + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + return_value += std::exp(-x_minus_xi.square() / + (this->width * this->width)); + } + + return return_value / + Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); + } + + + + template + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const + { + Tensor<1,dim> return_value; + + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += (-2 / (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width)) * + x_minus_xi); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * + this->width); + } + + + + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} + + virtual void vector_value (const Point &p, + Vector &v) const; + }; + + template + void SolutionAndGradient::vector_value (const Point &p, + Vector &v) const + { + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} + + virtual Tensor<1,dim> value (const Point &p) const; + }; + + + + template + Tensor<1,dim> + ConvectionVelocity::value(const Point &p) const + { + Tensor<1,dim> convection; + switch (dim) + { + case 1: + convection[0] = 1; + break; + case 2: + convection[0] = p[1]; + convection[1] = -p[0]; + break; + case 3: + convection[0] = p[1]; + convection[1] = -p[0]; + convection[2] = 1; + break; + default: + Assert(false, ExcNotImplemented()); + } + return convection; + } + + + + template + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + private: + const ConvectionVelocity convection_velocity; + }; + + + template + double RightHandSide::value (const Point &p, + const unsigned int) const + { + Tensor<1,dim> convection = convection_velocity.value(p); + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += + ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/ + (this->width * this->width)) / + (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width))); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } + + + template + class HDG + { + public: + HDG (const unsigned int degree); + void run (); + + private: + + void setup_system (); + void assemble_system (const bool reconstruct_trace = false); + void solve (); + void postprocess (); + void refine_grid (const unsigned int cylce); + void output_results (const unsigned int cycle); + + struct PerTaskData; + struct ScratchData; + + struct PostProcessScratchData; + + void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data); + + void copy_local_to_global(const PerTaskData &data); + + void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &empty_data); + + + Triangulation triangulation; + + FESystem fe_local; + DoFHandler dof_handler_local; + Vector solution_local; + + FE_FaceP fe; + DoFHandler dof_handler; + Vector solution; + Vector system_rhs; + + FE_DGP fe_u_post; + DoFHandler dof_handler_u_post; + Vector solution_u_post; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConvergenceTable convergence_table; + }; + + + template + HDG::HDG (const unsigned int degree) : + fe_local (FE_DGP(degree), dim+1), + dof_handler_local (triangulation), + fe (degree), + dof_handler (triangulation), + fe_u_post (degree+1), + dof_handler_u_post (triangulation) + {} + + + + template + void + HDG::setup_system () + { + dof_handler_local.distribute_dofs(fe_local); + dof_handler.distribute_dofs(fe); + dof_handler_u_post.distribute_dofs(fe_u_post); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + solution_local.reinit (dof_handler_local.n_dofs()); + solution_u_post.reinit (dof_handler_u_post.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + typename FunctionMap::type boundary_functions; + Solution solution_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, + boundary_functions, + QGauss(fe.degree+1), + constraints); + constraints.close (); + + { + CompressedSimpleSparsityPattern csp (dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + sparsity_pattern.copy_from(csp); + } + system_matrix.reinit (sparsity_pattern); + } + + + + template + struct HDG::PerTaskData + { + FullMatrix cell_matrix; + Vector cell_vector; + std::vector dof_indices; + + bool trace_reconstruct; + + PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) + : cell_matrix(n_dofs, n_dofs), + cell_vector(n_dofs), + dof_indices(n_dofs), + trace_reconstruct(trace_reconstruct) + {} + }; + + + + template + struct HDG::ScratchData + { + FEValues fe_values_local; + FEFaceValues fe_face_values_local; + FEFaceValues fe_face_values; + + FullMatrix ll_matrix; + FullMatrix lf_matrix; + FullMatrix fl_matrix; + FullMatrix tmp_matrix; + Vector l_rhs; + Vector tmp_rhs; + + std::vector > q_phi; + std::vector q_phi_div; + std::vector u_phi; + std::vector > u_phi_grad; + std::vector tr_phi; + std::vector trace_values; + + std::vector > fe_local_support_on_face; + std::vector > fe_support_on_face; + + ConvectionVelocity convection_velocity; + RightHandSide right_hand_side; + const Solution exact_solution; + + ScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const QGauss &face_quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags local_face_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags), + fe_face_values (fe, face_quadrature_formula, flags), + ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell), + lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell), + fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + l_rhs (fe_local.dofs_per_cell), + tmp_rhs (fe_local.dofs_per_cell), + q_phi (fe_local.dofs_per_cell), + q_phi_div (fe_local.dofs_per_cell), + u_phi (fe_local.dofs_per_cell), + u_phi_grad (fe_local.dofs_per_cell), + tr_phi (fe.dofs_per_cell), + trace_values(face_quadrature_formula.size()), + fe_local_support_on_face(GeometryInfo::faces_per_cell), + fe_support_on_face(GeometryInfo::faces_per_cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::faces_per_cell; ++face) + for (unsigned int i=0; i + struct HDG::PostProcessScratchData + { + FEValues fe_values_local; + FEValues fe_values; + + std::vector u_values; + std::vector > u_gradients; + FullMatrix cell_matrix; + + Vector cell_rhs; + Vector cell_sol; + + PostProcessScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_values (fe, quadrature_formula, flags), + u_values (quadrature_formula.size()), + u_gradients (quadrature_formula.size()), + cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell), + cell_rhs (fe.dofs_per_cell), + cell_sol (fe.dofs_per_cell) + {} + + PostProcessScratchData(const PostProcessScratchData &sd) + : + fe_values_local (sd.fe_values_local.get_fe(), + sd.fe_values_local.get_quadrature(), + sd.fe_values_local.get_update_flags()), + fe_values (sd.fe_values.get_fe(), + sd.fe_values.get_quadrature(), + sd.fe_values.get_update_flags()), + u_values (sd.u_values), + u_gradients (sd.u_gradients), + cell_matrix (sd.cell_matrix), + cell_rhs (sd.cell_rhs), + cell_sol (sd.cell_sol) + {} + }; + + + + template + void + HDG::assemble_system (const bool trace_reconstruct) + { + const QGauss quadrature_formula(fe.degree+1); + const QGauss face_quadrature_formula(fe.degree+1); + + const UpdateFlags local_flags (update_values | update_gradients | + update_JxW_values | update_quadrature_points); + + const UpdateFlags local_face_flags (update_values); + + const UpdateFlags flags ( update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + PerTaskData task_data (fe.dofs_per_cell, + trace_reconstruct); + ScratchData scratch (fe, fe_local, + quadrature_formula, + face_quadrature_formula, + local_flags, + local_face_flags, + flags); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &HDG::assemble_system_one_cell, + &HDG::copy_local_to_global, + scratch, + task_data); + } + + + + template + void + HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data) + { + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size(); + const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size(); + + const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; + + const FEValuesExtractors::Vector fluxes (0); + const FEValuesExtractors::Scalar scalar (dim); + + scratch.ll_matrix = 0; + scratch.l_rhs = 0; + if (!task_data.trace_reconstruct) + { + scratch.lf_matrix = 0; + scratch.fl_matrix = 0; + task_data.cell_matrix = 0; + task_data.cell_vector = 0; + } + scratch.fe_values_local.reinit (loc_cell); + + for (unsigned int q=0; q convection + = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); + const double JxW = scratch.fe_values_local.JxW(q); + for (unsigned int k=0; k::faces_per_cell; ++face) + { + scratch.fe_face_values_local.reinit(loc_cell, face); + scratch.fe_face_values.reinit(cell, face); + + if (task_data.trace_reconstruct) + scratch.fe_face_values.get_function_values (solution, scratch.trace_values); + + for (unsigned int q=0; q quadrature_point = + scratch.fe_face_values.quadrature_point(q); + const Point normal = scratch.fe_face_values.normal_vector(q); + const Tensor<1,dim> convection + = scratch.convection_velocity.value(quadrature_point); + + const double tau_stab = (5. + + std::abs(convection * normal)); + + for (unsigned int k=0; kface(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) + { + const double neumann_value = + - scratch.exact_solution.gradient (quadrature_point) * normal + + convection * normal * scratch.exact_solution.value(quadrature_point); + for (unsigned int i=0; iget_dof_indices(task_data.dof_indices); + } + else + { + scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); + loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); + } + } + + + + template + void HDG::copy_local_to_global(const PerTaskData &data) + { + if (data.trace_reconstruct == false) + constraints.distribute_local_to_global (data.cell_matrix, + data.cell_vector, + data.dof_indices, + system_matrix, system_rhs); + } + + + + template + void HDG::solve () + { + SolverControl solver_control (system_matrix.m()*10, + 1e-11*system_rhs.l2_norm()); + SolverBicgstab<> solver (solver_control, false); + solver.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); + + system_matrix.clear(); + + constraints.distribute(solution); + + assemble_system(true); + } + + + + + template + void + HDG::postprocess() + { + { + const QGauss quadrature_formula(fe_u_post.degree+1); + const UpdateFlags local_flags (update_values); + const UpdateFlags flags ( update_values | update_gradients | + update_JxW_values); + + PostProcessScratchData scratch (fe_u_post, fe_local, + quadrature_formula, + local_flags, + flags); + + WorkStream::run(dof_handler_u_post.begin_active(), + dof_handler_u_post.end(), + std_cxx1x::bind (&HDG::postprocess_one_cell, + std_cxx1x::ref(*this), + std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), + std_cxx1x::function(), + scratch, + 0U); + } + + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &value_select); + const double L2_error = difference_per_cell.l2_norm(); + + ComponentSelectFunction gradient_select (std::pair(0, dim), + dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &gradient_select); + const double grad_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler_u_post, + solution_u_post, + Solution(), + difference_per_cell, + QGauss(fe.degree+3), + VectorTools::L2_norm); + const double post_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", triangulation.n_active_cells()); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("val L2", L2_error); + convergence_table.add_value("grad L2", grad_error); + convergence_table.add_value("val L2-post", post_error); + } + + + + template + void + HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &) + { + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + scratch.fe_values_local.reinit (loc_cell); + scratch.fe_values.reinit(cell); + + FEValuesExtractors::Vector fluxes(0); + FEValuesExtractors::Scalar scalar(dim); + + const unsigned int n_q_points = scratch.fe_values.get_quadrature().size(); + const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell; + + scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values); + scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients); + + double sum = 0; + for (unsigned int i=1; idistribute_local_to_global(scratch.cell_sol, solution_u_post); + } + + + + template + void HDG::output_results (const unsigned int cycle) + { + // not included in test + } + + + + template + void HDG::refine_grid (const unsigned int cycle) + { + // only global refinement + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); + + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + for (unsigned int d=0; dface(face)->center()(d) - (1)) < 1e-12)) + cell->face(face)->set_boundary_indicator (1); + } + + + + template + void HDG::run () + { + for (unsigned int cycle=0; cycle<14-4*dim; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + refine_grid (cycle); + setup_system (); + assemble_system (false); + solve (); + postprocess(); + output_results (cycle); + } + + convergence_table.set_precision("val L2", 3); + convergence_table.set_scientific("val L2", true); + convergence_table.set_precision("grad L2", 3); + convergence_table.set_scientific("grad L2", true); + convergence_table.set_precision("val L2-post", 3); + convergence_table.set_scientific("val L2-post", true); + + convergence_table + .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table.write_text(deallog.get_file_stream()); + } + +} // end of namespace Step51 + + + + +int main () +{ + deallog << std::setprecision(6); + logfile << std::setprecision(6); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + { + Step51::HDG<1> hdg_problem (1); + hdg_problem.run (); + } + + { + Step51::HDG<1> hdg_problem (3); + hdg_problem.run (); + } + + { + Step51::HDG<2> hdg_problem (1); + hdg_problem.run (); + } + + { + Step51::HDG<2> hdg_problem (2); + hdg_problem.run (); + } + + return 0; +} diff --git a/tests/bits/step-51p.output b/tests/bits/step-51p.output new file mode 100644 index 0000000000..a0f2912453 --- /dev/null +++ b/tests/bits/step-51p.output @@ -0,0 +1,165 @@ + +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 9 +DEAL:Bicgstab::Starting value 8.19443 +DEAL:Bicgstab::Convergence step 8 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 13 +DEAL:Bicgstab::Starting value 13.9539 +DEAL:Bicgstab::Convergence step 12 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 17 +DEAL:Bicgstab::Starting value 13.0731 +DEAL:Bicgstab::Convergence step 16 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 25 +DEAL:Bicgstab::Starting value 12.1492 +DEAL:Bicgstab::Convergence step 24 value 1.07658e-10 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 33 +DEAL:Bicgstab::Starting value 11.0492 +DEAL:Bicgstab::Convergence step 34 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 49 +DEAL:Bicgstab::Starting value 9.35310 +DEAL:Bicgstab::Convergence step 51 value 0 +DEAL::Cycle 6: +DEAL:: Number of degrees of freedom: 65 +DEAL:Bicgstab::Starting value 8.20475 +DEAL:Bicgstab::Convergence step 75 value 0 +DEAL::Cycle 7: +DEAL:: Number of degrees of freedom: 97 +DEAL:Bicgstab::Starting value 6.76136 +DEAL:Bicgstab::Convergence step 104 value 0 +DEAL::Cycle 8: +DEAL:: Number of degrees of freedom: 129 +DEAL:Bicgstab::Starting value 5.87456 +DEAL:Bicgstab::Convergence step 142 value 0 +DEAL::Cycle 9: +DEAL:: Number of degrees of freedom: 193 +DEAL:Bicgstab::Starting value 4.80772 +DEAL:Bicgstab::Convergence step 211 value 0 +cells dofs val L2 grad L2 val L2-post +8 9 4.080e-01 - 1.722e+00 - 3.362e-01 - +12 13 1.659e-01 2.22 8.766e-01 1.67 2.151e-02 6.78 +16 17 7.996e-02 2.54 4.531e-01 2.29 1.125e-02 2.25 +24 25 3.705e-02 1.90 2.028e-01 1.98 3.372e-03 2.97 +32 33 2.125e-02 1.93 1.140e-01 2.00 1.430e-03 2.98 +48 49 9.610e-03 1.96 5.048e-02 2.01 4.242e-04 3.00 +64 65 5.449e-03 1.97 2.831e-02 2.01 1.787e-04 3.00 +96 97 2.440e-03 1.98 1.253e-02 2.01 5.281e-05 3.01 +128 129 1.377e-03 1.99 7.035e-03 2.01 2.224e-05 3.01 +192 193 6.141e-04 1.99 3.119e-03 2.01 6.573e-06 3.01 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 9 +DEAL:Bicgstab::Starting value 7.96477 +DEAL:Bicgstab::Convergence step 8 value 0 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 13 +DEAL:Bicgstab::Starting value 13.4076 +DEAL:Bicgstab::Convergence step 12 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 17 +DEAL:Bicgstab::Starting value 12.9549 +DEAL:Bicgstab::Convergence step 16 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 25 +DEAL:Bicgstab::Starting value 12.1266 +DEAL:Bicgstab::Convergence step 25 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 33 +DEAL:Bicgstab::Starting value 11.0426 +DEAL:Bicgstab::Convergence step 33 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 49 +DEAL:Bicgstab::Starting value 9.35200 +DEAL:Bicgstab::Convergence step 49 value 0 +DEAL::Cycle 6: +DEAL:: Number of degrees of freedom: 65 +DEAL:Bicgstab::Starting value 8.20444 +DEAL:Bicgstab::Convergence step 71 value 0 +DEAL::Cycle 7: +DEAL:: Number of degrees of freedom: 97 +DEAL:Bicgstab::Starting value 6.76131 +DEAL:Bicgstab::Convergence step 111 value 0 +DEAL::Cycle 8: +DEAL:: Number of degrees of freedom: 129 +DEAL:Bicgstab::Starting value 5.87455 +DEAL:Bicgstab::Convergence step 144 value 0 +DEAL::Cycle 9: +DEAL:: Number of degrees of freedom: 193 +DEAL:Bicgstab::Starting value 4.80772 +DEAL:Bicgstab::Convergence step 206 value 0 +cells dofs val L2 grad L2 val L2-post +8 9 2.141e-02 - 1.398e-01 - 3.487e-03 - +12 13 6.701e-03 2.87 3.811e-02 3.21 4.976e-04 4.80 +16 17 1.693e-03 4.78 9.999e-03 4.65 1.085e-04 5.29 +24 25 3.501e-04 3.89 2.009e-03 3.96 1.440e-05 4.98 +32 33 1.126e-04 3.94 6.373e-04 3.99 3.412e-06 5.01 +48 49 2.256e-05 3.97 1.259e-04 4.00 4.472e-07 5.01 +64 65 7.179e-06 3.98 3.980e-05 4.00 1.057e-07 5.01 +96 97 1.425e-06 3.99 7.848e-06 4.00 1.388e-08 5.01 +128 129 4.519e-07 3.99 2.481e-06 4.00 3.400e-09 4.89 +192 193 8.946e-08 3.99 4.894e-07 4.00 1.008e-09 3.00 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 80 +DEAL:Bicgstab::Starting value 17.6593 +DEAL:Bicgstab::Convergence step 31 value 1.30726e-10 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 168 +DEAL:Bicgstab::Starting value 11.4720 +DEAL:Bicgstab::Convergence step 44 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 288 +DEAL:Bicgstab::Starting value 16.7917 +DEAL:Bicgstab::Convergence step 63 value 1.47629e-10 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 624 +DEAL:Bicgstab::Starting value 12.5049 +DEAL:Bicgstab::Convergence step 90 value 0 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 1088 +DEAL:Bicgstab::Starting value 9.93628 +DEAL:Bicgstab::Convergence step 120 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 2400 +DEAL:Bicgstab::Starting value 6.94024 +DEAL:Bicgstab::Convergence step 177 value 0 +cells dofs val L2 grad L2 val L2-post +16 80 1.043e+01 - 1.616e+01 - 1.034e+01 - +36 168 3.699e+00 2.56 8.081e+00 1.71 3.637e+00 2.58 +64 288 8.022e-01 5.31 4.671e+00 1.91 4.622e-01 7.17 +144 624 2.947e-01 2.47 2.255e+00 1.80 6.403e-02 4.88 +256 1088 1.638e-01 2.04 1.343e+00 1.80 2.937e-02 2.71 +576 2400 7.616e-02 1.89 6.522e-01 1.78 9.203e-03 2.86 +DEAL::Cycle 0: +DEAL:: Number of degrees of freedom: 120 +DEAL:Bicgstab::Starting value 19.9021 +DEAL:Bicgstab::Convergence step 39 value 1.15234e-10 +DEAL::Cycle 1: +DEAL:: Number of degrees of freedom: 252 +DEAL:Bicgstab::Starting value 15.8809 +DEAL:Bicgstab::Convergence step 56 value 0 +DEAL::Cycle 2: +DEAL:: Number of degrees of freedom: 432 +DEAL:Bicgstab::Starting value 16.3149 +DEAL:Bicgstab::Convergence step 71 value 0 +DEAL::Cycle 3: +DEAL:: Number of degrees of freedom: 936 +DEAL:Bicgstab::Starting value 12.6721 +DEAL:Bicgstab::Convergence step 109 value 1.21847e-10 +DEAL::Cycle 4: +DEAL:: Number of degrees of freedom: 1632 +DEAL:Bicgstab::Starting value 10.0470 +DEAL:Bicgstab::Convergence step 144 value 0 +DEAL::Cycle 5: +DEAL:: Number of degrees of freedom: 3600 +DEAL:Bicgstab::Starting value 6.97340 +DEAL:Bicgstab::Convergence step 198 value 0 +cells dofs val L2 grad L2 val L2-post +16 120 8.296e-01 - 6.543e+00 - 4.572e-01 - +36 252 5.045e-01 1.23 3.062e+00 1.87 3.570e-01 0.61 +64 432 1.736e-01 3.71 1.503e+00 2.47 6.468e-02 5.94 +144 936 5.584e-02 2.80 5.042e-01 2.69 1.006e-02 4.59 +256 1632 2.731e-02 2.49 2.274e-01 2.77 3.120e-03 4.07 +576 3600 8.582e-03 2.85 7.254e-02 2.82 6.213e-04 3.98 diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc index 7d00b83826..a9e16506a1 100644 --- a/tests/fe/fe_data_test.cc +++ b/tests/fe/fe_data_test.cc @@ -44,20 +44,6 @@ template void test_2d_3d (std::vector *> &fe_datas) { - // Face Q elements - fe_datas.push_back(new FE_FaceQ (0)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; - fe_datas.push_back(new FE_FaceQ (1)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; - fe_datas.push_back(new FE_FaceQ (3)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; - // Face P elements - fe_datas.push_back(new FE_FaceP (0)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; - fe_datas.push_back(new FE_FaceP (1)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; - fe_datas.push_back(new FE_FaceP (3)); - deallog << (*fe_datas.rbegin())->get_name() << std::endl; // Vector DG elements fe_datas.push_back( new FE_DGRaviartThomas(0)); @@ -136,6 +122,21 @@ void test_fe_datas() FE_Q (2), 1)); deallog << (*fe_datas.rbegin())->get_name() << std::endl; + // Face Q elements + fe_datas.push_back(new FE_FaceQ (0)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + fe_datas.push_back(new FE_FaceQ (1)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + fe_datas.push_back(new FE_FaceQ (3)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + // Face P elements + fe_datas.push_back(new FE_FaceP (0)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + fe_datas.push_back(new FE_FaceP (1)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + fe_datas.push_back(new FE_FaceP (3)); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + // Check vector elements in 2d and higher only test_2d_3d (fe_datas);