<h3>Specific improvements</h3>
<ol>
+ <li>New: FE_FaceQ and FE_FaceP now also work in 1D (with a single dof
+ on each vertex).
+ <br>
+ (Martin Kronbichler, 2014/02/11)
+
+ <li>Fixed: FE_DGQ::has_support_on_face returned a wrong number for element
+ degree larger than 1 in 1D. This is now fixed.
+ <br>
+ (Martin Kronbichler, 2014/02/10)
+
<li>Changed: DerivativeApproximation used to be a class that only had
static members. It is now a namespace.
<br>
FiniteElementDomination::Domination
compare_for_face_domination (const FiniteElement<dim,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.
+/**
+ * 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 <int spacedim>
+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 <tt>FE_FaceQ<dim>(degree)</tt>, with <tt>dim</tt> and
+ * <tt>degree</tt> 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
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source,
+ FullMatrix<double> &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
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &source,
+ const unsigned int subface,
+ FullMatrix<double> &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 <tt>flags</tt>.
+ *
+ * For the purpuse of this function, refer to the documentation in
+ * FiniteElement.
+ *
+ * This class assumes that shape functions of this FiniteElement do
+ * <em>not</em> depend on the actual shape of the cells in real
+ * space. Therefore, the effect in this element is as follows: if
+ * <tt>update_values</tt> is set in <tt>flags</tt>, 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 <tt>flags</tt>.
+ *
+ * For the purpuse of this function, refer to the documentation in
+ * FiniteElement.
+ *
+ * This class assumes that shape functions of this FiniteElement do
+ * <em>not</em> depend on the actual shape of the cells in real space.
+ *
+ * The effect in this element is as follows:
+ * <ul>
+ * <li> if <tt>update_gradients</tt> is set, the result will contain
+ * <tt>update_gradients</tt> and <tt>update_covariant_transformation</tt>.
+ * The latter is required to transform the gradient on the unit cell to the
+ * real cell. Remark, that the action required by
+ * <tt>update_covariant_transformation</tt> is actually performed by the
+ * Mapping object used in conjunction with this finite element. <li> if
+ * <tt>update_hessians</tt> is set, the result will contain
+ * <tt>update_hessians</tt> and <tt>update_covariant_transformation</tt>.
+ * 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.
+ *
+ * </ul>
+ */
+ virtual UpdateFlags update_each (const UpdateFlags flags) const;
+
+private:
+ /**
+ * Return vector with dofs per vertex, line, quad, hex.
+ */
+ static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
+};
+
+
/**
* A finite element, which is a Legendre element of complete polynomials on
};
+
+/**
+ * FE_FaceP in 1D, i.e., with degrees of freedom on the element vertices.
+ */
+template <int spacedim>
+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
UpdateFlags
FE_PolyFace<POLY,dim,spacedim>::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;
}
+// ----------------------------- FE_FaceQ<1,spacedim> ------------------------
+
+template <int spacedim>
+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<bool>(1,true),
+ std::vector<ComponentMask> (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 <int spacedim>
+FiniteElement<1,spacedim> *
+FE_FaceQ<1,spacedim>::clone() const
+{
+ return new FE_FaceQ<1,spacedim>(this->degree);
+}
+
+
+
+template <int spacedim>
+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 <int spacedim>
+void
+FE_FaceQ<1,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
+}
+
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> &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 <int spacedim>
+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 <int spacedim>
+std::vector<unsigned int>
+FE_FaceQ<1,spacedim>::get_dpo_vector (const unsigned int)
+{
+ std::vector<unsigned int> dpo(2, 0U);
+ dpo[0] = 1;
+ return dpo;
+}
+
+
+
+template <int spacedim>
+bool
+FE_FaceQ<1,spacedim>::hp_constraints_are_implemented () const
+{
+ return true;
+}
+
+
+
+template <int spacedim>
+FiniteElementDomination::Domination
+FE_FaceQ<1,spacedim>::
+compare_for_face_domination (const FiniteElement<1,spacedim> &fe_other) const
+{
+ return FiniteElementDomination::no_requirements;
+}
+
+
+
+template <int spacedim>
+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; i<this->dofs_per_cell; ++i)
+ constant_modes(0,i) = true;
+ return constant_modes;
+}
+
+
+
+template <int spacedim>
+UpdateFlags
+FE_FaceQ<1,spacedim>::update_once (const UpdateFlags) const
+{
+ return update_default;
+}
+
+
+
+template <int spacedim>
+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 <int spacedim>
+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 <int spacedim>
+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 <int spacedim>
+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 <int spacedim>
+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 <int spacedim>
+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; k<this->dofs_per_cell; ++k)
+ data.shape_values(k,0) = 0.;
+ data.shape_values(foffset,0) = 1;
+ }
+}
+
+
+template <int spacedim>
+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 <int dim, int spacedim>
{}
+
template <int dim, int spacedim>
FiniteElement<dim,spacedim> *
FE_FaceP<dim,spacedim>::clone() const
}
+
template <int dim, int spacedim>
std::string
FE_FaceP<dim,spacedim>::get_name () const
+template <int spacedim>
+FE_FaceP<1,spacedim>::FE_FaceP (const unsigned int degree)
+ :
+ FE_FaceQ<1,spacedim> (degree)
+{}
+
+
+
+template <int spacedim>
+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"
template class FE_PolyFace<TensorProductPolynomials<deal_II_dimension-1> >;
template class FE_PolyFace<PolynomialSpace<deal_II_dimension-1>, deal_II_dimension>;
//template class FE_PolyFace<PolynomialsP<deal_II_dimension>, deal_II_dimension>;
- template class FE_FaceQ<deal_II_dimension>;
- template class FE_FaceP<deal_II_dimension>;
#endif
+ template class FE_FaceQ<deal_II_dimension,deal_II_dimension>;
+ template class FE_FaceP<deal_II_dimension,deal_II_dimension>;
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/logstream.h>
+#include <fstream>
+std::ofstream logfile("output");
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+
+namespace Step51
+{
+
+ using namespace dealii;
+
+ template <int dim>
+ class SolutionBase
+ {
+ protected:
+ static const unsigned int n_source_centers = 3;
+ static const Point<dim> 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 <int dim>
+ const double SolutionBase<dim>::width = 1./5.;
+
+
+ template <int dim>
+ class Solution : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ Solution () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
+
+
+
+ template <int dim>
+ double Solution<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2. * numbers::PI) * this->width);
+ }
+
+
+
+ template <int dim>
+ Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> return_value;
+
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2 * numbers::PI) *
+ this->width);
+ }
+
+
+
+ template <int dim>
+ class SolutionAndGradient : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ SolutionAndGradient () : Function<dim>(dim) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const;
+ };
+
+ template <int dim>
+ void SolutionAndGradient<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ AssertDimension(v.size(), dim+1);
+ Solution<dim> solution;
+ Tensor<1,dim> grad = solution.gradient(p);
+ for (unsigned int d=0; d<dim; ++d)
+ v[d] = -grad[d];
+ v[dim] = solution.value(p);
+ }
+
+
+
+ template <int dim>
+ class ConvectionVelocity : public TensorFunction<1,dim>
+ {
+ public:
+ ConvectionVelocity() : TensorFunction<1,dim>() {}
+
+ virtual Tensor<1,dim> value (const Point<dim> &p) const;
+ };
+
+
+
+ template <int dim>
+ Tensor<1,dim>
+ ConvectionVelocity<dim>::value(const Point<dim> &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 <int dim>
+ class RightHandSide : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ private:
+ const ConvectionVelocity<dim> convection_velocity;
+ };
+
+
+ template <int dim>
+ double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> convection = convection_velocity.value(p);
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2 * numbers::PI)
+ * this->width);
+ }
+
+
+ template <int dim>
+ 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<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data);
+
+ void copy_local_to_global(const PerTaskData &data);
+
+ void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &empty_data);
+
+
+ Triangulation<dim> triangulation;
+
+ FESystem<dim> fe_local;
+ DoFHandler<dim> dof_handler_local;
+ Vector<double> solution_local;
+
+ FE_FaceQ<dim> fe;
+ DoFHandler<dim> dof_handler;
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ FE_DGQ<dim> fe_u_post;
+ DoFHandler<dim> dof_handler_u_post;
+ Vector<double> solution_u_post;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ ConvergenceTable convergence_table;
+ };
+
+
+ template <int dim>
+ HDG<dim>::HDG (const unsigned int degree) :
+ fe_local (FE_DGQ<dim>(degree), dim+1),
+ dof_handler_local (triangulation),
+ fe (degree),
+ dof_handler (triangulation),
+ fe_u_post (degree+1),
+ dof_handler_u_post (triangulation)
+ {}
+
+
+
+ template <int dim>
+ void
+ HDG<dim>::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<dim>::type boundary_functions;
+ Solution<dim> solution_function;
+ boundary_functions[0] = &solution_function;
+ VectorTools::project_boundary_values (dof_handler,
+ boundary_functions,
+ QGauss<dim-1>(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 <int dim>
+ struct HDG<dim>::PerTaskData
+ {
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_vector;
+ std::vector<types::global_dof_index> 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 <int dim>
+ struct HDG<dim>::ScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEFaceValues<dim> fe_face_values_local;
+ FEFaceValues<dim> fe_face_values;
+
+ FullMatrix<double> ll_matrix;
+ FullMatrix<double> lf_matrix;
+ FullMatrix<double> fl_matrix;
+ FullMatrix<double> tmp_matrix;
+ Vector<double> l_rhs;
+ Vector<double> tmp_rhs;
+
+ std::vector<Tensor<1,dim> > q_phi;
+ std::vector<double> q_phi_div;
+ std::vector<double> u_phi;
+ std::vector<Tensor<1,dim> > u_phi_grad;
+ std::vector<double> tr_phi;
+ std::vector<double> trace_values;
+
+ std::vector<std::vector<unsigned int> > fe_local_support_on_face;
+ std::vector<std::vector<unsigned int> > fe_support_on_face;
+
+ ConvectionVelocity<dim> convection_velocity;
+ RightHandSide<dim> right_hand_side;
+ const Solution<dim> exact_solution;
+
+ ScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &quadrature_formula,
+ const QGauss<dim-1> &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<dim>::faces_per_cell),
+ fe_support_on_face(GeometryInfo<dim>::faces_per_cell)
+ {
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe_local.dofs_per_cell; ++i)
+ {
+ if (fe_local.has_support_on_face(i,face))
+ {
+ fe_local_support_on_face[face].push_back(i);
+ }
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ if (fe.has_support_on_face(i,face))
+ {
+ fe_support_on_face[face].push_back(i);
+ }
+ }
+ }
+
+ ScratchData(const ScratchData &sd)
+ :
+ fe_values_local (sd.fe_values_local.get_fe(),
+ sd.fe_values_local.get_quadrature(),
+ sd.fe_values_local.get_update_flags()),
+ fe_face_values_local (sd.fe_face_values_local.get_fe(),
+ sd.fe_face_values_local.get_quadrature(),
+ sd.fe_face_values_local.get_update_flags()),
+ fe_face_values (sd.fe_face_values.get_fe(),
+ sd.fe_face_values.get_quadrature(),
+ sd.fe_face_values.get_update_flags()),
+ ll_matrix (sd.ll_matrix),
+ lf_matrix (sd.lf_matrix),
+ fl_matrix (sd.fl_matrix),
+ tmp_matrix (sd.tmp_matrix),
+ l_rhs (sd.l_rhs),
+ tmp_rhs (sd.tmp_rhs),
+ q_phi (sd.q_phi),
+ q_phi_div (sd.q_phi_div),
+ u_phi (sd.u_phi),
+ u_phi_grad (sd.u_phi_grad),
+ tr_phi (sd.tr_phi),
+ trace_values(sd.trace_values),
+ fe_local_support_on_face(sd.fe_local_support_on_face),
+ fe_support_on_face(sd.fe_support_on_face)
+ {}
+ };
+
+
+
+ template <int dim>
+ struct HDG<dim>::PostProcessScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEValues<dim> fe_values;
+
+ std::vector<double> u_values;
+ std::vector<Tensor<1,dim> > u_gradients;
+ FullMatrix<double> cell_matrix;
+
+ Vector<double> cell_rhs;
+ Vector<double> cell_sol;
+
+ PostProcessScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &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 <int dim>
+ void
+ HDG<dim>::assemble_system (const bool trace_reconstruct)
+ {
+ const QGauss<dim> quadrature_formula(fe.degree+1);
+ const QGauss<dim-1> 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<dim>::assemble_system_one_cell,
+ &HDG<dim>::copy_local_to_global,
+ scratch,
+ task_data);
+ }
+
+
+
+ template <int dim>
+ void
+ HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data)
+ {
+ typename DoFHandler<dim>::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<n_q_points; ++q)
+ {
+ const double rhs_value
+ = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
+ const Tensor<1,dim> 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<loc_dofs_per_cell; ++k)
+ {
+ scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
+ scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
+ scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
+ scratch.u_phi_grad[k] = scratch.fe_values_local[scalar].gradient(k,q);
+ }
+ for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
+ scratch.ll_matrix(i,j) += (
+ scratch.q_phi[i] * scratch.q_phi[j]
+ -
+ scratch.q_phi_div[i] * scratch.u_phi[j]
+ +
+ scratch.u_phi[i] * scratch.q_phi_div[j]
+ -
+ (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
+ ) * JxW;
+ scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
+ }
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::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<n_face_q_points; ++q)
+ {
+ const double JxW = scratch.fe_face_values.JxW(q);
+ const Point<dim> quadrature_point =
+ scratch.fe_face_values.quadrature_point(q);
+ const Point<dim> 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; k<scratch.fe_local_support_on_face[face].size(); ++k)
+ {
+ const unsigned int kk=scratch.fe_local_support_on_face[face][k];
+ scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
+ scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
+ }
+
+ if (!task_data.trace_reconstruct)
+ {
+ for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
+ scratch.tr_phi[k] =
+ scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ scratch.lf_matrix(ii,jj) += (
+ (scratch.q_phi[i] * normal
+ +
+ (convection * normal -
+ tau_stab) * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+
+ scratch.fl_matrix(jj,ii) -= (
+ (scratch.q_phi[i] * normal
+ +
+ tau_stab * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+ }
+
+ for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ task_data.cell_matrix(ii,jj) += (
+ (convection * normal - tau_stab) *
+ scratch.tr_phi[i] * scratch.tr_phi[j]
+ ) * JxW;
+ }
+
+ if (cell->face(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; i<scratch.fe_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_support_on_face[face][i];
+ task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+ }
+ }
+ }
+
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_local_support_on_face[face][j];
+ scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+ }
+
+ if (task_data.trace_reconstruct)
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
+ +
+ scratch.u_phi[i] * (convection * normal - tau_stab)
+ ) * scratch.trace_values[q] * JxW;
+ }
+ }
+ }
+
+ scratch.ll_matrix.gauss_jordan();
+
+ if (task_data.trace_reconstruct == false)
+ {
+ scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
+ scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
+ scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
+ cell->get_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 <int dim>
+ void HDG<dim>::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 <int dim>
+ void HDG<dim>::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 <int dim>
+ void
+ HDG<dim>::postprocess()
+ {
+ {
+ const QGauss<dim> 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<dim>::postprocess_one_cell,
+ std_cxx1x::ref(*this),
+ std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+ std_cxx1x::function<void(const unsigned int &)>(),
+ scratch,
+ 0U);
+ }
+
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+ ComponentSelectFunction<dim> value_select (dim, dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(fe.degree+2),
+ VectorTools::L2_norm,
+ &value_select);
+ const double L2_error = difference_per_cell.l2_norm();
+
+ ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
+ dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(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<dim>(),
+ difference_per_cell,
+ QGauss<dim>(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 <int dim>
+ void
+ HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &)
+ {
+ typename DoFHandler<dim>::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; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += (scratch.fe_values.shape_grad(i,q) *
+ scratch.fe_values.shape_grad(j,q)
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(i,j) = sum;
+ }
+
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(i) = sum;
+ }
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(0,j) = sum;
+ }
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(0) = sum;
+ }
+
+ scratch.cell_matrix.gauss_jordan();
+ scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
+ cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::output_results (const unsigned int cycle)
+ {
+ // not included in test
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::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<dim>::cell_iterator
+ cell = triangulation.begin (),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary())
+ for (unsigned int d=0; d<dim; ++d)
+ if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+ cell->face(face)->set_boundary_indicator (1);
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/logstream.h>
+#include <fstream>
+std::ofstream logfile("output");
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+
+namespace Step51
+{
+
+ using namespace dealii;
+
+ template <int dim>
+ class SolutionBase
+ {
+ protected:
+ static const unsigned int n_source_centers = 3;
+ static const Point<dim> 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 <int dim>
+ const double SolutionBase<dim>::width = 1./5.;
+
+
+ template <int dim>
+ class Solution : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ Solution () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
+
+
+
+ template <int dim>
+ double Solution<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2. * numbers::PI) * this->width);
+ }
+
+
+
+ template <int dim>
+ Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> return_value;
+
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2 * numbers::PI) *
+ this->width);
+ }
+
+
+
+ template <int dim>
+ class SolutionAndGradient : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ SolutionAndGradient () : Function<dim>(dim) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const;
+ };
+
+ template <int dim>
+ void SolutionAndGradient<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ AssertDimension(v.size(), dim+1);
+ Solution<dim> solution;
+ Tensor<1,dim> grad = solution.gradient(p);
+ for (unsigned int d=0; d<dim; ++d)
+ v[d] = -grad[d];
+ v[dim] = solution.value(p);
+ }
+
+
+
+ template <int dim>
+ class ConvectionVelocity : public TensorFunction<1,dim>
+ {
+ public:
+ ConvectionVelocity() : TensorFunction<1,dim>() {}
+
+ virtual Tensor<1,dim> value (const Point<dim> &p) const;
+ };
+
+
+
+ template <int dim>
+ Tensor<1,dim>
+ ConvectionVelocity<dim>::value(const Point<dim> &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 <int dim>
+ class RightHandSide : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ private:
+ const ConvectionVelocity<dim> convection_velocity;
+ };
+
+
+ template <int dim>
+ double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> convection = convection_velocity.value(p);
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> 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<dim>(std::sqrt(2 * numbers::PI)
+ * this->width);
+ }
+
+
+ template <int dim>
+ 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<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data);
+
+ void copy_local_to_global(const PerTaskData &data);
+
+ void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &empty_data);
+
+
+ Triangulation<dim> triangulation;
+
+ FESystem<dim> fe_local;
+ DoFHandler<dim> dof_handler_local;
+ Vector<double> solution_local;
+
+ FE_FaceP<dim> fe;
+ DoFHandler<dim> dof_handler;
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ FE_DGP<dim> fe_u_post;
+ DoFHandler<dim> dof_handler_u_post;
+ Vector<double> solution_u_post;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ ConvergenceTable convergence_table;
+ };
+
+
+ template <int dim>
+ HDG<dim>::HDG (const unsigned int degree) :
+ fe_local (FE_DGP<dim>(degree), dim+1),
+ dof_handler_local (triangulation),
+ fe (degree),
+ dof_handler (triangulation),
+ fe_u_post (degree+1),
+ dof_handler_u_post (triangulation)
+ {}
+
+
+
+ template <int dim>
+ void
+ HDG<dim>::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<dim>::type boundary_functions;
+ Solution<dim> solution_function;
+ boundary_functions[0] = &solution_function;
+ VectorTools::project_boundary_values (dof_handler,
+ boundary_functions,
+ QGauss<dim-1>(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 <int dim>
+ struct HDG<dim>::PerTaskData
+ {
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_vector;
+ std::vector<types::global_dof_index> 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 <int dim>
+ struct HDG<dim>::ScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEFaceValues<dim> fe_face_values_local;
+ FEFaceValues<dim> fe_face_values;
+
+ FullMatrix<double> ll_matrix;
+ FullMatrix<double> lf_matrix;
+ FullMatrix<double> fl_matrix;
+ FullMatrix<double> tmp_matrix;
+ Vector<double> l_rhs;
+ Vector<double> tmp_rhs;
+
+ std::vector<Tensor<1,dim> > q_phi;
+ std::vector<double> q_phi_div;
+ std::vector<double> u_phi;
+ std::vector<Tensor<1,dim> > u_phi_grad;
+ std::vector<double> tr_phi;
+ std::vector<double> trace_values;
+
+ std::vector<std::vector<unsigned int> > fe_local_support_on_face;
+ std::vector<std::vector<unsigned int> > fe_support_on_face;
+
+ ConvectionVelocity<dim> convection_velocity;
+ RightHandSide<dim> right_hand_side;
+ const Solution<dim> exact_solution;
+
+ ScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &quadrature_formula,
+ const QGauss<dim-1> &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<dim>::faces_per_cell),
+ fe_support_on_face(GeometryInfo<dim>::faces_per_cell)
+ {
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe_local.dofs_per_cell; ++i)
+ {
+ if (fe_local.has_support_on_face(i,face))
+ {
+ fe_local_support_on_face[face].push_back(i);
+ }
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ if (fe.has_support_on_face(i,face))
+ {
+ fe_support_on_face[face].push_back(i);
+ }
+ }
+ }
+
+ ScratchData(const ScratchData &sd)
+ :
+ fe_values_local (sd.fe_values_local.get_fe(),
+ sd.fe_values_local.get_quadrature(),
+ sd.fe_values_local.get_update_flags()),
+ fe_face_values_local (sd.fe_face_values_local.get_fe(),
+ sd.fe_face_values_local.get_quadrature(),
+ sd.fe_face_values_local.get_update_flags()),
+ fe_face_values (sd.fe_face_values.get_fe(),
+ sd.fe_face_values.get_quadrature(),
+ sd.fe_face_values.get_update_flags()),
+ ll_matrix (sd.ll_matrix),
+ lf_matrix (sd.lf_matrix),
+ fl_matrix (sd.fl_matrix),
+ tmp_matrix (sd.tmp_matrix),
+ l_rhs (sd.l_rhs),
+ tmp_rhs (sd.tmp_rhs),
+ q_phi (sd.q_phi),
+ q_phi_div (sd.q_phi_div),
+ u_phi (sd.u_phi),
+ u_phi_grad (sd.u_phi_grad),
+ tr_phi (sd.tr_phi),
+ trace_values(sd.trace_values),
+ fe_local_support_on_face(sd.fe_local_support_on_face),
+ fe_support_on_face(sd.fe_support_on_face)
+ {}
+ };
+
+
+
+ template <int dim>
+ struct HDG<dim>::PostProcessScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEValues<dim> fe_values;
+
+ std::vector<double> u_values;
+ std::vector<Tensor<1,dim> > u_gradients;
+ FullMatrix<double> cell_matrix;
+
+ Vector<double> cell_rhs;
+ Vector<double> cell_sol;
+
+ PostProcessScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &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 <int dim>
+ void
+ HDG<dim>::assemble_system (const bool trace_reconstruct)
+ {
+ const QGauss<dim> quadrature_formula(fe.degree+1);
+ const QGauss<dim-1> 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<dim>::assemble_system_one_cell,
+ &HDG<dim>::copy_local_to_global,
+ scratch,
+ task_data);
+ }
+
+
+
+ template <int dim>
+ void
+ HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data)
+ {
+ typename DoFHandler<dim>::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<n_q_points; ++q)
+ {
+ const double rhs_value
+ = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
+ const Tensor<1,dim> 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<loc_dofs_per_cell; ++k)
+ {
+ scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
+ scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
+ scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
+ scratch.u_phi_grad[k] = scratch.fe_values_local[scalar].gradient(k,q);
+ }
+ for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
+ scratch.ll_matrix(i,j) += (
+ scratch.q_phi[i] * scratch.q_phi[j]
+ -
+ scratch.q_phi_div[i] * scratch.u_phi[j]
+ +
+ scratch.u_phi[i] * scratch.q_phi_div[j]
+ -
+ (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
+ ) * JxW;
+ scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
+ }
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::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<n_face_q_points; ++q)
+ {
+ const double JxW = scratch.fe_face_values.JxW(q);
+ const Point<dim> quadrature_point =
+ scratch.fe_face_values.quadrature_point(q);
+ const Point<dim> 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; k<scratch.fe_local_support_on_face[face].size(); ++k)
+ {
+ const unsigned int kk=scratch.fe_local_support_on_face[face][k];
+ scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
+ scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
+ }
+
+ if (!task_data.trace_reconstruct)
+ {
+ for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
+ scratch.tr_phi[k] =
+ scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ scratch.lf_matrix(ii,jj) += (
+ (scratch.q_phi[i] * normal
+ +
+ (convection * normal -
+ tau_stab) * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+
+ scratch.fl_matrix(jj,ii) -= (
+ (scratch.q_phi[i] * normal
+ +
+ tau_stab * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+ }
+
+ for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ task_data.cell_matrix(ii,jj) += (
+ (convection * normal - tau_stab) *
+ scratch.tr_phi[i] * scratch.tr_phi[j]
+ ) * JxW;
+ }
+
+ if (cell->face(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; i<scratch.fe_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_support_on_face[face][i];
+ task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+ }
+ }
+ }
+
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_local_support_on_face[face][j];
+ scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+ }
+
+ if (task_data.trace_reconstruct)
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
+ +
+ scratch.u_phi[i] * (convection * normal - tau_stab)
+ ) * scratch.trace_values[q] * JxW;
+ }
+ }
+ }
+
+ scratch.ll_matrix.gauss_jordan();
+
+ if (task_data.trace_reconstruct == false)
+ {
+ scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
+ scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
+ scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
+ cell->get_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 <int dim>
+ void HDG<dim>::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 <int dim>
+ void HDG<dim>::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 <int dim>
+ void
+ HDG<dim>::postprocess()
+ {
+ {
+ const QGauss<dim> 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<dim>::postprocess_one_cell,
+ std_cxx1x::ref(*this),
+ std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+ std_cxx1x::function<void(const unsigned int &)>(),
+ scratch,
+ 0U);
+ }
+
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+ ComponentSelectFunction<dim> value_select (dim, dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(fe.degree+2),
+ VectorTools::L2_norm,
+ &value_select);
+ const double L2_error = difference_per_cell.l2_norm();
+
+ ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
+ dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(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<dim>(),
+ difference_per_cell,
+ QGauss<dim>(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 <int dim>
+ void
+ HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &)
+ {
+ typename DoFHandler<dim>::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; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += (scratch.fe_values.shape_grad(i,q) *
+ scratch.fe_values.shape_grad(j,q)
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(i,j) = sum;
+ }
+
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(i) = sum;
+ }
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(0,j) = sum;
+ }
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(0) = sum;
+ }
+
+ scratch.cell_matrix.gauss_jordan();
+ scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
+ cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::output_results (const unsigned int cycle)
+ {
+ // not included in test
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::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<dim>::cell_iterator
+ cell = triangulation.begin (),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary())
+ for (unsigned int d=0; d<dim; ++d)
+ if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+ cell->face(face)->set_boundary_indicator (1);
+ }
+
+
+
+ template <int dim>
+ void HDG<dim>::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;
+}
--- /dev/null
+
+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
template <int dim>
void test_2d_3d (std::vector<FiniteElement<dim> *> &fe_datas)
{
- // Face Q elements
- fe_datas.push_back(new FE_FaceQ<dim> (0));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
- fe_datas.push_back(new FE_FaceQ<dim> (1));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
- fe_datas.push_back(new FE_FaceQ<dim> (3));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
- // Face P elements
- fe_datas.push_back(new FE_FaceP<dim> (0));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
- fe_datas.push_back(new FE_FaceP<dim> (1));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
- fe_datas.push_back(new FE_FaceP<dim> (3));
- deallog << (*fe_datas.rbegin())->get_name() << std::endl;
// Vector DG elements
fe_datas.push_back(
new FE_DGRaviartThomas<dim>(0));
FE_Q<dim> (2), 1));
deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ // Face Q elements
+ fe_datas.push_back(new FE_FaceQ<dim> (0));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ fe_datas.push_back(new FE_FaceQ<dim> (1));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ fe_datas.push_back(new FE_FaceQ<dim> (3));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ // Face P elements
+ fe_datas.push_back(new FE_FaceP<dim> (0));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ fe_datas.push_back(new FE_FaceP<dim> (1));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+ fe_datas.push_back(new FE_FaceP<dim> (3));
+ deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+
// Check vector elements in 2d and higher only
test_2d_3d (fe_datas);
projects / dealii-svn.git / commitdiff