--- /dev/null
+#ifndef __deal2__tensor_product_polynomials_const_h
+#define __deal2__tensor_product_polynomials_const_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @addtogroup Polynomials
+ * @{
+ */
+
+/**
+ * Tensor product of given polynomials.
+ *
+ * Given a vector of <i>n</i> one-dimensional polynomials
+ * <i>P<sub>1</sub></i> to <i>P<sub>n</sub></i>, this class generates
+ * <i>n<sup>dim</sup></i> polynomials of the form
+ * <i>Q<sub>ijk</sub>(x,y,z) =
+ * P<sub>i</sub>(x)P<sub>j</sub>(y)P<sub>k</sub>(z)</i>and a locally
+ * constant function. If the base polynomials are mutually orthogonal
+ * on the interval [-1,1] or [0,1], then the tensor product
+ * polynomials are orthogonal on [-1,1]<sup>dim</sup> or [0,1]
+ * <sup>dim</sup>, respectively.
+ *
+ * Indexing is as follows: the order of dim-dimensional polynomials is
+ * x-coordinates running fastest, then y-coordinate, etc. The first
+ * few polynomials are thus <i>P<sub>1</sub>(x)P<sub>1</sub>(y),
+ * P<sub>2</sub>(x)P<sub>1</sub>(y), P<sub>3</sub>(x)P<sub>1</sub>(y),
+ * ..., P<sub>1</sub>(x)P<sub>2</sub>(y),
+ * P<sub>2</sub>(x)P<sub>2</sub>(y), P<sub>3</sub>(x)P<sub>2</sub>(y),
+ * ...</i> and likewise in 3d. The locally constant function has the last index
+ *
+ * The output_indices() function prints the ordering of the
+ * dim-dimensional polynomials, i.e. for each polynomial in the
+ * polynomial space it gives the indices i,j,k of the one-dimensional
+ * polynomials in x,y and z direction. The ordering of the
+ * dim-dimensional polynomials can be changed by using the
+ * set_numbering() function.
+ *
+ * @author Ralf Hartmann, 2000, 2004, Guido Kanschat, 2000, Wolfgang Bangerth 2003
+ */
+template <int dim>
+class TensorProductPolynomialsConst
+{
+public:
+ /**
+ * Access to the dimension of
+ * this object, for checking and
+ * automatic setting of dimension
+ * in other classes.
+ */
+ static const unsigned int dimension = dim;
+
+ /**
+ * Constructor. <tt>pols</tt> is
+ * a vector of objects that
+ * should be derived or otherwise
+ * convertible to one-dimensional
+ * polynomial objects. It will be
+ * copied element by element into
+ * a private variable.
+ */
+ template <class Pol>
+ TensorProductPolynomialsConst (const std::vector<Pol> &pols);
+
+ /**
+ * Prints the list of the indices
+ * to <tt>out</tt>.
+ */
+ void output_indices(std::ostream &out) const;
+
+ /**
+ * Sets the ordering of the
+ * polynomials. Requires
+ * <tt>renumber.size()==n()</tt>.
+ * Stores a copy of
+ * <tt>renumber</tt>.
+ */
+ void set_numbering(const std::vector<unsigned int> &renumber);
+
+ /**
+ * Gives read access to the
+ * renumber vector.
+ */
+ const std::vector<unsigned int> &get_numbering() const;
+
+ /**
+ * Gives read access to the
+ * inverse renumber vector.
+ */
+ const std::vector<unsigned int> &get_numbering_inverse() const;
+
+ /**
+ * Computes the value and the
+ * first and second derivatives
+ * of each tensor product
+ * polynomial at <tt>unit_point</tt>.
+ *
+ * The size of the vectors must
+ * either be equal 0 or equal
+ * n(). In the first case, the
+ * function will not compute
+ * these values.
+ *
+ * If you need values or
+ * derivatives of all tensor
+ * product polynomials then use
+ * this function, rather than
+ * using any of the
+ * compute_value(),
+ * compute_grad() or
+ * compute_grad_grad()
+ * functions, see below, in a
+ * loop over all tensor product
+ * polynomials.
+ */
+ void compute (const Point<dim> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &grad_grads) const;
+
+ /**
+ * Computes the value of the
+ * <tt>i</tt>th tensor product
+ * polynomial at
+ * <tt>unit_point</tt>. Here <tt>i</tt> is
+ * given in tensor product
+ * numbering.
+ *
+ * Note, that using this function
+ * within a loop over all tensor
+ * product polynomials is not
+ * efficient, because then each
+ * point value of the underlying
+ * (one-dimensional) polynomials
+ * is (unnecessarily) computed
+ * several times. Instead use
+ * the compute() function with
+ * <tt>values.size()==</tt>n()
+ * to get the point values of all
+ * tensor polynomials all at once
+ * and in a much more efficient
+ * way.
+ */
+ double compute_value (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Computes the grad of the
+ * <tt>i</tt>th tensor product
+ * polynomial at
+ * <tt>unit_point</tt>. Here <tt>i</tt> is
+ * given in tensor product
+ * numbering.
+ *
+ * Note, that using this function
+ * within a loop over all tensor
+ * product polynomials is not
+ * efficient, because then each
+ * derivative value of the
+ * underlying (one-dimensional)
+ * polynomials is (unnecessarily)
+ * computed several times.
+ * Instead use the compute()
+ * function, see above, with
+ * <tt>grads.size()==</tt>n()
+ * to get the point value of all
+ * tensor polynomials all at once
+ * and in a much more efficient
+ * way.
+ */
+ Tensor<1,dim> compute_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Computes the second
+ * derivative (grad_grad) of the
+ * <tt>i</tt>th tensor product
+ * polynomial at
+ * <tt>unit_point</tt>. Here <tt>i</tt> is
+ * given in tensor product
+ * numbering.
+ *
+ * Note, that using this function
+ * within a loop over all tensor
+ * product polynomials is not
+ * efficient, because then each
+ * derivative value of the
+ * underlying (one-dimensional)
+ * polynomials is (unnecessarily)
+ * computed several times.
+ * Instead use the compute()
+ * function, see above, with
+ * <tt>grad_grads.size()==</tt>n()
+ * to get the point value of all
+ * tensor polynomials all at once
+ * and in a much more efficient
+ * way.
+ */
+ Tensor<2,dim> compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Returns the number of tensor
+ * product polynomials plus the constant
+ * function. For <i>n</i> 1d polynomials
+ * this is <i>n<sup>dim</sup>+1</i>.
+ */
+ unsigned int n () const;
+
+
+private:
+ /**
+ * Copy of the vector <tt>pols</tt> of
+ * polynomials given to the
+ * constructor.
+ */
+ std::vector<Polynomials::Polynomial<double> > polynomials;
+
+ /**
+ * Number of tensor product
+ * polynomials. See n().
+ */
+ unsigned int n_tensor_pols;
+
+ /**
+ * Index map for reordering the
+ * polynomials.
+ */
+ std::vector<unsigned int> index_map;
+
+ /**
+ * Index map for reordering the
+ * polynomials.
+ */
+ std::vector<unsigned int> index_map_inverse;
+
+ /**
+ * Each tensor product polynomial
+ * <i>i</i> is a product of
+ * one-dimensional polynomials in
+ * each space direction. Compute
+ * the indices of these
+ * one-dimensional polynomials
+ * for each space direction,
+ * given the index <i>i</i>.
+ */
+ // fix to avoid compiler warnings about zero
+ // length arrays
+ void compute_index (const unsigned int i,
+ unsigned int (&indices)[(dim>0?dim:1)]) const;
+
+ /**
+ * Computes
+ * <i>x<sup>dim</sup></i> for
+ * unsigned int <i>x</i>. Used in
+ * the constructor.
+ */
+ static
+ unsigned int x_to_the_dim (const unsigned int x);
+};
+
+#ifndef DOXYGEN
+
+template <int dim>
+inline
+const std::vector<unsigned int> &
+TensorProductPolynomialsConst<dim>::get_numbering() const
+{
+ return index_map;
+}
+
+
+template <int dim>
+inline
+const std::vector<unsigned int> &
+TensorProductPolynomialsConst<dim>::get_numbering_inverse() const
+{
+ return index_map_inverse;
+}
+
+
+#endif // DOXYGEN
+
+#ifndef DOXYGEN
+
+/* -------------- declaration of explicit specializations --- */
+
+template <>
+void
+TensorProductPolynomialsConst<1>::compute_index(const unsigned int n,
+ unsigned int (&index)[1]) const;
+template <>
+void
+TensorProductPolynomialsConst<2>::compute_index(const unsigned int n,
+ unsigned int (&index)[2]) const;
+template <>
+void
+TensorProductPolynomialsConst<3>::compute_index(const unsigned int n,
+ unsigned int (&index)[3]) const;
+
+
+/* ---------------- template and inline functions ---------- */
+
+template <int dim>
+inline
+unsigned int
+TensorProductPolynomialsConst<dim>::
+x_to_the_dim (const unsigned int x)
+{
+ unsigned int y = 1;
+ for (int d=0; d<dim; ++d)
+ y *= x;
+ return y;
+}
+
+
+
+template <int dim>
+template <class Pol>
+TensorProductPolynomialsConst<dim>::
+TensorProductPolynomialsConst(const std::vector<Pol> &pols)
+ :
+ polynomials (pols.begin(), pols.end()),
+ n_tensor_pols(x_to_the_dim(pols.size())+1),
+ index_map(n_tensor_pols),
+ index_map_inverse(n_tensor_pols)
+{
+ // per default set this index map
+ // to identity. This map can be
+ // changed by the user through the
+ // set_numbering() function
+ for (unsigned int i=0; i<n_tensor_pols; ++i)
+ {
+ index_map[i]=i;
+ index_map_inverse[i]=i;
+ }
+}
+
+#endif // DOXYGEN
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+#ifndef __deal2__fe_q_dg0_h
+#define __deal2__fe_q_dg0_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/tensor_product_polynomials_const.h>
+#include <deal.II/fe/fe_poly.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*!@addtogroup fe */
+/*@{*/
+
+/**
+ * Implementation of a scalar Lagrange finite element @p Qp+DG0 that yields the
+ * finite element space of continuous, piecewise polynomials of degree @p p in
+ * each coordinate direction plus the space of locally constant functions.
+ * This class is realized using tensor product polynomials based on equidistant
+ * or given support points.
+ *
+ * The standard constructor of this class takes the degree @p p of this finite
+ * element. Alternatively, it can take a quadrature formula @p points defining
+ * the support points of the Lagrange interpolation in one coordinate direction.
+ *
+ * For more information about the <tt>spacedim</tt> template parameter check
+ * the documentation of FiniteElement or the one of Triangulation.
+ *
+ * <h3>Implementation</h3>
+ *
+ * The constructor creates a TensorProductPolynomials object that includes the
+ * tensor product of @p LagrangeEquidistant polynomials of degree @p p plus the
+ * locally constant function. This @p TensorProductPolynomialsConst object
+ * provides all values and derivatives of the shape functions. In case a
+ * quadrature rule is given, the constructor creates a
+ * TensorProductPolynomialsConst object that includes the tensor product of @p
+ * Lagrange polynomials with the * support points from @p points and a locally
+ * constant function.
+ *
+ * <h3>Numbering of the degrees of freedom (DoFs)</h3>
+ *
+ * The original ordering of the shape functions represented by the
+ * TensorProductPolynomialsConst is a tensor product
+ * numbering. However, the shape functions on a cell are renumbered
+ * beginning with the shape functions whose support points are at the
+ * vertices, then on the line, on the quads, and finally (for 3d) on
+ * the hexes. Finally there is a support point for the discontinuous shape
+ * function in the middle of the cell. To be explicit, these numberings are
+ * listed in the following:
+ *
+ * <h4>Q1 elements</h4>
+ * <ul>
+ * <li> 1D case:
+ * @verbatim
+ * 0---2---1
+ * @endverbatim
+ *
+ * <li> 2D case:
+ * @verbatim
+ * 2-------3
+ * | |
+ * | 5 |
+ * | |
+ * 0-------1
+ * @endverbatim
+ *
+ * <li> 3D case:
+ * @verbatim
+ * 6-------7 6-------7
+ * /| | / /|
+ * / | | / / |
+ * / | | / / |
+ * 4 | 8 | 4-------5 |
+ * | 2-------3 | | 3
+ * | / / | | /
+ * | / / | | /
+ * |/ / | |/
+ * 0-------1 0-------1
+ * @endverbatim
+ *
+ * The respective coordinate values of the support points of the degrees
+ * of freedom are as follows:
+ * <ul>
+ * <li> Index 0: <tt>[ 0, 0, 0]</tt>;
+ * <li> Index 1: <tt>[ 1, 0, 0]</tt>;
+ * <li> Index 2: <tt>[ 0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[ 1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[ 0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[ 1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[ 0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[ 1, 1, 1]</tt>;
+ * <li> Index 8: <tt>[1/2, 1/2, 1/2]</tt>;
+ * </ul>
+ * </ul>
+ * <h4>Q2 elements</h4>
+ * <ul>
+ * <li> 1D case:
+ * @verbatim
+ * 0---2---1
+ * @endverbatim
+ * Index 3 has the same coordinates as index 2
+ *
+ * <li> 2D case:
+ * @verbatim
+ * 2---7---3
+ * | |
+ * 4 8 5
+ * | |
+ * 0---6---1
+ * @endverbatim
+ * Index 9 has the same coordinates as index 2
+ *
+ * <li> 3D case:
+ * @verbatim
+ * 6--15---7 6--15---7
+ * /| | / /|
+ * 12 | 19 12 1319
+ * / 18 | / / |
+ * 4 | | 4---14--5 |
+ * | 2---11--3 | | 3
+ * | / / | 17 /
+ * 16 8 9 16 | 9
+ * |/ / | |/
+ * 0---10--1 0---8---1
+ *
+ * *-------* *-------*
+ * /| | / /|
+ * / | 23 | / 25 / |
+ * / | | / / |
+ * * | | *-------* |
+ * |20 *-------* | |21 *
+ * | / / | 22 | /
+ * | / 24 / | | /
+ * |/ / | |/
+ * *-------* *-------*
+ * @endverbatim
+ * The center vertices have number 26 and 27.
+ *
+ * The respective coordinate values of the support points of the degrees
+ * of freedom are as follows:
+ * <ul>
+ * <li> Index 0: <tt>[0, 0, 0]</tt>;
+ * <li> Index 1: <tt>[1, 0, 0]</tt>;
+ * <li> Index 2: <tt>[0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[1, 1, 1]</tt>;
+ * <li> Index 8: <tt>[0, 1/2, 0]</tt>;
+ * <li> Index 9: <tt>[1, 1/2, 0]</tt>;
+ * <li> Index 10: <tt>[1/2, 0, 0]</tt>;
+ * <li> Index 11: <tt>[1/2, 1, 0]</tt>;
+ * <li> Index 12: <tt>[0, 1/2, 1]</tt>;
+ * <li> Index 13: <tt>[1, 1/2, 1]</tt>;
+ * <li> Index 14: <tt>[1/2, 0, 1]</tt>;
+ * <li> Index 15: <tt>[1/2, 1, 1]</tt>;
+ * <li> Index 16: <tt>[0, 0, 1/2]</tt>;
+ * <li> Index 17: <tt>[1, 0, 1/2]</tt>;
+ * <li> Index 18: <tt>[0, 1, 1/2]</tt>;
+ * <li> Index 19: <tt>[1, 1, 1/2]</tt>;
+ * <li> Index 20: <tt>[0, 1/2, 1/2]</tt>;
+ * <li> Index 21: <tt>[1, 1/2, 1/2]</tt>;
+ * <li> Index 22: <tt>[1/2, 0, 1/2]</tt>;
+ * <li> Index 23: <tt>[1/2, 1, 1/2]</tt>;
+ * <li> Index 24: <tt>[1/2, 1/2, 0]</tt>;
+ * <li> Index 25: <tt>[1/2, 1/2, 1]</tt>;
+ * <li> Index 26: <tt>[1/2, 1/2, 1/2]</tt>;
+ * <li> Index 27: <tt>[1/2, 1/2, 1/2]</tt>;
+ * </ul>
+ * </ul>
+ * <h4>Q3 elements</h4>
+ * <ul>
+ * <li> 1D case:
+ * @verbatim
+ * 0--2-4-3--1
+ * @endverbatim
+ *
+ * <li> 2D case:
+ * @verbatim
+ * 2--10-11-3
+ * | |
+ * 5 14 15 7
+ * | 16 |
+ * 4 12 13 6
+ * | |
+ * 0--8--9--1
+ * @endverbatim
+ * </ul>
+ * <h4>Q4 elements</h4>
+ * <ul>
+ * <li> 1D case:
+ * @verbatim
+ * 0--2--3--4--1
+ * @endverbatim
+ * Index 5 has the same coordinates as index 3
+ *
+ * <li> 2D case:
+ * @verbatim
+ * 2--13-14-15-3
+ * | |
+ * 6 22 23 24 9
+ * | |
+ * 5 19 20 21 8
+ * | |
+ * 4 16 17 18 7
+ * | |
+ * 0--10-11-12-1
+ * @endverbatim
+ * Index 21 has the same coordinates as index 20
+ * </ul>
+ *
+ */
+template <int dim, int spacedim=dim>
+class FE_Q_DG0 : public FE_Poly<TensorProductPolynomialsConst<dim>,dim,spacedim>
+{
+public:
+ /**
+ * Constructor for tensor product
+ * polynomials of degree @p p plus locally
+ * constant functions.
+ */
+ FE_Q_DG0 (const unsigned int p);
+
+ /**
+ * Constructor for tensor product
+ * polynomials with support
+ * points @p points plus locally constant
+ * functions based on a one-dimensional
+ * quadrature formula. The degree of the
+ * finite element is
+ * <tt>points.size()-1</tt>.
+ * Note that the first point has
+ * to be 0 and the last one 1.
+ */
+ FE_Q_DG0 (const Quadrature<1> &points);
+
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * <tt>FE_Q_DG0<dim>(degree)</tt>, with
+ * @p dim and @p degree
+ * replaced by appropriate
+ * values.
+ */
+ virtual std::string get_name () const;
+
+ /**
+ * Interpolate a set of scalar
+ * values, computed in the
+ * generalized support points.
+ */
+ virtual void interpolate(std::vector<double> &local_dofs,
+ const std::vector<double> &values) const;
+
+ /**
+ * Interpolate a set of vector
+ * values, computed in the
+ * generalized support points.
+ *
+ * Since a finite element often
+ * only interpolates part of a
+ * vector, <tt>offset</tt> is
+ * used to determine the first
+ * component of the vector to be
+ * interpolated. Maybe consider
+ * changing your data structures
+ * to use the next function.
+ */
+ virtual void interpolate(std::vector<double> &local_dofs,
+ const std::vector<Vector<double> > &values,
+ unsigned int offset = 0) const;
+
+ /**
+ * Interpolate a set of vector
+ * values, computed in the
+ * generalized support points.
+ */
+ virtual void interpolate(
+ std::vector<double> &local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > > &values) const;
+
+ /**
+ * Return the matrix interpolating from the
+ * given finite element to the present one.
+ * The size of the matrix is then @p
+ * dofs_per_cell times
+ * <tt>source.dofs_per_cell</tt>.
+ *
+ * These matrices are only available if the
+ * source element is also a @p FE_Q_DG0
+ * element. Otherwise, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented
+ * is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElement<dim,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>.
+ *
+ * Derived elements will have to
+ * implement this function. They
+ * may only provide interpolation
+ * matrices for certain source
+ * finite elements, for example
+ * those from the same family. If
+ * they don't implement
+ * interpolation from a given
+ * element, then they must throw
+ * an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+ */
+ virtual void
+ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Return the matrix
+ * interpolating from a face of
+ * of one element to the 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>.
+ *
+ * Derived elements will have to
+ * implement this function. They
+ * may only provide interpolation
+ * matrices for certain source
+ * finite elements, for example
+ * those from the same family. If
+ * they don't implement
+ * interpolation from a given
+ * element, then they must throw
+ * an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+ */
+ virtual void
+ get_subface_interpolation_matrix (const FiniteElement<dim,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;
+
+ /**
+ * @name Functions to support hp
+ * @{
+ */
+
+ /**
+ * Return whether this element
+ * implements its hanging node
+ * constraints in the new way,
+ * which has to be used to make
+ * elements "hp compatible".
+ *
+ * For the FE_Q_DG0 class the result is
+ * always true (independent of the degree
+ * of the element), as it implements the
+ * complete set of functions necessary
+ * for hp capability.
+ */
+ virtual bool hp_constraints_are_implemented () const;
+
+ /**
+ * If, on a vertex, several
+ * finite elements are active,
+ * the hp code first assigns the
+ * degrees of freedom of each of
+ * these FEs different global
+ * indices. It then calls this
+ * function to find out which of
+ * them should get identical
+ * values, and consequently can
+ * receive the same global DoF
+ * index. This function therefore
+ * returns a list of identities
+ * between DoFs of the present
+ * finite element object with the
+ * DoFs of @p fe_other, which is
+ * a reference to a finite
+ * element object representing
+ * one of the other finite
+ * elements active on this
+ * particular vertex. The
+ * function computes which of the
+ * degrees of freedom of the two
+ * finite element objects are
+ * equivalent, and returns a list
+ * of pairs of global dof indices
+ * in @p identities. The first
+ * index of each pair denotes one
+ * of the vertex dofs of the
+ * present element, whereas the
+ * second is the corresponding
+ * index of the other finite
+ * element.
+ */
+ virtual
+ std::vector<std::pair<unsigned int, unsigned int> >
+ hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+ /**
+ * Same as
+ * hp_vertex_dof_indices(),
+ * except that the function
+ * treats degrees of freedom on
+ * lines.
+ */
+ virtual
+ std::vector<std::pair<unsigned int, unsigned int> >
+ hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+ /**
+ * Same as
+ * hp_vertex_dof_indices(),
+ * except that the function
+ * treats degrees of freedom on
+ * quads.
+ */
+ virtual
+ std::vector<std::pair<unsigned int, unsigned int> >
+ hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
+
+ /**
+ * Return whether this element dominates
+ * the one given as argument when they
+ * meet at a common face,
+ * whether it is the other way around,
+ * whether neither dominates, or if
+ * either could dominate.
+ *
+ * For a definition of domination, see
+ * FiniteElementBase::Domination and in
+ * particular the @ref hp_paper "hp paper".
+ */
+ virtual
+ FiniteElementDomination::Domination
+ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+ //@}
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ *
+ * This function is made virtual,
+ * since finite element objects
+ * are usually accessed through
+ * pointers to their base class,
+ * rather than the class itself.
+ */
+ virtual std::size_t memory_consumption () const;
+
+protected:
+ /**
+ * @p clone function instead of
+ * a copy constructor.
+ *
+ * This function is needed by the
+ * constructors of @p FESystem.
+ */
+ virtual FiniteElement<dim,spacedim> *clone() const;
+
+private:
+ /**
+ * Returns the restriction_is_additive flags.
+ * Only the last component is true.
+ */
+ static std::vector<bool> get_riaf_vector(const unsigned int degree);
+
+ /**
+ * Only for internal use. Its
+ * full name is
+ * @p get_dofs_per_object_vector
+ * function and it creates the
+ * @p dofs_per_object vector that is
+ * needed within the constructor to
+ * be passed to the constructor of
+ * @p FiniteElementData.
+ */
+ static std::vector<unsigned int> get_dpo_vector(const unsigned int degree);
+
+ /**
+ * Only for internal use. Its
+ * full name is
+ * @p get_dofs_per_object_vector_q
+ * function and it creates the
+ * @p dofs_per_object vector that is
+ * needed within the constructor to
+ * be passed to the constructor of
+ * @p FiniteElementData for FE_Q_DG0
+ * objects.
+ */
+ static std::vector<unsigned int> get_dpo_vector_q(const unsigned int degree);
+
+ /**
+ * This is an analogon to the
+ * FETools::lexicographic_to_hierarchic_numbering
+ * function, but working on
+ * faces. Called from the
+ * constructor.
+ */
+ static
+ std::vector<unsigned int>
+ face_lexicographic_to_hierarchic_numbering (const unsigned int degree);
+
+ /**
+ * Initialize the hanging node
+ * constraints matrices. Called
+ * from the constructor.
+ */
+ void initialize_constraints ();
+
+ /**
+ * Initialize the
+ * @p unit_support_points field
+ * of the FiniteElement
+ * class. Called from the
+ * constructor.
+ */
+ void initialize_unit_support_points ();
+
+ /**
+ * Initialize the @p unit_support_points
+ * field of the FiniteElement
+ * class. Called from the constructor in
+ * case the finite element is based on
+ * quadrature points.
+ */
+ void initialize_unit_support_points (const Quadrature<1> &points);
+
+ /**
+ * Initialize the
+ * @p unit_face_support_points field
+ * of the FiniteElement
+ * class. Called from the
+ * constructor.
+ */
+ void initialize_unit_face_support_points ();
+
+ /**
+ * Initialize the @p
+ * unit_face_support_points field of the
+ * FiniteElement class. Called from the
+ * constructor in case the finite element
+ * is based on quadrature points.
+ */
+ void initialize_unit_face_support_points (const Quadrature<1> &points);
+
+ /**
+ * Initialize the
+ * @p adjust_quad_dof_index_for_face_orientation_table field
+ * of the FiniteElement
+ * class. Called from the
+ * constructor.
+ */
+ void initialize_quad_dof_index_permutation ();
+
+ /**
+ * Mapping from hierarchic to
+ * lexicographic numbering on
+ * first face. Hierarchic is the
+ * numbering of the shape
+ * functions.
+ */
+ const std::vector<unsigned int> face_index_map;
+
+
+ /**
+ * Forward declaration of a class
+ * into which we put significant
+ * parts of the implementation.
+ *
+ * See the .cc file for more
+ * information.
+ */
+ struct Implementation;
+
+ /**
+ * Allow access from other
+ * dimensions. We need this since
+ * we want to call the functions
+ * @p get_dpo_vector and
+ * @p lexicographic_to_hierarchic_numbering
+ * for the faces of the finite
+ * element of dimension dim+1.
+ */
+ template <int, int> friend class FE_Q_DG0;
+
+ friend struct FE_Q_DG0<dim,spacedim>::Implementation;
+};
+
+
+
+/*@}*/
+
+/* -------------- declaration of explicit specializations ------------- */
+
+template <>
+void FE_Q_DG0<1>::initialize_unit_face_support_points ();
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
+
+
+template <>
+void FE_Q_DG0<1,2>::initialize_unit_face_support_points ();
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1,2>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
+
+template <>
+void FE_Q_DG0<1,3>::initialize_unit_face_support_points ();
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1,3>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+#include <deal.II/base/tensor_product_polynomials_const.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/table.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/* ------------------- TensorProductPolynomialsConst -------------- */
+
+template <>
+inline
+void
+TensorProductPolynomialsConst<1>::
+compute_index (const unsigned int i,
+ unsigned int (&indices)[1]) const
+{
+ Assert (i<polynomials.size(), ExcInternalError());
+ indices[0] = index_map[i];
+}
+
+template <>
+inline
+void
+TensorProductPolynomialsConst<2>::
+compute_index (const unsigned int i,
+ unsigned int (&indices)[2]) const
+{
+ const unsigned int n_pols = polynomials.size();
+ Assert (i<n_pols*n_pols, ExcInternalError());
+ const unsigned int n=index_map[i];
+
+ indices[0] = n % n_pols;
+ indices[1] = n / n_pols;
+}
+
+template <>
+inline
+void
+TensorProductPolynomialsConst<3>::
+compute_index (const unsigned int i,
+ unsigned int (&indices)[3]) const
+{
+ const unsigned int n_pols = polynomials.size();
+ Assert (i<n_pols*n_pols*n_pols, ExcInternalError());
+ const unsigned int n=index_map[i];
+
+ indices[0] = n % n_pols;
+ indices[1] = (n/n_pols) % n_pols;
+ indices[2] = n / (n_pols*n_pols);
+}
+
+
+template <int dim>
+void
+TensorProductPolynomialsConst<dim>::output_indices(std::ostream &out) const
+{
+ unsigned int ix[dim];
+ //ouput without constant locally constant function
+ for (unsigned int i=0; i<n_tensor_pols-1; ++i)
+ {
+ compute_index(i,ix);
+ out << i << "\t";
+ for (unsigned int d=0; d<dim; ++d)
+ out << ix[d] << " ";
+ out << std::endl;
+ }
+}
+
+
+
+template <int dim>
+void
+TensorProductPolynomialsConst<dim>::set_numbering(
+ const std::vector<unsigned int> &renumber)
+{
+ Assert(renumber.size()==index_map.size(),
+ ExcDimensionMismatch(renumber.size(), index_map.size()));
+
+ index_map=renumber;
+ for (unsigned int i=0; i<index_map.size(); ++i)
+ index_map_inverse[index_map[i]]=i;
+}
+
+template <int dim>
+double
+TensorProductPolynomialsConst<dim>::compute_value (const unsigned int i,
+ const Point<dim> &p) const
+{
+ unsigned int max_indices=pow(polynomials.size(),dim);
+
+ Assert (i<=max_indices, ExcInternalError());
+ if (i<max_indices)
+ {
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ double value=1.;
+ for (unsigned int d=0; d<dim; ++d)
+ value *= polynomials[indices[d]].value(p(d));
+
+ return value;
+ }
+ else
+ return 1.; //return 1 for q0 node
+}
+
+template <>
+double
+TensorProductPolynomialsConst<0>::compute_value (const unsigned int ,
+ const Point<0> &) const
+{
+ Assert (false, ExcNotImplemented());
+ return 0.;
+}
+
+template <int dim>
+Tensor<1,dim>
+TensorProductPolynomialsConst<dim>::compute_grad (const unsigned int i,
+ const Point<dim> &p) const
+{
+ unsigned int max_indices=pow(polynomials.size(),dim);
+
+ Assert (i<=max_indices, ExcInternalError());
+ Tensor<1,dim> grad;
+ if (i<max_indices)
+ {
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ // compute values and
+ // uni-directional derivatives at
+ // the given point in each
+ // co-ordinate direction
+ double v [dim][2];
+ {
+ std::vector<double> tmp (2);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ polynomials[indices[d]].value (p(d), tmp);
+ v[d][0] = tmp[0];
+ v[d][1] = tmp[1];
+ }
+ }
+
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ grad[d] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grad[d] *= v[x][d==x];
+ }
+ }
+ return grad; //return 0 for q0 node
+}
+
+template <int dim>
+Tensor<2,dim>
+TensorProductPolynomialsConst<dim>::compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const
+{
+ unsigned int max_indices=pow(polynomials.size(),dim);
+
+ Assert (i<=max_indices, ExcInternalError());
+ Tensor<2,dim> grad_grad;
+ if (i<max_indices)
+ {
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ double v [dim][3];
+ {
+ std::vector<double> tmp (3);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ polynomials[indices[d]].value (p(d), tmp);
+ v[d][0] = tmp[0];
+ v[d][1] = tmp[1];
+ v[d][2] = tmp[2];
+ }
+ }
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ grad_grad[d1][d2] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x || d2==x)
+ {
+ if (d1==d2)
+ derivative=2;
+ else
+ derivative=1;
+ }
+ grad_grad[d1][d2] *= v[x][derivative];
+ }
+ }
+ }
+ return grad_grad; //return 0 for q0 node
+}
+
+template <int dim>
+void
+TensorProductPolynomialsConst<dim>::
+compute (const Point<dim> &p,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &grad_grads) const
+{
+ Assert (values.size()==n_tensor_pols || values.size()==0,
+ ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+ Assert (grads.size()==n_tensor_pols || grads.size()==0,
+ ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+ Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
+ ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
+
+ const bool update_values = (values.size() == n_tensor_pols),
+ update_grads = (grads.size()==n_tensor_pols),
+ update_grad_grads = (grad_grads.size()==n_tensor_pols);
+
+ // check how many
+ // values/derivatives we have to
+ // compute
+ unsigned int n_values_and_derivatives = 0;
+ if (update_values)
+ n_values_and_derivatives = 1;
+ if (update_grads)
+ n_values_and_derivatives = 2;
+ if (update_grad_grads)
+ n_values_and_derivatives = 3;
+
+
+ // compute the values (and derivatives, if
+ // necessary) of all polynomials at this
+ // evaluation point. to avoid many
+ // reallocation, use one std::vector for
+ // polynomial evaluation and store the
+ // result as Tensor<1,3> (that has enough
+ // fields for any evaluation of values and
+ // derivatives)
+ Table<2,Tensor<1,3> > v(dim, polynomials.size()+1);
+ {
+ std::vector<double> tmp (n_values_and_derivatives);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<polynomials.size(); ++i)
+ {
+ polynomials[i].value(p(d), tmp);
+ for (unsigned int e=0; e<n_values_and_derivatives; ++e)
+ v(d,i)[e] = tmp[e];
+ };
+ }
+
+ for (unsigned int i=0; i<n_tensor_pols-1; ++i)
+ {
+ // first get the
+ // one-dimensional indices of
+ // this particular tensor
+ // product polynomial
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ if (update_values)
+ {
+ values[i] = 1;
+ for (unsigned int x=0; x<dim; ++x)
+ values[i] *= v(x,indices[x])[0];
+ }
+
+ if (update_grads)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ grads[i][d] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grads[i][d] *= v(x,indices[x])[d==x];
+ }
+
+ if (update_grad_grads)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ grad_grads[i][d1][d2] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x || d2==x)
+ {
+ if (d1==d2)
+ derivative=2;
+ else
+ derivative=1;
+ }
+ grad_grads[i][d1][d2]
+ *= v(x,indices[x])[derivative];
+ }
+ }
+ }
+ //for dgq node: values =1, grads=0, grads_grads=0
+ if (update_values)
+ values[n_tensor_pols-1]=1;
+}
+
+
+
+template <int dim>
+unsigned int
+TensorProductPolynomialsConst<dim>::n() const
+{
+ return n_tensor_pols;
+}
+
+
+
+template <>
+unsigned int
+TensorProductPolynomialsConst<0>::n() const
+{
+ return numbers::invalid_unsigned_int;
+}
+
+/* ------------------- explicit instantiations -------------- */
+template class TensorProductPolynomialsConst<1>;
+template class TensorProductPolynomialsConst<2>;
+template class TensorProductPolynomialsConst<3>;
+
+DEAL_II_NAMESPACE_CLOSE
#include <deal.II/base/qprojector.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/tensor_product_polynomials_const.h>
#include <deal.II/base/polynomials_p.h>
#include <deal.II/fe/fe_poly.h>
#include <deal.II/fe/fe_values.h>
{
#if deal_II_dimension <= deal_II_space_dimension
template class FE_Poly<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
+ template class FE_Poly<TensorProductPolynomialsConst<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
#endif
--- /dev/null
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+
+#include <vector>
+#include <sstream>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+// namespace for some functions that are used in this file. they are
+// specific to numbering conventions used for the FE_Q element, and
+// are thus not very interesting to the outside world. we'd like to
+// simply put them into an anonymous namespace, but that triggers an
+// odd error with icc which can't compile this small snippet if the
+// function is static:
+// --------------------
+// template <int> struct int2type {};
+//
+// namespace {
+// static void SYMBOL (const int2type<1> & ) {}
+// }
+//
+// template <int dim, int spacedim> void g() {
+// SYMBOL(int2type<dim>());
+// }
+//
+// template void g<1>();
+// --------------------
+// the function needs to be static because of another icc bug, though.
+// work around this by packing everything into a namespace of its own
+// and have the anonymous namespace inside
+//
+// this is now intel icc issue 216082
+namespace FE_Q_DG0_Helper
+{
+ namespace
+ {
+ // given a permutation array,
+ // compute and return the inverse
+ // permutation
+#ifdef DEAL_II_ANON_NAMESPACE_BUG
+ static
+#endif
+ inline
+ std::vector<unsigned int>
+ invert_numbering (const std::vector<unsigned int> &in)
+ {
+ std::vector<unsigned int> out (in.size());
+ for (unsigned int i=0; i<in.size(); ++i)
+ out[in[i]]=i;
+ return out;
+ }
+ }
+}
+
+
+
+/**
+ * A class with the same purpose as the similarly named class of the
+ * Triangulation class. See there for more information.
+ */
+template <int xdim, int xspacedim>
+struct FE_Q_DG0<xdim,xspacedim>::Implementation
+{
+ /**
+ * Initialize the hanging node
+ * constraints matrices. Called from the
+ * constructor in case the finite element
+ * is based on quadrature points.
+ */
+ template <int spacedim>
+ static
+ void initialize_constraints (const Quadrature<1> &,
+ FE_Q_DG0<1,spacedim> &)
+ {
+ // no constraints in 1d
+ }
+
+
+ template <int spacedim>
+ static
+ void initialize_constraints (const Quadrature<1> &points,
+ FE_Q_DG0<2,spacedim> &fe)
+ {
+ const unsigned int dim = 2;
+
+ // restricted to each face, the
+ // traces of the shape functions is
+ // an element of P_{k} (in 2d), or
+ // Q_{k} (in 3d), where k is the
+ // degree of the element. from
+ // this, we interpolate between
+ // mother and cell face.
+
+ // the interpolation process works
+ // as follows: on each subface,
+ // we want that finite element
+ // solutions from both sides
+ // coincide. i.e. if a and b are
+ // expansion coefficients for the
+ // shape functions from both sides,
+ // we seek a relation between a and
+ // b such that
+ // sum_j a_j phi^c_j(x)
+ // == sum_j b_j phi_j(x)
+ // for all points x on the
+ // interface. here, phi^c_j are the
+ // shape functions on the small
+ // cell on one side of the face,
+ // and phi_j those on the big cell
+ // on the other side. To get this
+ // relation, it suffices to look at
+ // a sufficient number of points
+ // for which this has to hold. if
+ // there are n functions, then we
+ // need n evaluation points, and we
+ // choose them equidistantly.
+ //
+ // we obtain the matrix system
+ // A a == B b
+ // where
+ // A_ij = phi^c_j(x_i)
+ // B_ij = phi_j(x_i)
+ // and the relation we are looking
+ // for is
+ // a = A^-1 B b
+ //
+ // for the special case of Lagrange
+ // interpolation polynomials, A_ij
+ // reduces to delta_ij, and
+ // a_i = B_ij b_j
+ // Hence,
+ // interface_constraints(i,j)=B_ij.
+ //
+ // for the general case, where we
+ // don't have Lagrange
+ // interpolation polynomials, this
+ // is a little more
+ // complicated. Then we would
+ // evaluate at a number of points
+ // and invert the interpolation
+ // matrix A.
+ //
+ // Note, that we build up these
+ // matrices for all subfaces at
+ // once, rather than considering
+ // them separately. the reason is
+ // that we finally will want to
+ // have them in this order anyway,
+ // as this is the format we need
+ // inside deal.II
+
+ // In the following the points x_i
+ // are constructed in following
+ // order (n=degree-1)
+ // *----------*---------*
+ // 1..n 0 n+1..2n
+ // i.e. first the midpoint of the
+ // line, then the support points on
+ // subface 0 and on subface 1
+ std::vector<Point<dim-1> > constraint_points;
+ // Add midpoint
+ constraint_points.push_back (Point<dim-1> (0.5));
+
+ if (fe.degree>1)
+ {
+ const unsigned int n=fe.degree-1;
+ const double step=1./fe.degree;
+ // subface 0
+ for (unsigned int i=1; i<=n; ++i)
+ constraint_points.push_back (
+ GeometryInfo<dim-1>::child_to_cell_coordinates(Point<dim-1>(i*step),0));
+ // subface 1
+ for (unsigned int i=1; i<=n; ++i)
+ constraint_points.push_back (
+ GeometryInfo<dim-1>::child_to_cell_coordinates(Point<dim-1>(i*step),1));
+ }
+
+ // Now construct relation between
+ // destination (child) and source (mother)
+ // dofs.
+ const std::vector<Polynomials::Polynomial<double> > polynomials=
+ Polynomials::generate_complete_Lagrange_basis(points.get_points());
+
+ fe.interface_constraints
+ .TableBase<2,double>::reinit (fe.interface_constraints_size());
+
+ for (unsigned int i=0; i<constraint_points.size(); ++i)
+ for (unsigned int j=0; j<fe.degree+1; ++j)
+ {
+ fe.interface_constraints(i,j) =
+ polynomials[fe.face_index_map[j]].value (constraint_points[i](0));
+
+ // if the value is small up
+ // to round-off, then
+ // simply set it to zero to
+ // avoid unwanted fill-in
+ // of the constraint
+ // matrices (which would
+ // then increase the number
+ // of other DoFs a
+ // constrained DoF would
+ // couple to)
+ if (std::fabs(fe.interface_constraints(i,j)) < 1e-14)
+ fe.interface_constraints(i,j) = 0;
+ }
+ }
+
+
+ template <int spacedim>
+ static
+ void initialize_constraints (const Quadrature<1> &points,
+ FE_Q_DG0<3,spacedim> &fe)
+ {
+ const unsigned int dim = 3;
+
+ // For a detailed documentation of
+ // the interpolation see the
+ // FE_Q<2>::initialize_constraints
+ // function.
+
+ // In the following the points x_i
+ // are constructed in the order as
+ // described in the documentation
+ // of the FiniteElement class
+ // (fe_base.h), i.e.
+ // *--15--4--16--*
+ // | | |
+ // 10 19 6 20 12
+ // | | |
+ // 1--7---0--8---2
+ // | | |
+ // 9 17 5 18 11
+ // | | |
+ // *--13--3--14--*
+ std::vector<Point<dim-1> > constraint_points;
+
+ // Add midpoint
+ constraint_points.push_back (Point<dim-1> (0.5, 0.5));
+
+ // Add midpoints of lines of
+ // "mother-face"
+ constraint_points.push_back (Point<dim-1> (0, 0.5));
+ constraint_points.push_back (Point<dim-1> (1, 0.5));
+ constraint_points.push_back (Point<dim-1> (0.5, 0));
+ constraint_points.push_back (Point<dim-1> (0.5, 1));
+
+ if (fe.degree>1)
+ {
+ const unsigned int n=fe.degree-1;
+ const double step=1./fe.degree;
+ std::vector<Point<dim-2> > line_support_points(n);
+ for (unsigned int i=0; i<n; ++i)
+ line_support_points[i](0)=(i+1)*step;
+ Quadrature<dim-2> qline(line_support_points);
+
+ // auxiliary points in 2d
+ std::vector<Point<dim-1> > p_line(n);
+
+ // Add nodes of lines interior
+ // in the "mother-face"
+
+ // line 5: use line 9
+ QProjector<dim-1>::project_to_subface(qline, 0, 0, p_line);
+ for (unsigned int i=0; i<n; ++i)
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0.5, 0));
+ // line 6: use line 10
+ QProjector<dim-1>::project_to_subface(qline, 0, 1, p_line);
+ for (unsigned int i=0; i<n; ++i)
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0.5, 0));
+ // line 7: use line 13
+ QProjector<dim-1>::project_to_subface(qline, 2, 0, p_line);
+ for (unsigned int i=0; i<n; ++i)
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
+ // line 8: use line 14
+ QProjector<dim-1>::project_to_subface(qline, 2, 1, p_line);
+ for (unsigned int i=0; i<n; ++i)
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
+
+ // DoFs on bordering lines
+ // lines 9-16
+ for (unsigned int face=0; face<GeometryInfo<dim-1>::faces_per_cell; ++face)
+ for (unsigned int subface=0;
+ subface<GeometryInfo<dim-1>::max_children_per_face; ++subface)
+ {
+ QProjector<dim-1>::project_to_subface(qline, face, subface, p_line);
+ constraint_points.insert(constraint_points.end(),
+ p_line.begin(), p_line.end());
+ }
+
+ // Create constraints for
+ // interior nodes
+ std::vector<Point<dim-1> > inner_points(n*n);
+ for (unsigned int i=0, iy=1; iy<=n; ++iy)
+ for (unsigned int ix=1; ix<=n; ++ix)
+ inner_points[i++] = Point<dim-1> (ix*step, iy*step);
+
+ // at the moment do this for
+ // isotropic face refinement only
+ for (unsigned int child=0;
+ child<GeometryInfo<dim-1>::max_children_per_cell; ++child)
+ for (unsigned int i=0; i<inner_points.size(); ++i)
+ constraint_points.push_back (
+ GeometryInfo<dim-1>::child_to_cell_coordinates(inner_points[i], child));
+ }
+
+ // Now construct relation between
+ // destination (child) and source (mother)
+ // dofs.
+ const unsigned int pnts=(fe.degree+1)*(fe.degree+1);
+ const std::vector<Polynomials::Polynomial<double> > polynomial_basis=
+ Polynomials::generate_complete_Lagrange_basis(points.get_points());
+
+ const TensorProductPolynomialsConst<dim-1> face_polynomials(polynomial_basis);
+
+ fe.interface_constraints
+ .TableBase<2,double>::reinit (fe.interface_constraints_size());
+
+ for (unsigned int i=0; i<constraint_points.size(); ++i)
+ {
+ const double interval = (double) (fe.degree * 2);
+ bool mirror[dim - 1];
+ Point<dim-1> constraint_point;
+
+ // Eliminate FP errors in constraint
+ // points. Due to their origin, they
+ // must all be fractions of the unit
+ // interval. If we have polynomial
+ // degree 4, the refined element has 8
+ // intervals. Hence the coordinates
+ // must be 0, 0.125, 0.25, 0.375 etc.
+ // Now the coordinates of the
+ // constraint points will be multiplied
+ // by the inverse of the interval size
+ // (in the example by 8). After that
+ // the coordinates must be integral
+ // numbers. Hence a normal truncation
+ // is performed and the coordinates
+ // will be scaled back. The equal
+ // treatment of all coordinates should
+ // eliminate any FP errors.
+ for (unsigned int k=0; k<dim-1; ++k)
+ {
+ const int coord_int =
+ static_cast<int> (constraint_points[i](k) * interval + 0.25);
+ constraint_point(k) = 1.*coord_int / interval;
+
+ // The following lines of code
+ // should eliminate the problems
+ // with the Constraint-Matrix,
+ // which appeared for P>=4. The
+ // Constraint-Matrix class
+ // complained about different
+ // constraints for the same entry
+ // of the Constraint-Matrix.
+ // Actually this difference could
+ // be attributed to FP errors, as
+ // it was in the range of
+ // 1.0e-16. These errors originate
+ // in the loss of symmetry in the
+ // FP approximation of the
+ // shape-functions. Considering a
+ // 3rd order shape function in 1D,
+ // we have N0(x)=N3(1-x) and
+ // N1(x)=N2(1-x). For higher order
+ // polynomials the FP
+ // approximations of the shape
+ // functions do not satisfy these
+ // equations any more! Thus in the
+ // following code everything is
+ // computed in the interval x \in
+ // [0..0.5], which is sufficient to
+ // express all values that could
+ // come out from a computation of
+ // any shape function in the full
+ // interval [0..1]. If x > 0.5 the
+ // computation is done for 1-x with
+ // the shape function N_{p-n}
+ // instead of N_n. Hence symmetry
+ // is preserved and everything
+ // works fine...
+ //
+ // For a different explanation of
+ // the problem, see the discussion
+ // in the FiniteElement class
+ // for constraint matrices in 3d.
+ mirror[k] = (constraint_point(k) > 0.5);
+ if (mirror[k])
+ constraint_point(k) = 1.0 - constraint_point(k);
+ }
+
+ for (unsigned int j=0; j<pnts; ++j)
+ {
+ unsigned int indices[2]
+ = { fe.face_index_map[j] % (fe.degree + 1),
+ fe.face_index_map[j] / (fe.degree + 1)
+ };
+
+ for (unsigned int k = 0; k<2; ++k)
+ if (mirror[k])
+ indices[k] = fe.degree - indices[k];
+
+ const unsigned int
+ new_index = indices[1] * (fe.degree + 1) + indices[0];
+
+ fe.interface_constraints(i,j) =
+ face_polynomials.compute_value (new_index, constraint_point);
+
+ // if the value is small up
+ // to round-off, then
+ // simply set it to zero to
+ // avoid unwanted fill-in
+ // of the constraint
+ // matrices (which would
+ // then increase the number
+ // of other DoFs a
+ // constrained DoF would
+ // couple to)
+ if (std::fabs(fe.interface_constraints(i,j)) < 1e-14)
+ fe.interface_constraints(i,j) = 0;
+ }
+ }
+ }
+};
+
+
+template <int dim, int spacedim>
+FE_Q_DG0<dim,spacedim>::FE_Q_DG0 (const unsigned int degree)
+ :
+ FE_Poly<TensorProductPolynomialsConst<dim>, dim, spacedim> (
+ TensorProductPolynomialsConst<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
+ FiniteElementData<dim>(get_dpo_vector(degree),
+ 1, degree,
+ FiniteElementData<dim>::L2),
+ get_riaf_vector(degree),
+ std::vector<ComponentMask>(1ll , std::vector<bool>(1,true))),
+ face_index_map(FE_Q_DG0_Helper::invert_numbering(
+ face_lexicographic_to_hierarchic_numbering (degree)))
+{
+ Assert (degree > 0,
+ ExcMessage ("This element can only be used for polynomial degrees "
+ "greater than zero"));
+
+ std::vector<unsigned int> renumber (this->dofs_per_cell-1);
+ const FiniteElementData<dim> fe(get_dpo_vector_q(degree),1,degree);
+ FETools::hierarchic_to_lexicographic_numbering (fe, renumber);
+ renumber.push_back(this->dofs_per_cell-1);
+ this->poly_space.set_numbering(renumber);
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points ();
+ initialize_unit_face_support_points ();
+
+ // reinit constraints
+ initialize_constraints ();
+
+ //TODO not working
+ // Reinit the vectors of restriction and
+ // prolongation matrices to the right sizes
+ // and compute the matrices
+ /*this->reinit_restriction_and_prolongation_matrices();
+ if (dim == spacedim)
+ {
+ FETools::compute_embedding_matrices (*this, this->prolongation);
+ FETools::compute_projection_matrices (*this, this->restriction);
+ }
+ else
+ {
+ FE_Q_DG0<dim> tmp (degree);
+ this->prolongation = tmp.prolongation;
+ this->restriction = tmp.restriction;
+ }*/
+
+ initialize_quad_dof_index_permutation();
+}
+
+
+
+template <int dim, int spacedim>
+FE_Q_DG0<dim,spacedim>::FE_Q_DG0 (const Quadrature<1> &points)
+ :
+ FE_Poly<TensorProductPolynomialsConst<dim>, dim, spacedim> (
+ TensorProductPolynomialsConst<dim>(Polynomials::generate_complete_Lagrange_basis(points.get_points())),
+ FiniteElementData<dim>(get_dpo_vector(points.size()-1),
+ 1, points.size()-1,
+ FiniteElementData<dim>::L2),
+ get_riaf_vector(points.size()-1),
+ std::vector<ComponentMask>(1, std::vector<bool>(1,true))),
+
+ face_index_map(FE_Q_DG0_Helper::invert_numbering(
+ face_lexicographic_to_hierarchic_numbering (points.size()-1)))
+{
+ const int degree = points.size()-1;
+
+ Assert (degree > 0,
+ ExcMessage ("This element can only be used for polynomial degrees "
+ "at least zero"));
+ Assert (points.point(0)(0) == 0,
+ ExcMessage ("The first support point has to be zero."));
+ Assert (points.point(degree)(0) == 1,
+ ExcMessage ("The last support point has to be one."));
+
+ std::vector<unsigned int> renumber (this->dofs_per_cell-1);
+ const FiniteElementData<dim> fe(get_dpo_vector_q(degree),1,degree);
+ FETools::hierarchic_to_lexicographic_numbering (fe, renumber);
+ renumber.push_back(this->dofs_per_cell-1);
+ this->poly_space.set_numbering(renumber);
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points (points);
+ initialize_unit_face_support_points (points);
+
+ // reinit constraints
+ Implementation::initialize_constraints (points, *this);
+
+ // Reinit the vectors of restriction and
+ // prolongation matrices to the right sizes
+ // and compute the matrices
+ this->reinit_restriction_and_prolongation_matrices();
+
+ initialize_quad_dof_index_permutation();
+}
+
+
+
+template <int dim, int spacedim>
+std::string
+FE_Q_DG0<dim,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;
+ bool type = true;
+ const unsigned int n_points = this->degree +1;
+ std::vector<double> points(n_points);
+ const unsigned int dofs_per_cell = this->dofs_per_cell;
+ const std::vector<Point<dim> > &unit_support_points = this->unit_support_points;
+ unsigned int index = 0;
+
+ // Decode the support points
+ // in one coordinate direction.
+ for (unsigned int j=0; j<dofs_per_cell; j++)
+ {
+ if ((dim>1) ? (unit_support_points[j](1)==0 &&
+ ((dim>2) ? unit_support_points[j](2)==0: true)) : true)
+ {
+ if (index == 0)
+ points[index] = unit_support_points[j](0);
+ else if (index == 1)
+ points[n_points-1] = unit_support_points[j](0);
+ else
+ points[index-1] = unit_support_points[j](0);
+
+ index++;
+ }
+ }
+ //Do not consider the discontinuous node for dimension 1
+ Assert (index == n_points || (dim==1 && index == n_points+1),
+ ExcMessage ("Could not decode support points in one coordinate direction."));
+
+ // Check whether the support
+ // points are equidistant.
+ for (unsigned int j=0; j<n_points; j++)
+ if (std::fabs(points[j] - (double)j/this->degree) > 1e-15)
+ {
+ type = false;
+ break;
+ }
+
+ if (type == true)
+ namebuf << "FE_Q_DG0<" << dim << ">(" << this->degree << ")";
+ else
+ {
+
+ // Check whether the support
+ // points come from QGaussLobatto.
+ const QGaussLobatto<1> points_gl(n_points);
+ type = true;
+ for (unsigned int j=0; j<n_points; j++)
+ if (points[j] != points_gl.point(j)(0))
+ {
+ type = false;
+ break;
+ }
+ if (type == true)
+ namebuf << "FE_Q_DG0<" << dim << ">(QGaussLobatto(" << this->degree+1 << "))";
+ else
+ namebuf << "FE_Q_DG0<" << dim << ">(QUnknownNodes(" << this->degree << "))";
+ }
+ return namebuf.str();
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_Q_DG0<dim,spacedim>::clone() const
+{
+ return new FE_Q_DG0<dim,spacedim>(*this);
+}
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::interpolate(std::vector<double> &local_dofs,
+ const std::vector<double> &values) const
+{
+ Assert (values.size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (this->n_components() == 1,
+ ExcDimensionMismatch(this->n_components(), 1));
+
+ std::copy(values.begin(), values.end(), local_dofs.begin());
+ //We don't need the discontinuous function for local interpolation
+ local_dofs[local_dofs.size()-1]=0.;
+}
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::interpolate(std::vector<double> &local_dofs,
+ const std::vector<Vector<double> > &values,
+ unsigned int offset) const
+{
+ Assert (values.size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values[0].size() >= offset+this->n_components(),
+ ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
+
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[i](offset+index.first);
+ }
+ //We don't need the discontinuous function for local interpolation
+ local_dofs[local_dofs.size()-1]=0.;
+}
+
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::interpolate(
+ std::vector<double> &local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > > &values) const
+{
+ Assert (values[0].size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values.size() == this->n_components(),
+ ExcDimensionMismatch(values.size(), this->n_components()));
+
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[index.first][i];
+ }
+ //We don't need the discontinuous function for local interpolation
+ local_dofs[local_dofs.size()-1]=0.;
+}
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::
+get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this is only implemented, if the
+ // source FE is also a
+ // Q_DG0 element
+ typedef FE_Q_DG0<dim,spacedim> FEQDG0;
+ typedef FiniteElement<dim,spacedim> FEL;
+
+ AssertThrow ((x_source_fe.get_name().find ("FE_Q_DG0<") == 0)
+ ||
+ (dynamic_cast<const FEQDG0 *>(&x_source_fe) != 0),
+ typename FEL::
+ ExcInterpolationNotImplemented());
+
+ Assert (interpolation_matrix.m() == this->dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ this->dofs_per_cell));
+ Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_cell));
+
+ // ok, source is a Q_DG0 element, so
+ // we will be able to do the work
+ const FE_Q_DG0<dim,spacedim> &source_fe
+ = dynamic_cast<const FE_Q_DG0<dim,spacedim>&>(x_source_fe);
+
+ // compute the interpolation
+ // matrices in much the same way as
+ // we do for the embedding matrices
+ // from mother to child.
+ FullMatrix<double> cell_interpolation (this->dofs_per_cell,
+ this->dofs_per_cell);
+ FullMatrix<double> source_interpolation (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ FullMatrix<double> tmp (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ for (unsigned int j=0; j<this->dofs_per_cell-1; ++j)
+ {
+ // read in a point on this
+ // cell and evaluate the
+ // shape functions there
+ const Point<dim> p = this->unit_support_points[j];
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ cell_interpolation(j,i) = this->poly_space.compute_value (i, p);
+
+ for (unsigned int i=0; i<source_fe.dofs_per_cell-1; ++i)
+ source_interpolation(j,i) = source_fe.poly_space.compute_value (i,p);
+
+ }
+ //the discontinuous node is not transformed
+ cell_interpolation(this->dofs_per_cell-1,
+ this->dofs_per_cell-1)=1.;
+ source_interpolation(this->dofs_per_cell-1,
+ source_fe.dofs_per_cell-1)=1.;
+
+ // then compute the
+ // interpolation matrix
+ // for this coordinate
+ // direction
+ cell_interpolation.gauss_jordan ();
+ cell_interpolation.mmult (interpolation_matrix,
+ source_interpolation);
+
+ const double eps = 2e-13*this->degree*dim;
+
+ // cut off very small values
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ if (std::fabs(interpolation_matrix(i,j)) < eps)
+ interpolation_matrix(i,j) = 0.;
+
+ // make sure that the row sum of
+ // each of the matrices is 1 at
+ // this point. this must be so
+ // since the shape functions sum up
+ // to 1
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ {
+ double sum = 0.;
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ sum += interpolation_matrix(i,j);
+
+ Assert (std::fabs(sum-1) < eps, ExcInternalError());
+ }
+}
+
+
+
+template <>
+void
+FE_Q_DG0<1>::
+get_face_interpolation_matrix (const FiniteElement<1,1> &/*x_source_fe*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+}
+
+
+
+template <>
+void
+FE_Q_DG0<1>::
+get_subface_interpolation_matrix (const FiniteElement<1,1> &/*x_source_fe*/,
+ const unsigned int /*subface*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+}
+
+
+template <>
+void
+FE_Q_DG0<1,2>::
+get_face_interpolation_matrix (const FiniteElement<1,2> &/*x_source_fe*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ typedef FiniteElement<1,2> FEL;
+ Assert (false,
+ FEL::
+ ExcInterpolationNotImplemented ());
+}
+
+
+template <>
+void
+FE_Q_DG0<1,2>::
+get_subface_interpolation_matrix (const FiniteElement<1,2> &/*x_source_fe*/,
+ const unsigned int /*subface*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ typedef FiniteElement<1,2> FEL;
+ Assert (false,
+ FEL::
+ ExcInterpolationNotImplemented ());
+}
+
+
+template <>
+void
+FE_Q_DG0<1,3>::
+get_face_interpolation_matrix (const FiniteElement<1,3> &/*x_source_fe*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ typedef FiniteElement<1,3> FEL;
+ Assert (false,
+ FEL::
+ ExcInterpolationNotImplemented ());
+}
+
+
+template <>
+void
+FE_Q_DG0<1,3>::
+get_subface_interpolation_matrix (const FiniteElement<1,3> &/*x_source_fe*/,
+ const unsigned int /*subface*/,
+ FullMatrix<double> &/*interpolation_matrix*/) const
+{
+ typedef FiniteElement<1,3> FEL;
+ Assert (false,
+ FEL::
+ ExcInterpolationNotImplemented ());
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this is only implemented, if the
+ // source FE is also a
+ // Q_DG0 element
+ typedef FE_Q_DG0<dim,spacedim> FEQDG0;
+ typedef FiniteElement<dim,spacedim> FEL;
+ AssertThrow ((x_source_fe.get_name().find ("FE_Q_DG0<") == 0)
+ ||
+ (dynamic_cast<const FEQDG0 *>(&x_source_fe) != 0),
+ typename FEL::
+ ExcInterpolationNotImplemented());
+
+ 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));
+
+ // ok, source is a Q_DG0 element, so
+ // we will be able to do the work
+ const FE_Q_DG0<dim,spacedim> &source_fe
+ = dynamic_cast<const FE_Q_DG0<dim,spacedim>&>(x_source_fe);
+
+ // Make sure, that the element,
+ // for which the DoFs should be
+ // constrained is the one with
+ // the higher polynomial degree.
+ // Actually the procedure will work
+ // also if this assertion is not
+ // satisfied. But the matrices
+ // produced in that case might
+ // lead to problems in the
+ // hp procedures, which use this
+ // method.
+ Assert (this->dofs_per_face <= source_fe.dofs_per_face,
+ typename FEL::
+ ExcInterpolationNotImplemented ());
+
+ // generate a quadrature
+ // with the unit support points.
+ // This is later based as a
+ // basis for the QProjector,
+ // which returns the support
+ // points on the face.
+ Quadrature<dim-1> quad_face_support (source_fe.get_unit_face_support_points ());
+
+ // Rule of thumb for FP accuracy,
+ // that can be expected for a
+ // given polynomial degree.
+ // This value is used to cut
+ // off values close to zero.
+ const double eps = 2e-13*this->degree*(dim-1);
+
+ // compute the interpolation
+ // matrix by simply taking the
+ // value at the support points.
+//TODO: Verify that all faces are the same with respect to
+// these support points. Furthermore, check if something has to
+// be done for the face orientation flag in 3D.
+ const Quadrature<dim> face_quadrature
+ = QProjector<dim>::project_to_face (quad_face_support, 0);
+ for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
+ {
+ const Point<dim> &p = face_quadrature.point (i);
+
+ for (unsigned int j=0; j<this->dofs_per_face; ++j)
+ {
+ double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
+
+ // Correct the interpolated
+ // value. I.e. if it is close
+ // to 1 or 0, make it exactly
+ // 1 or 0. Unfortunately, this
+ // is required to avoid problems
+ // with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i,j) = matrix_entry;
+ }
+ }
+
+ // make sure that the row sum of
+ // each of the matrices is 1 at
+ // this point. this must be so
+ // since the shape functions sum up
+ // to 1
+ for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ sum += interpolation_matrix(j,i);
+
+ Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1),
+ ExcInternalError());
+ }
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_face));
+
+ // see if source is a Q_DG0 element
+ if (const FE_Q_DG0<dim,spacedim> *source_fe
+ = dynamic_cast<const FE_Q_DG0<dim,spacedim> *>(&x_source_fe))
+ {
+ // have this test in here since
+ // a table of size 2x0 reports
+ // its size as 0x0
+ Assert (interpolation_matrix.n() == this->dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.n(),
+ this->dofs_per_face));
+
+ // Make sure, that the element,
+ // for which the DoFs should be
+ // constrained is the one with
+ // the higher polynomial degree.
+ // Actually the procedure will work
+ // also if this assertion is not
+ // satisfied. But the matrices
+ // produced in that case might
+ // lead to problems in the
+ // hp procedures, which use this
+ // method.
+ Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented ()));
+
+ // generate a point on this
+ // cell and evaluate the
+ // shape functions there
+ const Quadrature<dim-1>
+ quad_face_support (source_fe->get_unit_face_support_points ());
+
+ // Rule of thumb for FP accuracy,
+ // that can be expected for a
+ // given polynomial degree.
+ // This value is used to cut
+ // off values close to zero.
+ double eps = 2e-13*this->degree*(dim-1);
+
+ // compute the interpolation
+ // matrix by simply taking the
+ // value at the support points.
+//TODO: Verify that all faces are the same with respect to
+// these support points. Furthermore, check if something has to
+// be done for the face orientation flag in 3D.
+ const Quadrature<dim> subface_quadrature
+ = QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
+ for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
+ {
+ const Point<dim> &p = subface_quadrature.point (i);
+
+ for (unsigned int j=0; j<this->dofs_per_face; ++j)
+ {
+ double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
+
+ // Correct the interpolated
+ // value. I.e. if it is close
+ // to 1 or 0, make it exactly
+ // 1 or 0. Unfortunately, this
+ // is required to avoid problems
+ // with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i,j) = matrix_entry;
+ }
+ }
+
+ // make sure that the row sum of
+ // each of the matrices is 1 at
+ // this point. this must be so
+ // since the shape functions sum up
+ // to 1
+ for (unsigned int j=0; j<source_fe->dofs_per_face; ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ sum += interpolation_matrix(j,i);
+
+ Assert (std::fabs(sum-1) < 2e-13*this->degree*dim,
+ ExcInternalError());
+ }
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
+ {
+ // nothing to do here, the
+ // FE_Nothing has no degrees of
+ // freedom anyway
+ }
+ else
+ AssertThrow (false,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented()));
+}
+
+
+
+template <int dim, int spacedim>
+bool
+FE_Q_DG0<dim,spacedim>::hp_constraints_are_implemented () const
+{
+ return true;
+}
+
+
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Q_DG0<dim,spacedim>::
+hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ // we can presently only compute
+ // these identities if both FEs are
+ // FE_Q_DG0s or if the other one is an
+ // FE_Nothing. in the first case,
+ // there should be exactly one
+ // single DoF of each FE at a
+ // vertex, and they should have
+ // identical value
+ if (dynamic_cast<const FE_Q_DG0<dim,spacedim>*>(&fe_other) != 0)
+ {
+ return
+ std::vector<std::pair<unsigned int, unsigned int> >
+ (1, std::make_pair (0U, 0U));
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no
+ // degrees of freedom, so there
+ // are no equivalencies to be
+ // recorded
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+ else
+ {
+ Assert (false, ExcNotImplemented());
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Q_DG0<dim,spacedim>::
+hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ // we can presently only compute
+ // these identities if both FEs are
+ // FE_Q_DG0s or if the other one is an
+ // FE_Nothing
+ if (const FE_Q_DG0<dim,spacedim> *fe_q_dg0_other = dynamic_cast<const
+FE_Q_DG0<dim,spacedim>*>(&fe_other))
+ {
+ // dofs are located along lines, so two
+ // dofs are identical if they are
+ // located at identical positions. if
+ // we had only equidistant points, we
+ // could simple check for similarity
+ // like (i+1)*q == (j+1)*p, but we
+ // might have other support points
+ // (e.g. Gauss-Lobatto
+ // points). Therefore, read the points
+ // in unit_support_points for the first
+ // coordinate direction. We take the
+ // lexicographic ordering of the points
+ // in the first direction (i.e.,
+ // x-direction), which we access
+ // between index 1 and p-1 (index 0 and
+ // p are vertex dofs).
+ const unsigned int p = this->degree;
+ const unsigned int q = fe_q_dg0_other->degree;
+
+ std::vector<std::pair<unsigned int, unsigned int> > identities;
+
+ const std::vector<unsigned int> &index_map_inverse=
+ this->poly_space.get_numbering_inverse();
+ const std::vector<unsigned int> &index_map_inverse_other=
+ fe_q_dg0_other->poly_space.get_numbering_inverse();
+
+ for (unsigned int i=0; i<p-1; ++i)
+ for (unsigned int j=0; j<q-1; ++j)
+ if (std::fabs(this->unit_support_points[index_map_inverse[i+1]][0]-
+
+fe_q_dg0_other->unit_support_points[index_map_inverse_other[j+1]][0])
+ < 1e-14)
+ identities.push_back (std::make_pair(i,j));
+
+ return identities;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no
+ // degrees of freedom, so there
+ // are no equivalencies to be
+ // recorded
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+ else
+ {
+ Assert (false, ExcNotImplemented());
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_Q_DG0<dim,spacedim>::
+hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ // we can presently only compute
+ // these identities if both FEs are
+ // FE_Q_DG0s or if the other one is an
+ // FE_Nothing
+ if (const FE_Q_DG0<dim,spacedim> *fe_q_dg0_other = dynamic_cast<const
+FE_Q_DG0<dim,spacedim>*>(&fe_other))
+ {
+ // this works exactly like the line
+ // case above, except that now we have
+ // to have two indices i1, i2 and j1,
+ // j2 to characterize the dofs on the
+ // face of each of the finite
+ // elements. since they are ordered
+ // lexicographically along the first
+ // line and we have a tensor product,
+ // the rest is rather straightforward
+ const unsigned int p = this->degree;
+ const unsigned int q = fe_q_dg0_other->degree;
+
+ std::vector<std::pair<unsigned int, unsigned int> > identities;
+
+ const std::vector<unsigned int> &index_map_inverse=
+ this->poly_space.get_numbering_inverse();
+ const std::vector<unsigned int> &index_map_inverse_other=
+ fe_q_dg0_other->poly_space.get_numbering_inverse();
+
+ for (unsigned int i1=0; i1<p-1; ++i1)
+ for (unsigned int i2=0; i2<p-1; ++i2)
+ for (unsigned int j1=0; j1<q-1; ++j1)
+ for (unsigned int j2=0; j2<q-1; ++j2)
+ if ((std::fabs(this->unit_support_points[index_map_inverse[i1+1]][0]-
+
+fe_q_dg0_other->unit_support_points[index_map_inverse_other[j1+1]][0])
+ < 1e-14)
+ &&
+ (std::fabs(this->unit_support_points[index_map_inverse[i2+1]][0]-
+
+fe_q_dg0_other->unit_support_points[index_map_inverse_other[j2+1]][0])
+ < 1e-14))
+ identities.push_back (std::make_pair(i1*(p-1)+i2,
+ j1*(q-1)+j2));
+
+ return identities;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no
+ // degrees of freedom, so there
+ // are no equivalencies to be
+ // recorded
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+ else
+ {
+ Assert (false, ExcNotImplemented());
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_Q_DG0<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ if (const FE_Q_DG0<dim,spacedim> *fe_q_dg0_other
+ = dynamic_cast<const FE_Q_DG0<dim,spacedim>*>(&fe_other))
+ {
+ if (this->degree < fe_q_dg0_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_dg0_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no
+ // degrees of freedom and it is
+ // typically used in a context
+ // where we don't require any
+ // continuity along the
+ // interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert (false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+}
+
+
+//---------------------------------------------------------------------------
+// Auxiliary functions
+//---------------------------------------------------------------------------
+
+
+
+template <int dim, int spacedim>
+void FE_Q_DG0<dim,spacedim>::initialize_unit_support_points ()
+{
+ // number of points: (degree+1)^dim+1
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= this->degree+1;
+ n++;
+
+ this->unit_support_points.resize(n);
+
+ const std::vector<unsigned int> &index_map_inverse=
+ this->poly_space.get_numbering_inverse();
+
+ const double step = 1./this->degree;
+ Point<dim> p;
+
+ unsigned int k=0;
+ for (unsigned int iz=0; iz <= ((dim>2) ? this->degree : 0) ; ++iz)
+ for (unsigned int iy=0; iy <= ((dim>1) ? this->degree : 0) ; ++iy)
+ for (unsigned int ix=0; ix<=this->degree; ++ix)
+ {
+ p(0) = ix * step;
+ if (dim>1)
+ p(1) = iy * step;
+ if (dim>2)
+ p(2) = iz * step;
+
+ this->unit_support_points[index_map_inverse[k++]] = p;
+ }
+ //dg0 support point in the center of the cell
+ p(0) = .5;
+ if (dim>1)
+ {
+ p(1) = .5;
+ if (dim>2)
+ p(2) = .5;
+ }
+
+ this->unit_support_points[index_map_inverse[k]] = p;
+}
+
+
+
+template <int dim, int spacedim>
+void FE_Q_DG0<dim,spacedim>::initialize_unit_support_points (const Quadrature<1> &points)
+{
+ // number of points: (degree+1)^dim+1
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= this->degree+1;
+ n++;
+
+ this->unit_support_points.resize(n);
+
+ const std::vector<unsigned int> &index_map_inverse=
+ this->poly_space.get_numbering_inverse();
+
+ Quadrature<dim> support_quadrature(points);
+
+ Point<dim> p;
+
+ for (unsigned int k=0; k<n ; k++)
+ {
+ this->unit_support_points[index_map_inverse[k]] = support_quadrature.point(k);
+ }
+}
+
+
+
+template <>
+void FE_Q_DG0<1>::initialize_unit_face_support_points ()
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <>
+void FE_Q_DG0<1>::initialize_unit_face_support_points (const Quadrature<1> &/*points*/)
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <>
+void FE_Q_DG0<1,2>::initialize_unit_face_support_points ()
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <>
+void FE_Q_DG0<1,2>::initialize_unit_face_support_points (const Quadrature<1> &/*points*/)
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <>
+void FE_Q_DG0<1,3>::initialize_unit_face_support_points ()
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <>
+void FE_Q_DG0<1,3>::initialize_unit_face_support_points (const Quadrature<1> &/*points*/)
+{
+ // no faces in 1d, so nothing to do
+}
+
+template <int dim, int spacedim>
+void FE_Q_DG0<dim,spacedim>::initialize_unit_face_support_points ()
+{
+ const unsigned int codim = dim-1;
+
+ // number of points: (degree+1)^codim
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<codim; ++i)
+ n *= this->degree+1;
+
+ this->unit_face_support_points.resize(n);
+
+ const std::vector<unsigned int> &face_index_map_inverse=
+ FE_Q_DG0_Helper::invert_numbering(face_index_map);
+
+ const double step = 1./this->degree;
+ Point<codim> p;
+
+ unsigned int k=0;
+ for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
+ for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
+ for (unsigned int ix=0; ix<=this->degree; ++ix)
+ {
+ p(0) = ix * step;
+ if (codim>1)
+ p(1) = iy * step;
+ if (codim>2)
+ p(2) = iz * step;
+
+ this->unit_face_support_points[face_index_map_inverse[k++]] = p;
+ }
+}
+
+
+
+template <int dim, int spacedim>
+void FE_Q_DG0<dim,spacedim>::initialize_unit_face_support_points (const Quadrature<1> &points)
+{
+ const unsigned int codim = dim-1;
+
+ // number of points: (degree+1)^codim
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<codim; ++i)
+ n *= this->degree+1;
+
+ this->unit_face_support_points.resize(n);
+
+ const std::vector< Point<1> > edge = points.get_points();
+
+ const std::vector<unsigned int> &face_index_map_inverse=
+ FE_Q_DG0_Helper::invert_numbering(face_index_map);
+
+ Point<codim> p;
+
+ unsigned int k=0;
+ for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
+ for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
+ for (unsigned int ix=0; ix<=this->degree; ++ix)
+ {
+ p(0) = edge[ix](0);
+ if (codim>1)
+ p(1) = edge[iy](0);
+ if (codim>2)
+ p(2) = edge[iz](0);
+
+ this->unit_face_support_points[face_index_map_inverse[k++]] = p;
+ }
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::initialize_quad_dof_index_permutation ()
+{
+ // general template for 1D and 2D, do nothing
+}
+
+
+
+template <>
+void
+FE_Q_DG0<3>::initialize_quad_dof_index_permutation ()
+{
+ Assert (adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad,
+ ExcInternalError());
+
+ const unsigned int n=this->degree-1;
+ Assert(n*n==this->dofs_per_quad, ExcInternalError());
+
+ // alias for the table to fill
+ Table<2,int> &data=this->adjust_quad_dof_index_for_face_orientation_table;
+
+ // the dofs on a face are connected to a n x
+ // n matrix. for example, for degree==4 we
+ // have the following dofs on a quad
+
+ // ___________
+ // | |
+ // | 6 7 8 |
+ // | |
+ // | 3 4 5 |
+ // | |
+ // | 0 1 2 |
+ // |___________|
+ //
+ // we have dof_no=i+n*j with index i in
+ // x-direction and index j in y-direction
+ // running from 0 to n-1. to extract i and j
+ // we can use i=dof_no%n and j=dof_no/n. i
+ // and j can then be used to construct the
+ // rotated and mirrored numbers.
+
+
+ for (unsigned int local=0; local<this->dofs_per_quad; ++local)
+ // face support points are in lexicographic
+ // ordering with x running fastest. invert
+ // that (y running fastest)
+ {
+ unsigned int i=local%n,
+ j=local/n;
+
+ // face_orientation=false, face_flip=false, face_rotation=false
+ data(local,0)=j + i *n - local;
+ // face_orientation=false, face_flip=false, face_rotation=true
+ data(local,1)=i + (n-1-j)*n - local;
+ // face_orientation=false, face_flip=true, face_rotation=false
+ data(local,2)=(n-1-j) + (n-1-i)*n - local;
+ // face_orientation=false, face_flip=true, face_rotation=true
+ data(local,3)=(n-1-i) + j *n - local;
+ // face_orientation=true, face_flip=false, face_rotation=false
+ data(local,4)=0;
+ // face_orientation=true, face_flip=false, face_rotation=true
+ data(local,5)=j + (n-1-i)*n - local;
+ // face_orientation=true, face_flip=true, face_rotation=false
+ data(local,6)=(n-1-i) + (n-1-j)*n - local;
+ // face_orientation=true, face_flip=true, face_rotation=true
+ data(local,7)=(n-1-j) + i *n - local;
+ }
+
+ // aditionally initialize reordering of line
+ // dofs
+ for (unsigned int i=0; i<this->dofs_per_line; ++i)
+ this->adjust_line_dof_index_for_line_orientation_table[i]=this->dofs_per_line-1-i - i;
+}
+
+template <int dim, int spacedim>
+std::vector<bool>
+FE_Q_DG0<dim,spacedim>::get_riaf_vector(const unsigned int deg)
+{
+ std::vector<bool> riaf
+ (FiniteElementData<dim> (get_dpo_vector(deg),1,deg).dofs_per_cell,
+ false);
+ riaf[riaf.size()-1]=true;
+ return riaf;
+}
+
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_Q_DG0<dim,spacedim>::get_dpo_vector(const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(deg-1);
+
+ dpo[dim]++;//we need an additional DG0-node for a dim-dimensional object
+ return dpo;
+}
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_Q_DG0<dim,spacedim>::get_dpo_vector_q(const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(deg-1);
+
+ return dpo;
+}
+
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_Q_DG0<dim,spacedim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+{
+ const FiniteElementData<dim-1> face_data(FE_Q_DG0<dim-1>::get_dpo_vector_q(degree),1,degree);
+ std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
+ FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
+ return face_renumber;
+}
+
+
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+ return std::vector<unsigned int>();
+}
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1,2>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+ return std::vector<unsigned int>();
+}
+
+
+template <>
+std::vector<unsigned int>
+FE_Q_DG0<1,3>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+ return std::vector<unsigned int>();
+}
+
+
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::initialize_constraints ()
+{
+ QTrapez<1> trapez;
+ QIterated<1> points (trapez,this->degree);
+ Implementation::initialize_constraints (points, *this);
+}
+
+//---------------------------------------------------------------------------
+// Data field initialization
+//---------------------------------------------------------------------------
+
+
+template <>
+bool
+FE_Q_DG0<1>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+ Assert (face_index < GeometryInfo<1>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell));
+
+
+ // in 1d, things are simple. since
+ // there is only one degree of
+ // freedom per vertex in this
+ // class, the first is on vertex 0
+ // (==face 0 in some sense), the
+ // second on face 1:
+ return (((shape_index == 0) && (face_index == 0)) ||
+ ((shape_index == 1) && (face_index == 1)));
+}
+
+
+template <>
+bool
+FE_Q_DG0<1,2>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+ Assert (face_index < GeometryInfo<1>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell));
+
+
+ // in 1d, things are simple. since
+ // there is only one degree of
+ // freedom per vertex in this
+ // class, the first is on vertex 0
+ // (==face 0 in some sense), the
+ // second on face 1:
+ return (((shape_index == 0) && (face_index == 0)) ||
+ ((shape_index == 1) && (face_index == 1)));
+}
+
+template <>
+bool
+FE_Q_DG0<1,3>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+ Assert (face_index < GeometryInfo<1>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell));
+
+
+ // in 1d, things are simple. since
+ // there is only one degree of
+ // freedom per vertex in this
+ // class, the first is on vertex 0
+ // (==face 0 in some sense), the
+ // second on face 1:
+ return (((shape_index == 0) && (face_index == 0)) ||
+ ((shape_index == 1) && (face_index == 1)));
+}
+
+
+template <int dim, int spacedim>
+bool
+FE_Q_DG0<dim,spacedim>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+ Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
+
+
+ // first, special-case interior
+ // shape functions, since they
+ // have no support no-where on
+ // the boundary
+ if (((dim==2) && (shape_index>=this->first_quad_index))
+ ||
+ ((dim==3) && (shape_index>=this->first_hex_index)))
+ return false;
+
+ // let's see whether this is a
+ // vertex
+ if (shape_index < this->first_line_index)
+ {
+ // for Q_DG0 elements, there is
+ // one dof per vertex, so
+ // shape_index==vertex_number. check
+ // whether this vertex is
+ // on the given face. thus,
+ // for each face, give a
+ // list of vertices
+ const unsigned int vertex_no = shape_index;
+ Assert (vertex_no < GeometryInfo<dim>::vertices_per_cell,
+ ExcInternalError());
+
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ if (GeometryInfo<dim>::face_to_cell_vertices(face_index, v) == vertex_no)
+ return true;
+
+ return false;
+ }
+ else if (shape_index < this->first_quad_index)
+ // ok, dof is on a line
+ {
+ const unsigned int line_index
+ = (shape_index - this->first_line_index) / this->dofs_per_line;
+ Assert (line_index < GeometryInfo<dim>::lines_per_cell,
+ ExcInternalError());
+
+ // in 2d, the line is the
+ // face, so get the line
+ // index
+ if (dim == 2)
+ return (line_index == face_index);
+ else if (dim == 3)
+ {
+ // see whether the
+ // given line is on the
+ // given face.
+ for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+ if (GeometryInfo<3>::face_to_cell_lines(face_index, l) == line_index)
+ return true;
+
+ return false;
+ }
+ else
+ Assert (false, ExcNotImplemented());
+ }
+ else if (shape_index < this->first_hex_index)
+ // dof is on a quad
+ {
+ const unsigned int quad_index
+ = (shape_index - this->first_quad_index) / this->dofs_per_quad;
+ Assert (static_cast<signed int>(quad_index) <
+ static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
+ ExcInternalError());
+
+ // in 2d, cell bubble are
+ // zero on all faces. but
+ // we have treated this
+ // case above already
+ Assert (dim != 2, ExcInternalError());
+
+ // in 3d,
+ // quad_index=face_index
+ if (dim == 3)
+ return (quad_index == face_index);
+ else
+ Assert (false, ExcNotImplemented());
+ }
+ else
+ // dof on hex
+ {
+ // can only happen in 3d,
+ // but this case has
+ // already been covered
+ // above
+ Assert (false, ExcNotImplemented());
+ return false;
+ }
+
+ // we should not have gotten here
+ Assert (false, ExcInternalError());
+ return false;
+}
+
+
+
+template <int dim, int spacedim>
+std::size_t
+FE_Q_DG0<dim,spacedim>::memory_consumption () const
+{
+ return (FiniteElement<dim,spacedim>::memory_consumption() +
+ MemoryConsumption::memory_consumption (face_index_map));
+}
+
+
+
+// explicit instantiations
+#include "fe_q_dg0.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: fe_q.inst.in 25229 2012-03-09 18:34:59Z pauletti $
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+ {
+#if deal_II_dimension <= deal_II_space_dimension
+ template class FE_Q_DG0<deal_II_dimension, deal_II_space_dimension>;
+#endif
+ }
+
+
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_abf.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_q_hierarchical.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_dgp.h>
--- /dev/null
+//----------------------------------------------------------------------
+// $Id: fe_prolongation_dgq.cc 26378 2012-09-14 12:58:04Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2007, 2008, 2010, 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+//
+// Compute support points
+
+#include "../tests.h"
+#include "fe_prolongation_common.h"
+
+
+
+int
+main()
+{
+ initlog(__FILE__);
+ deallog.threshold_double(1.e-10);
+
+ CHECK_ALL(Q_DG0,1,2);
+ CHECK_ALL(Q_DG0,2,2);
+ CHECK_ALL(Q_DG0,3,2);
+ CHECK_ALL(Q_DG0,4,2);
+
+ CHECK_ALL(Q_DG0,1,3);
+ CHECK_ALL(Q_DG0,2,3);
+}
--- /dev/null
+
+DEAL::Q1<2> constraint
+DEAL::0.5000000 0.5000000
+DEAL::Q1<2> prolongation 0
+DEAL::1.00000 ~ ~ ~
+DEAL::0.500000 0.500000 ~ ~
+DEAL::0.500000 ~ 0.500000 ~
+DEAL::0.250000 0.250000 0.250000 0.250000
+DEAL::Q1<2> prolongation 1
+DEAL::0.500000 0.500000 ~ ~
+DEAL::~ 1.00000 ~ ~
+DEAL::0.250000 0.250000 0.250000 0.250000
+DEAL::~ 0.500000 ~ 0.500000
+DEAL::Q1<2> prolongation 2
+DEAL::0.500000 ~ 0.500000 ~
+DEAL::0.250000 0.250000 0.250000 0.250000
+DEAL::~ ~ 1.00000 ~
+DEAL::~ ~ 0.500000 0.500000
+DEAL::Q1<2> prolongation 3
+DEAL::0.250000 0.250000 0.250000 0.250000
+DEAL::~ 0.500000 ~ 0.500000
+DEAL::~ ~ 0.500000 0.500000
+DEAL::~ ~ ~ 1.00000
+DEAL::Q2<2> constraint
+DEAL::~ ~ 1.000000
+DEAL::0.3750000 -0.1250000 0.7500000
+DEAL::-0.1250000 0.3750000 0.7500000
+DEAL::Q2<2> prolongation 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::0.375000 -0.125000 ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::0.140625 -0.0468750 -0.0468750 0.0156250 0.281250 -0.0937500 0.281250 -0.0937500 0.562500
+DEAL::Q2<2> prolongation 1
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::~ 0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~
+DEAL::-0.125000 0.375000 ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::-0.0468750 0.140625 0.0156250 -0.0468750 -0.0937500 0.281250 0.281250 -0.0937500 0.562500
+DEAL::Q2<2> prolongation 2
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::-0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::~ ~ 0.375000 -0.125000 ~ ~ ~ 0.750000 ~
+DEAL::-0.0468750 0.0156250 0.140625 -0.0468750 0.281250 -0.0937500 -0.0937500 0.281250 0.562500
+DEAL::Q2<2> prolongation 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::~ -0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::~ ~ -0.125000 0.375000 ~ ~ ~ 0.750000 ~
+DEAL::0.0156250 -0.0468750 -0.0468750 0.140625 -0.0937500 0.281250 -0.0937500 0.281250 0.562500
+DEAL::Q3<2> constraint
+DEAL::-0.06250000 -0.06250000 0.5625000 0.5625000
+DEAL::0.3125000 0.06250000 0.9375000 -0.3125000
+DEAL::~ ~ 1.000000 ~
+DEAL::~ ~ ~ 1.000000
+DEAL::0.06250000 0.3125000 -0.3125000 0.9375000
+DEAL::Q3<2> prolongation 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::-0.0625000 -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::-0.0625000 ~ -0.0625000 ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::0.00390625 0.00390625 0.00390625 0.00390625 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 0.316406 0.316406 0.316406 0.316406
+DEAL::0.312500 ~ 0.0625000 ~ 0.937500 -0.312500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::-0.0195312 -0.0195312 -0.00390625 -0.00390625 -0.0585937 0.0195312 -0.0585937 0.0195312 0.175781 0.175781 0.0351562 0.0351562 0.527344 0.527344 -0.175781 -0.175781
+DEAL::~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~
+DEAL::0.312500 0.0625000 ~ ~ ~ ~ ~ ~ 0.937500 -0.312500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::-0.0195312 -0.00390625 -0.0195312 -0.00390625 0.175781 0.175781 0.0351562 0.0351562 -0.0585937 0.0195312 -0.0585937 0.0195312 0.527344 -0.175781 0.527344 -0.175781
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ 0.562500 ~ 0.562500 ~
+DEAL::0.0976562 0.0195312 0.0195312 0.00390625 0.292969 -0.0976562 0.0585937 -0.0195312 0.292969 -0.0976562 0.0585937 -0.0195312 0.878906 -0.292969 -0.292969 0.0976562
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.312500 ~ 0.0625000 ~ 0.937500 ~ -0.312500 ~
+DEAL::~ ~ ~ ~ 0.312500 ~ 0.0625000 ~ ~ ~ ~ ~ 0.937500 -0.312500 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::Q3<2> prolongation 1
+DEAL::-0.0625000 -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::0.00390625 0.00390625 0.00390625 0.00390625 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 0.316406 0.316406 0.316406 0.316406
+DEAL::~ -0.0625000 ~ -0.0625000 ~ ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::-0.0195312 -0.0195312 -0.00390625 -0.00390625 -0.0585937 0.0195312 -0.0585937 0.0195312 0.175781 0.175781 0.0351562 0.0351562 0.527344 0.527344 -0.175781 -0.175781
+DEAL::~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~
+DEAL::~ 0.312500 ~ 0.0625000 ~ ~ 0.937500 -0.312500 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::0.0625000 0.312500 ~ ~ ~ ~ ~ ~ -0.312500 0.937500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ 0.562500 ~ 0.562500
+DEAL::-0.00390625 -0.0195312 -0.00390625 -0.0195312 0.0351562 0.0351562 0.175781 0.175781 0.0195312 -0.0585938 0.0195312 -0.0585937 -0.175781 0.527344 -0.175781 0.527344
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 0.312500 ~ 0.0625000 ~ 0.937500 ~ -0.312500
+DEAL::0.0195312 0.0976562 0.00390625 0.0195312 0.0585937 -0.0195312 0.292969 -0.0976562 -0.0976562 0.292969 -0.0195312 0.0585937 -0.292969 0.878906 0.0976562 -0.292969
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ 0.0625000 ~ 0.312500 ~ ~ ~ ~ ~ -0.312500 0.937500 ~ ~
+DEAL::Q3<2> prolongation 2
+DEAL::-0.0625000 ~ -0.0625000 ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::0.00390625 0.00390625 0.00390625 0.00390625 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 0.316406 0.316406 0.316406 0.316406
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ -0.0625000 -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::0.0625000 ~ 0.312500 ~ -0.312500 0.937500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500
+DEAL::-0.00390625 -0.00390625 -0.0195312 -0.0195312 0.0195312 -0.0585938 0.0195312 -0.0585937 0.0351562 0.0351562 0.175781 0.175781 -0.175781 -0.175781 0.527344 0.527344
+DEAL::-0.0195312 -0.00390625 -0.0195312 -0.00390625 0.175781 0.175781 0.0351562 0.0351562 -0.0585937 0.0195312 -0.0585937 0.0195312 0.527344 -0.175781 0.527344 -0.175781
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ 0.562500 ~ 0.562500 ~
+DEAL::~ ~ 0.312500 0.0625000 ~ ~ ~ ~ ~ ~ 0.937500 -0.312500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 0.312500 ~ 0.0625000 ~ ~ ~ ~ ~ ~ 0.937500 -0.312500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::0.0195312 0.00390625 0.0976562 0.0195312 -0.0976562 0.292969 -0.0195312 0.0585937 0.0585937 -0.0195312 0.292969 -0.0976562 -0.292969 0.0976562 0.878906 -0.292969
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.0625000 ~ 0.312500 ~ -0.312500 ~ 0.937500 ~
+DEAL::Q3<2> prolongation 3
+DEAL::0.00390625 0.00390625 0.00390625 0.00390625 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 -0.0351562 0.316406 0.316406 0.316406 0.316406
+DEAL::~ -0.0625000 ~ -0.0625000 ~ ~ 0.562500 0.562500 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ -0.0625000 -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ ~ ~ ~ ~ ~ 0.562500 0.562500
+DEAL::-0.00390625 -0.00390625 -0.0195312 -0.0195312 0.0195312 -0.0585938 0.0195312 -0.0585937 0.0351562 0.0351562 0.175781 0.175781 -0.175781 -0.175781 0.527344 0.527344
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 0.0625000 ~ 0.312500 ~ ~ -0.312500 0.937500 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0625000 ~ -0.0625000 ~ 0.562500 ~ 0.562500
+DEAL::-0.00390625 -0.0195312 -0.00390625 -0.0195312 0.0351562 0.0351562 0.175781 0.175781 0.0195312 -0.0585938 0.0195312 -0.0585937 -0.175781 0.527344 -0.175781 0.527344
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ 0.0625000 0.312500 ~ ~ ~ ~ ~ ~ -0.312500 0.937500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ 0.0625000 ~ 0.312500 ~ ~ ~ ~ ~ ~ -0.312500 0.937500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 0.0625000 ~ 0.312500 ~ -0.312500 ~ 0.937500
+DEAL::0.00390625 0.0195312 0.0195312 0.0976562 -0.0195312 0.0585938 -0.0976562 0.292969 -0.0195312 0.0585938 -0.0976562 0.292969 0.0976562 -0.292969 -0.292969 0.878906
+DEAL::Q1<3> constraint
+DEAL::0.2500000 0.2500000 0.2500000 0.2500000
+DEAL::0.5000000 ~ 0.5000000 ~
+DEAL::~ 0.5000000 ~ 0.5000000
+DEAL::0.5000000 0.5000000 ~ ~
+DEAL::~ ~ 0.5000000 0.5000000
+DEAL::Q1<3> prolongation 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::0.500000 0.500000 ~ ~ ~ ~ ~ ~
+DEAL::0.500000 ~ 0.500000 ~ ~ ~ ~ ~
+DEAL::0.250000 0.250000 0.250000 0.250000 ~ ~ ~ ~
+DEAL::0.500000 ~ ~ ~ 0.500000 ~ ~ ~
+DEAL::0.250000 0.250000 ~ ~ 0.250000 0.250000 ~ ~
+DEAL::0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000 ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::Q1<3> prolongation 1
+DEAL::0.500000 0.500000 ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::0.250000 0.250000 0.250000 0.250000 ~ ~ ~ ~
+DEAL::~ 0.500000 ~ 0.500000 ~ ~ ~ ~
+DEAL::0.250000 0.250000 ~ ~ 0.250000 0.250000 ~ ~
+DEAL::~ 0.500000 ~ ~ ~ 0.500000 ~ ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ 0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000
+DEAL::Q1<3> prolongation 2
+DEAL::0.500000 ~ 0.500000 ~ ~ ~ ~ ~
+DEAL::0.250000 0.250000 0.250000 0.250000 ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ 0.500000 0.500000 ~ ~ ~ ~
+DEAL::0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000 ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ ~ 0.500000 ~ ~ ~ 0.500000 ~
+DEAL::~ ~ 0.250000 0.250000 ~ ~ 0.250000 0.250000
+DEAL::Q1<3> prolongation 3
+DEAL::0.250000 0.250000 0.250000 0.250000 ~ ~ ~ ~
+DEAL::~ 0.500000 ~ 0.500000 ~ ~ ~ ~
+DEAL::~ ~ 0.500000 0.500000 ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ 0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000
+DEAL::~ ~ 0.250000 0.250000 ~ ~ 0.250000 0.250000
+DEAL::~ ~ ~ 0.500000 ~ ~ ~ 0.500000
+DEAL::Q1<3> prolongation 4
+DEAL::0.500000 ~ ~ ~ 0.500000 ~ ~ ~
+DEAL::0.250000 0.250000 ~ ~ 0.250000 0.250000 ~ ~
+DEAL::0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000 ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ 0.500000 0.500000 ~ ~
+DEAL::~ ~ ~ ~ 0.500000 ~ 0.500000 ~
+DEAL::~ ~ ~ ~ 0.250000 0.250000 0.250000 0.250000
+DEAL::Q1<3> prolongation 5
+DEAL::0.250000 0.250000 ~ ~ 0.250000 0.250000 ~ ~
+DEAL::~ 0.500000 ~ ~ ~ 0.500000 ~ ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ 0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000
+DEAL::~ ~ ~ ~ 0.500000 0.500000 ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ 0.250000 0.250000 0.250000 0.250000
+DEAL::~ ~ ~ ~ ~ 0.500000 ~ 0.500000
+DEAL::Q1<3> prolongation 6
+DEAL::0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000 ~
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ ~ 0.500000 ~ ~ ~ 0.500000 ~
+DEAL::~ ~ 0.250000 0.250000 ~ ~ 0.250000 0.250000
+DEAL::~ ~ ~ ~ 0.500000 ~ 0.500000 ~
+DEAL::~ ~ ~ ~ 0.250000 0.250000 0.250000 0.250000
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ 0.500000 0.500000
+DEAL::Q1<3> prolongation 7
+DEAL::0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::~ 0.250000 ~ 0.250000 ~ 0.250000 ~ 0.250000
+DEAL::~ ~ 0.250000 0.250000 ~ ~ 0.250000 0.250000
+DEAL::~ ~ ~ 0.500000 ~ ~ ~ 0.500000
+DEAL::~ ~ ~ ~ 0.250000 0.250000 0.250000 0.250000
+DEAL::~ ~ ~ ~ ~ 0.500000 ~ 0.500000
+DEAL::~ ~ ~ ~ ~ ~ 0.500000 0.500000
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::Q2<3> constraint
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.000000
+DEAL::~ ~ ~ ~ 1.000000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.000000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.000000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.000000 ~
+DEAL::~ ~ ~ ~ ~ ~ 0.3750000 -0.1250000 0.7500000
+DEAL::~ ~ ~ ~ ~ ~ -0.1250000 0.3750000 0.7500000
+DEAL::~ ~ ~ ~ 0.3750000 -0.1250000 ~ ~ 0.7500000
+DEAL::~ ~ ~ ~ -0.1250000 0.3750000 ~ ~ 0.7500000
+DEAL::0.3750000 ~ -0.1250000 ~ 0.7500000 ~ ~ ~ ~
+DEAL::-0.1250000 ~ 0.3750000 ~ 0.7500000 ~ ~ ~ ~
+DEAL::~ 0.3750000 ~ -0.1250000 ~ 0.7500000 ~ ~ ~
+DEAL::~ -0.1250000 ~ 0.3750000 ~ 0.7500000 ~ ~ ~
+DEAL::0.3750000 -0.1250000 ~ ~ ~ ~ 0.7500000 ~ ~
+DEAL::-0.1250000 0.3750000 ~ ~ ~ ~ 0.7500000 ~ ~
+DEAL::~ ~ 0.3750000 -0.1250000 ~ ~ ~ 0.7500000 ~
+DEAL::~ ~ -0.1250000 0.3750000 ~ ~ ~ 0.7500000 ~
+DEAL::0.1406250 -0.04687500 -0.04687500 0.01562500 0.2812500 -0.09375000 0.2812500 -0.09375000 0.5625000
+DEAL::-0.04687500 0.1406250 0.01562500 -0.04687500 -0.09375000 0.2812500 0.2812500 -0.09375000 0.5625000
+DEAL::-0.04687500 0.01562500 0.1406250 -0.04687500 0.2812500 -0.09375000 -0.09375000 0.2812500 0.5625000
+DEAL::0.01562500 -0.04687500 -0.04687500 0.1406250 -0.09375000 0.2812500 -0.09375000 0.2812500 0.5625000
+DEAL::Q2<3> prolongation 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::0.375000 ~ -0.125000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000
+DEAL::0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::0.140625 ~ -0.0468750 ~ -0.0468750 ~ 0.0156250 ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ -0.0937500 ~ 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 0.281250 -0.0937500 0.562500
+DEAL::0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ 0.281250 -0.0937500 ~ ~ ~ ~ 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 ~ ~ 0.281250 -0.0937500 0.562500
+DEAL::0.140625 -0.0468750 -0.0468750 0.0156250 ~ ~ ~ ~ 0.281250 -0.0937500 0.281250 -0.0937500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 -0.0468750 0.0156250 0.281250 -0.0937500 0.281250 -0.0937500 ~ ~ 0.562500
+DEAL::0.0527344 -0.0175781 -0.0175781 0.00585938 -0.0175781 0.00585938 0.00585938 -0.00195312 0.105469 -0.0351562 0.105469 -0.0351562 -0.0351562 0.0117188 -0.0351562 0.0117188 0.105469 -0.0351562 -0.0351562 0.0117188 0.210938 -0.0703125 0.210938 -0.0703125 0.210938 -0.0703125 0.421875
+DEAL::Q2<3> prolongation 1
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ 0.375000 ~ -0.125000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::-0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 0.281250 -0.0937500 0.562500
+DEAL::~ 0.140625 ~ -0.0468750 ~ -0.0468750 ~ 0.0156250 ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ -0.0937500 ~ 0.562500 ~ ~ ~ ~ ~
+DEAL::-0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ -0.0937500 0.281250 ~ ~ ~ ~ 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 ~ ~ 0.281250 -0.0937500 0.562500
+DEAL::-0.0468750 0.140625 0.0156250 -0.0468750 ~ ~ ~ ~ -0.0937500 0.281250 0.281250 -0.0937500 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 0.0156250 -0.0468750 -0.0937500 0.281250 0.281250 -0.0937500 ~ ~ 0.562500
+DEAL::-0.0175781 0.0527344 0.00585938 -0.0175781 0.00585938 -0.0175781 -0.00195312 0.00585938 -0.0351562 0.105469 0.105469 -0.0351562 0.0117188 -0.0351562 -0.0351562 0.0117188 -0.0351562 0.105469 0.0117188 -0.0351562 -0.0703125 0.210938 0.210938 -0.0703125 0.210938 -0.0703125 0.421875
+DEAL::Q2<3> prolongation 2
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::-0.125000 ~ 0.375000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::-0.0468750 ~ 0.140625 ~ 0.0156250 ~ -0.0468750 ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ ~ -0.0937500 ~ 0.281250 ~ 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 0.281250 -0.0937500 0.562500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 ~ ~ 0.281250 -0.0937500 0.562500
+DEAL::~ ~ 0.140625 -0.0468750 ~ ~ -0.0468750 0.0156250 ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ 0.281250 -0.0937500 ~ ~ ~ 0.562500 ~ ~ ~
+DEAL::-0.0468750 0.0156250 0.140625 -0.0468750 ~ ~ ~ ~ 0.281250 -0.0937500 -0.0937500 0.281250 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 0.140625 -0.0468750 0.281250 -0.0937500 -0.0937500 0.281250 ~ ~ 0.562500
+DEAL::-0.0175781 0.00585938 0.0527344 -0.0175781 0.00585938 -0.00195312 -0.0175781 0.00585938 0.105469 -0.0351562 -0.0351562 0.105469 -0.0351562 0.0117188 0.0117188 -0.0351562 -0.0351562 0.0117188 0.105469 -0.0351562 0.210938 -0.0703125 -0.0703125 0.210938 0.210938 -0.0703125 0.421875
+DEAL::Q2<3> prolongation 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ -0.125000 ~ 0.375000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~
+DEAL::~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ 0.375000 ~ ~ ~ -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 0.281250 -0.0937500 0.562500
+DEAL::~ -0.0468750 ~ 0.140625 ~ 0.0156250 ~ -0.0468750 ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ ~ -0.0937500 ~ 0.281250 ~ 0.562500 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 ~ ~ 0.281250 -0.0937500 0.562500
+DEAL::~ ~ -0.0468750 0.140625 ~ ~ 0.0156250 -0.0468750 ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ ~ -0.0937500 0.281250 ~ ~ ~ 0.562500 ~ ~ ~
+DEAL::0.0156250 -0.0468750 -0.0468750 0.140625 ~ ~ ~ ~ -0.0937500 0.281250 -0.0937500 0.281250 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 -0.0468750 0.140625 -0.0937500 0.281250 -0.0937500 0.281250 ~ ~ 0.562500
+DEAL::0.00585938 -0.0175781 -0.0175781 0.0527344 -0.00195312 0.00585938 0.00585938 -0.0175781 -0.0351562 0.105469 -0.0351562 0.105469 0.0117188 -0.0351562 0.0117188 -0.0351562 0.0117188 -0.0351562 -0.0351562 0.105469 -0.0703125 0.210938 -0.0703125 0.210938 0.210938 -0.0703125 0.421875
+DEAL::Q2<3> prolongation 4
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ 0.375000 ~ -0.125000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::-0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::-0.0468750 ~ 0.0156250 ~ 0.140625 ~ -0.0468750 ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ ~ 0.281250 ~ -0.0937500 ~ 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 -0.0937500 0.281250 0.562500
+DEAL::-0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ 0.281250 -0.0937500 ~ ~ ~ ~ 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 ~ ~ -0.0937500 0.281250 0.562500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.140625 -0.0468750 -0.0468750 0.0156250 0.281250 -0.0937500 0.281250 -0.0937500 ~ ~ 0.562500
+DEAL::~ ~ ~ ~ 0.140625 -0.0468750 -0.0468750 0.0156250 ~ ~ ~ ~ 0.281250 -0.0937500 0.281250 -0.0937500 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~
+DEAL::-0.0175781 0.00585938 0.00585938 -0.00195312 0.0527344 -0.0175781 -0.0175781 0.00585938 -0.0351562 0.0117188 -0.0351562 0.0117188 0.105469 -0.0351562 0.105469 -0.0351562 0.105469 -0.0351562 -0.0351562 0.0117188 0.210938 -0.0703125 0.210938 -0.0703125 -0.0703125 0.210938 0.421875
+DEAL::Q2<3> prolongation 5
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 ~ -0.125000 ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ 0.375000 ~ -0.125000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~
+DEAL::~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 -0.0937500 0.281250 0.562500
+DEAL::~ -0.0468750 ~ 0.0156250 ~ 0.140625 ~ -0.0468750 ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ ~ 0.281250 ~ -0.0937500 ~ 0.562500 ~ ~ ~ ~ ~
+DEAL::0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ -0.0937500 0.281250 ~ ~ ~ ~ 0.562500 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 ~ ~ -0.0937500 0.281250 0.562500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.140625 0.0156250 -0.0468750 -0.0937500 0.281250 0.281250 -0.0937500 ~ ~ 0.562500
+DEAL::~ ~ ~ ~ -0.0468750 0.140625 0.0156250 -0.0468750 ~ ~ ~ ~ -0.0937500 0.281250 0.281250 -0.0937500 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~
+DEAL::0.00585938 -0.0175781 -0.00195312 0.00585938 -0.0175781 0.0527344 0.00585938 -0.0175781 0.0117188 -0.0351562 -0.0351562 0.0117188 -0.0351562 0.105469 0.105469 -0.0351562 -0.0351562 0.105469 0.0117188 -0.0351562 -0.0703125 0.210938 0.210938 -0.0703125 -0.0703125 0.210938 0.421875
+DEAL::Q2<3> prolongation 6
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ ~ -0.125000 ~ 0.375000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ ~ 0.375000 -0.125000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::0.0156250 ~ -0.0468750 ~ -0.0468750 ~ 0.140625 ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ 0.281250 ~ 0.562500 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 -0.0937500 0.281250 0.562500
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ ~ ~ ~ 0.281250 -0.0937500 ~ ~ -0.0937500 0.281250 0.562500
+DEAL::~ ~ -0.0468750 0.0156250 ~ ~ 0.140625 -0.0468750 ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ 0.281250 -0.0937500 ~ ~ ~ 0.562500 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.0468750 0.0156250 0.140625 -0.0468750 0.281250 -0.0937500 -0.0937500 0.281250 ~ ~ 0.562500
+DEAL::~ ~ ~ ~ -0.0468750 0.0156250 0.140625 -0.0468750 ~ ~ ~ ~ 0.281250 -0.0937500 -0.0937500 0.281250 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~
+DEAL::0.00585938 -0.00195312 -0.0175781 0.00585938 -0.0175781 0.00585938 0.0527344 -0.0175781 -0.0351562 0.0117188 0.0117188 -0.0351562 0.105469 -0.0351562 -0.0351562 0.105469 -0.0351562 0.0117188 0.105469 -0.0351562 0.210938 -0.0703125 -0.0703125 0.210938 -0.0703125 0.210938 0.421875
+DEAL::Q2<3> prolongation 7
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ 0.375000 ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ -0.125000 ~ 0.375000 ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~
+DEAL::~ ~ ~ ~ ~ ~ -0.125000 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 0.375000 0.750000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~
+DEAL::~ ~ ~ -0.125000 ~ ~ ~ 0.375000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.750000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 -0.0937500 0.281250 0.562500
+DEAL::~ 0.0156250 ~ -0.0468750 ~ -0.0468750 ~ 0.140625 ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ ~ -0.0937500 ~ 0.281250 ~ 0.562500 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ ~ ~ ~ -0.0937500 0.281250 ~ ~ -0.0937500 0.281250 0.562500
+DEAL::~ ~ 0.0156250 -0.0468750 ~ ~ -0.0468750 0.140625 ~ ~ ~ -0.0937500 ~ ~ ~ 0.281250 ~ ~ -0.0937500 0.281250 ~ ~ ~ 0.562500 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.0156250 -0.0468750 -0.0468750 0.140625 -0.0937500 0.281250 -0.0937500 0.281250 ~ ~ 0.562500
+DEAL::~ ~ ~ ~ 0.0156250 -0.0468750 -0.0468750 0.140625 ~ ~ ~ ~ -0.0937500 0.281250 -0.0937500 0.281250 ~ ~ ~ ~ ~ ~ ~ ~ ~ 0.562500 ~
+DEAL::-0.00195312 0.00585938 0.00585938 -0.0175781 0.00585938 -0.0175781 -0.0175781 0.0527344 0.0117188 -0.0351562 0.0117188 -0.0351562 -0.0351562 0.105469 -0.0351562 0.105469 0.0117188 -0.0351562 -0.0351562 0.105469 -0.0703125 0.210938 -0.0703125 0.210938 -0.0703125 0.210938 0.421875
--- /dev/null
+// Version:
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+
+
+// test FE_Q_DG0 (modified step-22)
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_ilu.h>
+
+#include <deal.II/grid/tria.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/grid_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <fstream>
+#include <sstream>
+
+namespace Step22
+{
+ using namespace dealii;
+
+ template <int dim>
+ struct InnerPreconditioner;
+
+ template <>
+ struct InnerPreconditioner<2>
+ {
+ typedef SparseDirectUMFPACK type;
+ };
+
+ template <>
+ struct InnerPreconditioner<3>
+ {
+ typedef SparseILU<double> type;
+ };
+
+ template <int dim>
+ class StokesProblem
+ {
+ public:
+ StokesProblem (const unsigned int degree, FESystem<dim> & fe_);
+ void run ();
+
+ private:
+ void setup_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results (const unsigned int refinement_cycle);
+
+ void divergence_velocity(const BlockVector<double> &calc_solution,
+ Vector<double> &output_vector,
+ const Quadrature<dim> &quadrature, bool norm);
+
+ const unsigned int degree;
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> &fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+
+ std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>
+ A_preconditioner;
+
+ ConvergenceTable convergence_table;
+ };
+
+ template <int dim>
+ class ExactSolution : public Function<dim>
+ {
+ public:
+ ExactSolution () : Function<dim>(dim+1) {}
+
+ /*virtual*/ double value (const Point<dim> &p,
+ const unsigned int component) const;
+
+ /*virtual*/
+ Tensor<1,dim> gradient (const Point< dim> &p,
+ const unsigned int component) const;
+
+ /*virtual*/
+ double laplacian (const Point<dim> &p,
+ const unsigned int component) const;
+ };
+
+
+ template <int dim>
+ double
+ ExactSolution<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+
+ double x = p[0];
+ double y = p[1];
+
+ switch (component)
+ {
+ //velocity
+ case 0:
+ return 2*(x-1)*(x-1)*x*x*(y-1)*y*(2*y-1);
+ break;
+ case 1:
+ return -2*(y-1)*(y-1)*y*y*(x-1)*x*(2*x-1);
+ break;
+ //pressure
+ case 2:
+ //discontinuous Boffi
+ return y*(1-y)*exp((x-.5)*(x-.5))-.5+(x<.5);
+ //discontinuous simple
+ //return (x<.5)-.5;
+ //continuous
+ //return .5*x*x-1/6;
+
+ }
+ ExcNotImplemented();
+
+ return 0;
+ }
+
+ template <int dim>
+ Tensor<1,dim>
+ ExactSolution<dim>::gradient (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+
+ double x = p[0];
+ double y = p[1];
+
+ Tensor<1,dim> gradient;
+
+ switch (component)
+ {
+ //velocity
+ case 0:
+ gradient[0]=x*(x*(x*y*(y*(16*y-24)+8)+y*((36-24*y)*y-12))
+ +y*(y*(8*y-12)+4));
+ gradient[1]=x*x*(x*(x*(y*(12*y-12)+2)+(24-24*y)*y-4)+y*(12*y-12)+2);
+ break;
+ case 1:
+ gradient[0]=x*(x*((24-12*y)*y-12)*y*y+(y*(12*y-24)+12)*y*y)
+ +((4-2*y)*y-2)*y*y;
+ gradient[1]=x*(x*(x*y*((24-16*y)*y-8)+y*(y*(24*y-36)+12))
+ +y*((12-8*y)*y-4));
+ break;
+ //pressure
+ case 2:
+ //discontinuous Boffi
+ gradient[0]=-exp((x-.5)*(x-.5))*(2*x-1)*(y-1)*y;
+ gradient[1]=-exp((x-.5)*(x-.5))*(2*y-1);
+ //discontinuous simple
+ //gradient[0]=0;
+ //gradient[1]=0;
+ //continuous
+ //gradient[0]=x;
+ //gradient[1]=0;
+ }
+
+ return gradient;
+ }
+
+
+ template <int dim>
+ double
+ ExactSolution<dim>::laplacian (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+
+ double x = p[0];
+ double y = p[1];
+
+ switch (component)
+ {
+ //velocity
+ case 0:
+ return x*(x*(x*(x*(24*y-12)-48*y+24)
+ +y*(y*(48*y-72)+48)-12)
+ +y*((72-48*y)*y-24))
+ +y*(y*(8*y-12)+4);
+
+ case 1:
+ return x*(x*(x*((48-48*y)*y-8)
+ +y*(72*y-72)+12)
+ +y*(y*((48-24*y)*y-48)+24)-4)
+ +(y*(12*y-24)+12)*y*y;
+ }
+ ExcNotImplemented();
+ return 0;
+ }
+
+ template <int dim>
+ class JumpFunction : public Function<dim>
+ {
+ public:
+ JumpFunction () : Function<dim>(1) {}
+
+ double jump (const Point< dim> &p,
+ const Point<dim> &normal) const;
+ };
+
+ template <int dim>
+ double
+ JumpFunction<dim>::jump (const Point<dim> &p,
+ const Point<dim> &normal) const
+ {
+ double x = p[0];
+ double y = p[1];
+ //discontinuous
+ if (std::abs(x-.5)>1e-10)
+ return 0;
+ if (normal[0]>0)
+ return -1;
+ return 1;
+ //continuous
+ //return 0;
+ }
+
+ template <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>(dim+1) {}
+
+ /*virtual*/ double value (const Point<dim> &p,
+ const unsigned int component) const;
+
+ const ExactSolution<dim> solution;
+ };
+
+ template <int dim>
+ double
+ RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+ if (component==dim)
+ return 0;
+ //grad p -laplace u
+ return solution.gradient(p,dim)[component]
+ -solution.laplacian(p,component);
+ }
+
+
+ template <class Matrix, class Preconditioner>
+ class InverseMatrix : public Subscriptor
+ {
+ public:
+ InverseMatrix (const Matrix &m,
+ const Preconditioner &preconditioner);
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ private:
+ const SmartPointer<const Matrix> matrix;
+ const SmartPointer<const Preconditioner> preconditioner;
+ };
+
+
+ template <class Matrix, class Preconditioner>
+ InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
+ const Preconditioner &preconditioner)
+ :
+ matrix (&m),
+ preconditioner (&preconditioner)
+ {}
+
+ template <class Matrix, class Preconditioner>
+ void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
+ SolverCG<> cg (solver_control);
+
+ dst = 0;
+
+ cg.solve (*matrix, dst, src, *preconditioner);
+
+ /*std::cout << " "
+ << solver_control.last_step()
+ << " inner CG iterations for pressure"
+ << std::endl;*/
+ }
+
+ template <class Preconditioner>
+ class SchurComplement : public Subscriptor
+ {
+ public:
+ SchurComplement (const BlockSparseMatrix<double> &system_matrix,
+ const InverseMatrix<SparseMatrix<double>, Preconditioner>
+ &A_inverse);
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ private:
+ const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+ const SmartPointer<const InverseMatrix<SparseMatrix<double>,
+ Preconditioner> > A_inverse;
+
+ mutable Vector<double> tmp1, tmp2;
+ };
+
+
+
+ template <class Preconditioner>
+ SchurComplement<Preconditioner>::
+ SchurComplement (const BlockSparseMatrix<double> &system_matrix,
+ const InverseMatrix<SparseMatrix<double>,Preconditioner>
+ &A_inverse)
+ :
+ system_matrix (&system_matrix),
+ A_inverse (&A_inverse),
+ tmp1 (system_matrix.block(0,0).m()),
+ tmp2 (system_matrix.block(0,0).m())
+ {}
+
+
+ template <class Preconditioner>
+ void SchurComplement<Preconditioner>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ system_matrix->block(0,1).vmult (tmp1, src);
+ A_inverse->vmult (tmp2, tmp1);
+ system_matrix->block(1,0).vmult (dst, tmp2);
+ }
+
+ template <int dim>
+ StokesProblem<dim>::StokesProblem (const unsigned int degree, FESystem<dim> & fe_)
+ :
+ degree (degree),
+ triangulation (Triangulation<dim>::maximum_smoothing),
+ fe (fe_),
+ dof_handler (triangulation)
+ {}
+
+ template <int dim>
+ void StokesProblem<dim>::setup_dofs ()
+ {
+ A_preconditioner.reset ();
+ system_matrix.clear ();
+
+ dof_handler.distribute_dofs (fe);
+ DoFRenumbering::Cuthill_McKee (dof_handler);
+
+ std::vector<unsigned int> block_component (dim+1,0);
+ block_component[dim] = 1;
+ DoFRenumbering::component_wise (dof_handler, block_component);
+
+ {
+ constraints.clear ();
+ std::vector<bool> component_mask (dim+1, true);
+ component_mask[dim] = false;
+
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ExactSolution<dim>(),
+ constraints,
+ component_mask);
+
+ /*std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
+
+ std::vector<bool>boundary_mask (dim+1, false);
+ boundary_mask[dim]=true;
+
+ DoFTools::extract_boundary_dofs (dof_handler,boundary_mask,boundary_dofs);
+
+ const unsigned int first_boundary_dof
+ = std::distance (boundary_dofs.begin(),std::find
+ (boundary_dofs.begin(),boundary_dofs.end(),true));
+
+ constraints.add_line (first_boundary_dof);
+ for (unsigned int i=first_boundary_dof+1; i<dof_handler.n_dofs();++i)
+ if (boundary_dofs[i] == true)
+ constraints.add_entry (first_boundary_dof, i, -1.);*/
+
+ /*const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ cell->get_dof_indices(local_dof_indices);
+ unsigned int first_disc_dof=local_dof_indices[dofs_per_cell-1];
+ constraints.add_line (first_disc_dof);
+
+ for (++cell; cell!=endc; ++cell)
+ {
+ cell->get_dof_indices(local_dof_indices);
+ if(cell->at_boundary())
+ constraints.add_entry (first_disc_dof,
+ local_dof_indices[dofs_per_cell-1],-1.);
+ }*/
+ }
+
+ constraints.close ();
+
+ std::vector<unsigned int> dofs_per_block (2);
+ DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
+ block_component);
+ const unsigned int n_u = dofs_per_block[0],
+ n_p = dofs_per_block[1];
+
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << " (" << n_u << '+' << n_p << ')'
+ << std::endl;
+
+ {
+ BlockCompressedSimpleSparsityPattern csp (2,2);
+
+ csp.block(0,0).reinit (n_u, n_u);
+ csp.block(1,0).reinit (n_p, n_u);
+ csp.block(0,1).reinit (n_u, n_p);
+ csp.block(1,1).reinit (n_p, n_p);
+
+ csp.collect_sizes();
+
+ DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+ sparsity_pattern.copy_from (csp);
+ }
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (2);
+ solution.block(0).reinit (n_u);
+ solution.block(1).reinit (n_p);
+ solution.collect_sizes ();
+
+ system_rhs.reinit (2);
+ system_rhs.block(0).reinit (n_u);
+ system_rhs.block(1).reinit (n_p);
+ system_rhs.collect_sizes ();
+ }
+
+ template <int dim>
+ void StokesProblem<dim>::assemble_system ()
+ {
+ system_matrix=0;
+ system_rhs=0;
+
+ QGauss<dim> quadrature_formula(degree+2);
+ QGauss<dim-1> quadrature_face(degree+2);
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values |
+ update_gradients);
+
+ FEFaceValues<dim> fe_v_face (fe, quadrature_face,
+ update_values | update_quadrature_points |
+ update_JxW_values | update_normal_vectors);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ const unsigned int n_q_points = quadrature_formula.size();
+ const unsigned int n_q_face = quadrature_face.size();
+
+ FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ const RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim+1));
+
+ const JumpFunction<dim> jumpfunction;
+
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+
+ std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
+ std::vector<double> div_phi_u (dofs_per_cell);
+ std::vector<double> phi_p (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+ local_matrix = 0;
+ local_rhs = 0;
+
+ right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+ rhs_values);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k,
+ q);
+ div_phi_u[k] = fe_values[velocities].divergence (k, q);
+ phi_p[k] = fe_values[pressure].value (k, q);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<=i; ++j)
+ {
+ local_matrix(i,j) += (2*symgrad_phi_u[i] * symgrad_phi_u[j]
+ - div_phi_u[i] * phi_p[j]
+ - phi_p[i] * div_phi_u[j]
+ + phi_p[i] * phi_p[j])
+ * fe_values.JxW(q);
+
+ }
+
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ local_rhs(i) += fe_values.shape_value(i,q) *
+ rhs_values[q](component_i) *
+ fe_values.JxW(q);
+ }
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+ local_matrix(i,j) = local_matrix(j,i);
+
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ typename DoFHandler<dim>::face_iterator face= cell->face(face_no);
+ if (face->at_boundary()==false)
+ {
+ typename DoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);
+ if (neighbor->index() > cell->index())
+ {
+ fe_v_face.reinit (cell, face_no);
+
+ const std::vector<Point<dim> > &normals =
+ fe_v_face.get_normal_vectors ();
+ const std::vector<Point<dim> > &quad_points=
+ fe_v_face.get_quadrature_points();
+
+ for (unsigned int q=0; q<n_q_face; ++q)
+ {
+ double jump=jumpfunction.jump(quad_points[q],normals[q]);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ if (component_i<dim)
+ local_rhs(i)+=fe_v_face.shape_value(i,q)
+ *jump
+ *normals[q][component_i]
+ *fe_v_face.JxW(q);
+ }
+ }
+ }
+ }
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global (local_matrix, local_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+
+ std::cout << " Computing preconditioner..." << std::endl << std::flush;
+
+ A_preconditioner
+ = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new
+ typename InnerPreconditioner<dim>::type());
+ A_preconditioner->initialize (system_matrix.block(0,0),
+ typename
+ InnerPreconditioner<dim>::type::AdditionalData());
+
+ }
+
+ template <int dim>
+ void StokesProblem<dim>::solve ()
+ {
+ const InverseMatrix<SparseMatrix<double>,
+ typename InnerPreconditioner<dim>::type>
+ A_inverse (system_matrix.block(0,0), *A_preconditioner);
+ Vector<double> tmp (solution.block(0).size());
+
+ {
+ Vector<double> schur_rhs (solution.block(1).size());
+ A_inverse.vmult (tmp, system_rhs.block(0));
+ system_matrix.block(1,0).vmult (schur_rhs, tmp);
+ schur_rhs -= system_rhs.block(1);
+
+ SchurComplement<typename InnerPreconditioner<dim>::type>
+ schur_complement (system_matrix, A_inverse);
+
+ SolverControl solver_control (solution.block(1).size(),
+ 1e-8*schur_rhs.l2_norm());
+ SolverCG<> cg (solver_control);
+
+ SparseILU<double> preconditioner;
+ preconditioner.initialize (system_matrix.block(1,1),
+ SparseILU<double>::AdditionalData());
+
+ InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+ m_inverse (system_matrix.block(1,1), preconditioner);
+
+ cg.solve (schur_complement, solution.block(1), schur_rhs,
+ m_inverse);
+
+ constraints.distribute (solution);
+
+ std::cout << " "
+ << solver_control.last_step()
+ << " outer CG Schur complement iterations for pressure"
+ << std::endl;
+
+ }
+
+ {
+ system_matrix.block(0,1).vmult (tmp, solution.block(1));
+ tmp *= -1;
+ tmp += system_rhs.block(0);
+
+ A_inverse.vmult (solution.block(0), tmp);
+
+ constraints.distribute (solution);
+ }
+ }
+
+ template <int dim>
+ void
+ StokesProblem<dim>::output_results (const unsigned int refinement_cycle)
+ {
+ const ComponentSelectFunction<dim> pressure_mask (dim, dim+1);
+ const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
+ dim+1);
+
+ ExactSolution<dim> exactsolution;
+
+ const unsigned int n_active_cells=triangulation.n_active_cells();
+
+ Vector<double> difference_per_cell (n_active_cells);
+
+ QGauss<dim> quadrature (degree+3);
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::L2_norm,&velocity_mask);
+ const double L2_error_velocity = difference_per_cell.l2_norm();
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::H1_seminorm,&velocity_mask);
+ const double H1_error_velocity = difference_per_cell.l2_norm();
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::Linfty_norm,&velocity_mask);
+ const double Linfty_error_velocity = difference_per_cell.linfty_norm();
+
+ divergence_velocity(solution, difference_per_cell, quadrature, true);
+
+ const double L2_div_velocity =
+ sqrt(difference_per_cell.mean_value()*n_active_cells);
+
+ divergence_velocity(solution, difference_per_cell, quadrature, false);
+ std::cout<<"maximum divergence per cell: "
+ <<difference_per_cell.linfty_norm()<<std::endl;
+
+ //int_\Omega (f(x)-c)^2 dx is minimized for c=1/|\Omega|\int_\Omega f(x)dx.
+ //That gives \int_\Omega f(x)^2 dx - 1/|\Omega|(\int_\Omega f(x) dx)^2
+ //=\int_\Omega f(x)^2 dx - |\Omega|c^2=\int_\Omega f(x)^2 dx - c^2
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::L2_norm, &pressure_mask);
+ /*std::cout<<"Maximal difference per cell:"
+ <<difference_per_cell.linfty_norm()<<std::endl;*/
+
+ double L2_error_pressure = difference_per_cell.l2_norm() *
+ difference_per_cell.l2_norm();
+
+ std::cout<<"l2 difference "<<L2_error_pressure<<std::endl;
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::mean,&pressure_mask);
+ double integral = difference_per_cell.mean_value()*n_active_cells;
+ std::cout<<"mean difference "<<integral *integral<<std::endl;
+ L2_error_pressure -= integral*integral;
+
+ std::cout<<"Pressure error squared: "<<L2_error_pressure<<std::endl;
+ L2_error_pressure = (L2_error_pressure>0)?sqrt(L2_error_pressure):0;
+
+ /*Vector<double> difference(dim+1);
+ const Point<dim> point(.25,0);
+ VectorTools::point_difference (dof_handler, solution, exactsolution,
+ difference, point);
+ std::cout<<"difference at 0.25,0 "<<difference[dim]<<std::endl;*/
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::H1_seminorm,&pressure_mask);
+ const double H1_error_pressure = difference_per_cell.l2_norm();
+
+ VectorTools::integrate_difference (dof_handler, solution, exactsolution,
+ difference_per_cell, quadrature,
+ VectorTools::Linfty_norm,&pressure_mask);
+ const double Linfty_error_pressure = difference_per_cell.linfty_norm();
+
+ convergence_table.add_value("cycle", refinement_cycle);
+ convergence_table.add_value("cells", n_active_cells);
+ convergence_table.add_value("dofs", dof_handler.n_dofs());
+ convergence_table.add_value("L2_v", L2_error_velocity);
+ convergence_table.add_value("H1_v", H1_error_velocity);
+ convergence_table.add_value("Linfty_v", Linfty_error_velocity);
+ convergence_table.add_value("L2_div_v", L2_div_velocity);
+ convergence_table.add_value("L2_p", L2_error_pressure);
+ convergence_table.add_value("H1_p", H1_error_pressure);
+ convergence_table.add_value("Linfty_p", Linfty_error_pressure);
+
+ std::vector<std::string> solution_names (dim, "velocity");
+ solution_names.push_back ("pressure");
+
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation
+ (dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_component_interpretation
+ .push_back (DataComponentInterpretation::component_is_scalar);
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, solution_names,
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ data_out.build_patches ();
+
+ }
+
+ template <int dim>
+ void StokesProblem<dim>::run ()
+ {
+ Assert(dim==2, ExcNotImplemented());
+ GridGenerator::hyper_cube(triangulation);
+ triangulation.refine_global (1);
+
+ for (unsigned int refinement_cycle = 0; refinement_cycle<5;
+ ++refinement_cycle)
+ {
+ std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+
+ if (refinement_cycle > 0)
+ triangulation.refine_global(1);
+
+ setup_dofs ();
+
+ std::cout << " Assembling..." << std::endl << std::flush;
+ assemble_system ();
+
+ std::cout << " Solving..." << std::flush;
+ solve ();
+
+ output_results (refinement_cycle);
+
+ std::cout << std::endl;
+ }
+
+ convergence_table.set_precision("L2_v", 3);
+ convergence_table.set_precision("H1_v", 3);
+ convergence_table.set_precision("Linfty_v", 3);
+ convergence_table.set_precision("L2_div_v", 3);
+ convergence_table.set_precision("L2_p", 3);
+ convergence_table.set_precision("H1_p", 3);
+ convergence_table.set_precision("Linfty_p", 3);
+
+ convergence_table.set_scientific("L2_v", true);
+ convergence_table.set_scientific("H1_v", true);
+ convergence_table.set_scientific("Linfty_v", true);
+ convergence_table.set_scientific("L2_div_v", true);
+ convergence_table.set_scientific("L2_p", true);
+ convergence_table.set_scientific("H1_p", true);
+ convergence_table.set_scientific("Linfty_p", true);
+
+ convergence_table.set_tex_caption("cells", "\\# cells");
+ convergence_table.set_tex_caption("dofs", "\\# dofs");
+ convergence_table.set_tex_caption("L2_v", "$L^2$-error velocity");
+ convergence_table.set_tex_caption("H1_v", "$H^1$-error velocity");
+ convergence_table.set_tex_caption("Linfty_v",
+ "$L^\\infty$-error velocity");
+ convergence_table.set_tex_caption("L2_div_v",
+ "$L^2$-error divergence velocity");
+ convergence_table.set_tex_caption("L2_p", "$L^2$-error pressure");
+ convergence_table.set_tex_caption("H1_p", "$H^1$-error pressure");
+ convergence_table.set_tex_caption("Linfty_p",
+ "$L^\\infty$-error pressure");
+
+ convergence_table.set_tex_format("cells", "r");
+ convergence_table.set_tex_format("dofs", "r");
+
+ //std::cout << std::endl;
+ //convergence_table.write_text(std::cout);
+
+
+ std::ostringstream error_filename;
+ error_filename << "error-"
+ << Utilities::int_to_string (degree, 1)
+ << ".tex";
+ std::ofstream error_table_file(error_filename.str().c_str());
+ convergence_table.write_tex(error_table_file);
+
+
+ convergence_table.add_column_to_supercolumn("cycle", "n cells");
+ convergence_table.add_column_to_supercolumn("cells", "n cells");
+
+ std::vector<std::string> new_order;
+ new_order.push_back("n cells");
+ new_order.push_back("L2_v");
+ new_order.push_back("H1_v");
+ new_order.push_back("L2_div_v");
+ new_order.push_back("L2_p");
+ new_order.push_back("H1_p");
+ convergence_table.set_column_order (new_order);
+
+ convergence_table.evaluate_convergence_rates
+ ("L2_v",ConvergenceTable::reduction_rate_log2);
+ convergence_table.evaluate_convergence_rates
+ ("H1_v",ConvergenceTable::reduction_rate_log2);
+ convergence_table.evaluate_convergence_rates
+ ("L2_div_v",ConvergenceTable::reduction_rate_log2);
+ convergence_table.evaluate_convergence_rates
+ ("L2_p",ConvergenceTable::reduction_rate_log2);
+ convergence_table.evaluate_convergence_rates
+ ("H1_p",ConvergenceTable::reduction_rate_log2);
+
+ convergence_table.write_text(deallog.get_file_stream());
+ }
+
+ //squared l^2 norm of divergence of velocity
+ template<int dim>
+ void StokesProblem<dim>::divergence_velocity
+ (const BlockVector<double> &calc_solution,
+ Vector<double> &output_vector,
+ const Quadrature<dim> &quadrature, bool norm)
+ {
+ output_vector = 0;
+
+ FEValues<dim> fe_v(fe, quadrature, update_gradients | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature.size();
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ for (; cell!=endc; ++cell)
+ {
+ fe_v.reinit(cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ double div = 0;
+ for (unsigned int i=0; i < dofs_per_cell-1; ++i)
+ {
+ double tmp=0;
+ for (unsigned int d = 0; d < dim; ++d)
+ tmp+= fe_v.shape_grad_component(i,q,d)[d];
+
+ div+=tmp*calc_solution(local_dof_indices[i]);
+ }
+ if (norm)
+ output_vector(cell->index()) += div * div * fe_v.JxW(q); //L^2-Norm
+ else
+ output_vector(cell->index()) += div * fe_v.JxW(q); //Integral
+ }
+ }
+ }
+}
+
+int main ()
+{
+ using namespace dealii;
+ using namespace Step22;
+
+ initlog(__FILE__);
+
+ deallog.depth_file (1);
+
+ unsigned int degree;
+ const unsigned int dim=2;
+ {
+ degree=1;
+ FESystem<dim> fe(FE_Q<dim>(degree+1), dim, FE_Q<dim>(degree), 1);
+ deallog << fe.get_name() << ":" << std::endl;
+ StokesProblem<2> flow_problem(degree, fe);
+ flow_problem.run ();
+ }
+ {
+ degree=1;
+ FESystem<2> fe(FE_Q<dim>(degree+1), dim, FE_Q_DG0<dim>(degree), 1);
+ deallog << fe.get_name() << ":" << std::endl;
+ StokesProblem<2> flow_problem(degree, fe);
+ flow_problem.run ();
+ }
+ {
+ degree=2;
+ FESystem<2> fe(FE_Q<dim>(degree+1), dim, FE_Q_DG0<dim>(degree), 1);
+ deallog << fe.get_name() << ":" << std::endl;
+ StokesProblem<2> flow_problem(degree, fe);
+ flow_problem.run ();
+ }
+
+ return 0;
+}
--- /dev/null
+
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(1)]:
+n cells L2_v H1_v L2_div_v L2_p H1_p
+0 4 1.026e-02 - 1.188e-01 - 1.134e-01 - 2.512e-01 - 1.522e+00 -
+1 16 4.852e-03 1.08 9.199e-02 0.37 9.664e-02 0.23 1.891e-01 0.41 1.904e+00 -0.32
+2 64 1.820e-03 1.41 6.625e-02 0.47 9.883e-02 -0.03 1.343e-01 0.49 2.623e+00 -0.46
+3 256 6.537e-04 1.48 4.717e-02 0.49 1.556e-01 -0.65 9.498e-02 0.50 3.715e+00 -0.50
+4 1024 2.329e-04 1.49 3.347e-02 0.50 2.986e-01 -0.94 6.716e-02 0.50 5.259e+00 -0.50
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q_DG0<2>(1)]:
+n cells L2_v H1_v L2_div_v L2_p H1_p
+0 4 1.267e-03 - 1.783e-02 - 1.017e-02 - 2.085e-02 - 3.231e-01 -
+1 16 1.714e-04 2.89 4.500e-03 1.99 2.991e-03 1.77 5.231e-03 2.00 1.616e-01 1.00
+2 64 2.151e-05 2.99 1.118e-03 2.01 7.756e-04 1.95 1.304e-03 2.00 8.076e-02 1.00
+3 256 2.687e-06 3.00 2.787e-04 2.00 1.960e-04 1.98 3.258e-04 2.00 4.038e-02 1.00
+4 1024 3.357e-07 3.00 6.962e-05 2.00 4.916e-05 2.00 8.146e-05 2.00 2.019e-02 1.00
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q_DG0<2>(2)]:
+n cells L2_v H1_v L2_div_v L2_p H1_p
+0 4 1.002e-04 - 1.968e-03 - 1.795e-03 - 3.744e-04 - 7.466e-03 -
+1 16 6.154e-06 4.03 2.351e-04 3.07 2.281e-04 2.98 3.963e-05 3.24 1.990e-03 1.91
+2 64 3.815e-07 4.01 2.898e-05 3.02 2.872e-05 2.99 5.187e-06 2.93 5.144e-04 1.95
+3 256 2.378e-08 4.00 3.609e-06 3.01 3.600e-06 3.00 6.610e-07 2.97 1.295e-04 1.99
+4 1024 1.486e-09 4.00 4.506e-07 3.00 4.503e-07 3.00 8.372e-08 2.98 3.244e-05 2.00
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_q_hierarchical.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_dgp.h>
--- /dev/null
+//----------------------------------------------------------------------
+// $Id: fe_restriction_q_dg0.cc 26377 2012-09-14 12:53:49Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2007, 2008, 2010, 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+//
+// Compute support points
+
+#include "../tests.h"
+#include "fe_restriction_common.h"
+
+
+
+int
+main()
+{
+ initlog(__FILE__);
+ deallog.threshold_double(1.e-10);
+
+ CHECK_ALL(Q_DG0,1,2);
+ CHECK_ALL(Q_DG0,2,2);
+ CHECK_ALL(Q_DG0,3,2);
+ CHECK_ALL(Q_DG0,4,2);
+
+ CHECK_ALL(Q_DG0,1,3);
+ CHECK_ALL(Q_DG0,2,3);
+}
--- /dev/null
+
+DEAL::Q1<2> constraint
+DEAL::0.5000000 0.5000000
+DEAL::Q1<2> restriction 0
+DEAL::1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::Q1<2> restriction 1
+DEAL::~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::Q1<2> restriction 2
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~
+DEAL::Q1<2> restriction 3
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000
+DEAL::Q2<2> constraint
+DEAL::~ ~ 1.000000
+DEAL::0.3750000 -0.1250000 0.7500000
+DEAL::-0.1250000 0.3750000 0.7500000
+DEAL::Q2<2> restriction 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::Q2<2> restriction 1
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::Q2<2> restriction 2
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<2> restriction 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q3<2> constraint
+DEAL::-0.06250000 -0.06250000 0.5625000 0.5625000
+DEAL::0.3125000 0.06250000 0.9375000 -0.3125000
+DEAL::~ ~ 1.000000 ~
+DEAL::~ ~ ~ 1.000000
+DEAL::0.06250000 0.3125000 -0.3125000 0.9375000
+DEAL::Q3<2> restriction 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q3<2> restriction 1
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q3<2> restriction 2
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q3<2> restriction 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::Q1<3> constraint
+DEAL::0.2500000 0.2500000 0.2500000 0.2500000
+DEAL::0.5000000 ~ 0.5000000 ~
+DEAL::~ 0.5000000 ~ 0.5000000
+DEAL::0.5000000 0.5000000 ~ ~
+DEAL::~ ~ 0.5000000 0.5000000
+DEAL::Q1<3> restriction 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 1
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 2
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 4
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 5
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 6
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q1<3> restriction 7
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000
+DEAL::Q2<3> constraint
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ 1.000000
+DEAL::~ ~ ~ ~ 1.000000 ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.000000 ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.000000 ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.000000 ~
+DEAL::~ ~ ~ ~ ~ ~ 0.3750000 -0.1250000 0.7500000
+DEAL::~ ~ ~ ~ ~ ~ -0.1250000 0.3750000 0.7500000
+DEAL::~ ~ ~ ~ 0.3750000 -0.1250000 ~ ~ 0.7500000
+DEAL::~ ~ ~ ~ -0.1250000 0.3750000 ~ ~ 0.7500000
+DEAL::0.3750000 ~ -0.1250000 ~ 0.7500000 ~ ~ ~ ~
+DEAL::-0.1250000 ~ 0.3750000 ~ 0.7500000 ~ ~ ~ ~
+DEAL::~ 0.3750000 ~ -0.1250000 ~ 0.7500000 ~ ~ ~
+DEAL::~ -0.1250000 ~ 0.3750000 ~ 0.7500000 ~ ~ ~
+DEAL::0.3750000 -0.1250000 ~ ~ ~ ~ 0.7500000 ~ ~
+DEAL::-0.1250000 0.3750000 ~ ~ ~ ~ 0.7500000 ~ ~
+DEAL::~ ~ 0.3750000 -0.1250000 ~ ~ ~ 0.7500000 ~
+DEAL::~ ~ -0.1250000 0.3750000 ~ ~ ~ 0.7500000 ~
+DEAL::0.1406250 -0.04687500 -0.04687500 0.01562500 0.2812500 -0.09375000 0.2812500 -0.09375000 0.5625000
+DEAL::-0.04687500 0.1406250 0.01562500 -0.04687500 -0.09375000 0.2812500 0.2812500 -0.09375000 0.5625000
+DEAL::-0.04687500 0.01562500 0.1406250 -0.04687500 0.2812500 -0.09375000 -0.09375000 0.2812500 0.5625000
+DEAL::0.01562500 -0.04687500 -0.04687500 0.1406250 -0.09375000 0.2812500 -0.09375000 0.2812500 0.5625000
+DEAL::Q2<3> restriction 0
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 1
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 2
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 3
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 4
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 5
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 6
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::Q2<3> restriction 7
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::~ ~ ~ ~ 1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
+DEAL::1.00000 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_q_hierarchical.h>
#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_bdm.h>
--- /dev/null
+//----------------------------------------------------------------------
+// $Id: fe_support_points_q_dg0.cc 26376 2012-09-14 12:32:17Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2007, 2008, 2010, 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+//
+// Compute support points
+
+#include "../tests.h"
+#include "fe_support_points_common.h"
+
+
+int
+main()
+{
+ initlog(__FILE__);
+ deallog.threshold_double(1.e-10);
+
+ CHECK_ALL(Q_DG0,1,2);
+ CHECK_ALL(Q_DG0,2,2);
+ CHECK_ALL(Q_DG0,3,2);
+
+ CHECK_ALL(Q_DG0,1,3);
+ CHECK_ALL(Q_DG0,2,3);
+}
--- /dev/null
+
+DEAL::Q1<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::Q1<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::Q1<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::Q1<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::Q1<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::Q2<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.00 0.500
+DEAL::1.00 0.500
+DEAL::0.500 0.00
+DEAL::0.500 1.00
+DEAL::0.500 0.500
+DEAL::Q2<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::0.00 0.500
+DEAL::Q2<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::1.00 0.500
+DEAL::Q2<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.500 0.00
+DEAL::Q2<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.500 1.00
+DEAL::Q3<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.00 0.333
+DEAL::0.00 0.667
+DEAL::1.00 0.333
+DEAL::1.00 0.667
+DEAL::0.333 0.00
+DEAL::0.667 0.00
+DEAL::0.333 1.00
+DEAL::0.667 1.00
+DEAL::0.333 0.333
+DEAL::0.667 0.333
+DEAL::0.333 0.667
+DEAL::0.667 0.667
+DEAL::Q3<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::0.00 0.333
+DEAL::0.00 0.667
+DEAL::Q3<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::1.00 0.333
+DEAL::1.00 0.667
+DEAL::Q3<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.333 0.00
+DEAL::0.667 0.00
+DEAL::Q3<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.333 1.00
+DEAL::0.667 1.00
+DEAL::Q1<3> cell support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::Q1<3> face 0 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::Q1<3> face 1 support points
+DEAL::1.00 0.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::Q1<3> face 2 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::Q1<3> face 3 support points
+DEAL::0.00 1.00 0.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 1.00 1.00
+DEAL::Q1<3> face 4 support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::Q1<3> face 5 support points
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::Q2<3> cell support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::0.00 0.500 0.00
+DEAL::1.00 0.500 0.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 1.00 0.00
+DEAL::0.00 0.500 1.00
+DEAL::1.00 0.500 1.00
+DEAL::0.500 0.00 1.00
+DEAL::0.500 1.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::1.00 0.00 0.500
+DEAL::0.00 1.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::0.00 0.500 0.500
+DEAL::1.00 0.500 0.500
+DEAL::0.500 0.00 0.500
+DEAL::0.500 1.00 0.500
+DEAL::0.500 0.500 0.00
+DEAL::0.500 0.500 1.00
+DEAL::0.500 0.500 0.500
+DEAL::Q2<3> face 0 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::0.00 1.00 0.500
+DEAL::0.00 0.500 0.00
+DEAL::0.00 0.500 1.00
+DEAL::0.00 0.500 0.500
+DEAL::Q2<3> face 1 support points
+DEAL::1.00 0.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::1.00 0.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::1.00 0.500 0.00
+DEAL::1.00 0.500 1.00
+DEAL::1.00 0.500 0.500
+DEAL::Q2<3> face 2 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 0.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::1.00 0.00 0.500
+DEAL::0.500 0.00 0.500
+DEAL::Q2<3> face 3 support points
+DEAL::0.00 1.00 0.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 1.00 1.00
+DEAL::0.500 1.00 0.00
+DEAL::0.500 1.00 1.00
+DEAL::0.00 1.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::0.500 1.00 0.500
+DEAL::Q2<3> face 4 support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.500 0.00
+DEAL::1.00 0.500 0.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 1.00 0.00
+DEAL::0.500 0.500 0.00
+DEAL::Q2<3> face 5 support points
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::0.00 0.500 1.00
+DEAL::1.00 0.500 1.00
+DEAL::0.500 0.00 1.00
+DEAL::0.500 1.00 1.00
+DEAL::0.500 0.500 1.00
--- /dev/null
+//----------------------------------------------------------------------
+// $Id: fe_support_points_q.cc 26376 2012-09-14 12:32:17Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2007, 2008, 2010, 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+//
+// Compute support points
+
+#include "../tests.h"
+#include "fe_support_points_common.h"
+
+
+
+int
+main()
+{
+ initlog(__FILE__);
+ deallog.threshold_double(1.e-10);
+
+ CHECK_ALL(Q_DG0,1,2);
+ CHECK_ALL(Q_DG0,2,2);
+ CHECK_ALL(Q_DG0,3,2);
+
+ CHECK_ALL(Q_DG0,1,3);
+ CHECK_ALL(Q_DG0,2,3);
+}
--- /dev/null
+
+DEAL::Q_DG01<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.500 0.500
+DEAL::Q_DG01<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::Q_DG01<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::Q_DG01<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::Q_DG01<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::Q_DG02<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.00 0.500
+DEAL::1.00 0.500
+DEAL::0.500 0.00
+DEAL::0.500 1.00
+DEAL::0.500 0.500
+DEAL::0.500 0.500
+DEAL::Q_DG02<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::0.00 0.500
+DEAL::Q_DG02<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::1.00 0.500
+DEAL::Q_DG02<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.500 0.00
+DEAL::Q_DG02<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.500 1.00
+DEAL::Q_DG03<2> cell support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.00 0.333
+DEAL::0.00 0.667
+DEAL::1.00 0.333
+DEAL::1.00 0.667
+DEAL::0.333 0.00
+DEAL::0.667 0.00
+DEAL::0.333 1.00
+DEAL::0.667 1.00
+DEAL::0.333 0.333
+DEAL::0.667 0.333
+DEAL::0.333 0.667
+DEAL::0.667 0.667
+DEAL::0.500 0.500
+DEAL::Q_DG03<2> face 0 support points
+DEAL::0.00 0.00
+DEAL::0.00 1.00
+DEAL::0.00 0.333
+DEAL::0.00 0.667
+DEAL::Q_DG03<2> face 1 support points
+DEAL::1.00 0.00
+DEAL::1.00 1.00
+DEAL::1.00 0.333
+DEAL::1.00 0.667
+DEAL::Q_DG03<2> face 2 support points
+DEAL::0.00 0.00
+DEAL::1.00 0.00
+DEAL::0.333 0.00
+DEAL::0.667 0.00
+DEAL::Q_DG03<2> face 3 support points
+DEAL::0.00 1.00
+DEAL::1.00 1.00
+DEAL::0.333 1.00
+DEAL::0.667 1.00
+DEAL::Q_DG01<3> cell support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::0.500 0.500 0.500
+DEAL::Q_DG01<3> face 0 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::Q_DG01<3> face 1 support points
+DEAL::1.00 0.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::Q_DG01<3> face 2 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::Q_DG01<3> face 3 support points
+DEAL::0.00 1.00 0.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 1.00 1.00
+DEAL::Q_DG01<3> face 4 support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::Q_DG01<3> face 5 support points
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::Q_DG02<3> cell support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::0.00 0.500 0.00
+DEAL::1.00 0.500 0.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 1.00 0.00
+DEAL::0.00 0.500 1.00
+DEAL::1.00 0.500 1.00
+DEAL::0.500 0.00 1.00
+DEAL::0.500 1.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::1.00 0.00 0.500
+DEAL::0.00 1.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::0.00 0.500 0.500
+DEAL::1.00 0.500 0.500
+DEAL::0.500 0.00 0.500
+DEAL::0.500 1.00 0.500
+DEAL::0.500 0.500 0.00
+DEAL::0.500 0.500 1.00
+DEAL::0.500 0.500 0.500
+DEAL::0.500 0.500 0.500
+DEAL::Q_DG02<3> face 0 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::0.00 1.00 0.500
+DEAL::0.00 0.500 0.00
+DEAL::0.00 0.500 1.00
+DEAL::0.00 0.500 0.500
+DEAL::Q_DG02<3> face 1 support points
+DEAL::1.00 0.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::1.00 0.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::1.00 0.500 0.00
+DEAL::1.00 0.500 1.00
+DEAL::1.00 0.500 0.500
+DEAL::Q_DG02<3> face 2 support points
+DEAL::0.00 0.00 0.00
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 0.00
+DEAL::1.00 0.00 1.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 0.00 1.00
+DEAL::0.00 0.00 0.500
+DEAL::1.00 0.00 0.500
+DEAL::0.500 0.00 0.500
+DEAL::Q_DG02<3> face 3 support points
+DEAL::0.00 1.00 0.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 0.00
+DEAL::1.00 1.00 1.00
+DEAL::0.500 1.00 0.00
+DEAL::0.500 1.00 1.00
+DEAL::0.00 1.00 0.500
+DEAL::1.00 1.00 0.500
+DEAL::0.500 1.00 0.500
+DEAL::Q_DG02<3> face 4 support points
+DEAL::0.00 0.00 0.00
+DEAL::1.00 0.00 0.00
+DEAL::0.00 1.00 0.00
+DEAL::1.00 1.00 0.00
+DEAL::0.00 0.500 0.00
+DEAL::1.00 0.500 0.00
+DEAL::0.500 0.00 0.00
+DEAL::0.500 1.00 0.00
+DEAL::0.500 0.500 0.00
+DEAL::Q_DG02<3> face 5 support points
+DEAL::0.00 0.00 1.00
+DEAL::1.00 0.00 1.00
+DEAL::0.00 1.00 1.00
+DEAL::1.00 1.00 1.00
+DEAL::0.00 0.500 1.00
+DEAL::1.00 0.500 1.00
+DEAL::0.500 0.00 1.00
+DEAL::0.500 1.00 1.00
+DEAL::0.500 0.500 1.00
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_q_hierarchical.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_system.h>
--- /dev/null
+//---------------------------- injection_dgq.cc ---------------------------
+// $Id: injection_q_dg0.cc 15183 2007-09-10 01:10:03Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2006, 2007 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- injection_dgq.cc ---------------------------
+
+
+char logname[] = "injection_dgq/output";
+
+
+#include "injection_common.h"
+
+template <int dim>
+void test ()
+{
+ for (unsigned int i=1; i<4; ++i)
+ for (unsigned int j=i; j<4; ++j)
+ do_check (FE_Q_DG0<dim>(i), FE_Q_DG0<dim>(j));
+}
#include <deal.II/base/function.h>
#include <deal.II/lac/vector.h>
#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <iomanip>
--- /dev/null
+//----------------------------------------------------------------------
+// $Id: interpolate_q_dg0.cc 23710 2011-05-17 04:50:10Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2005 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+
+#include "interpolate_common.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <fstream>
+
+// FE_Q_DG0<dim>::interpolate(...)
+
+template <int dim>
+void check_q_dg0(const Function<dim>& f,
+ const unsigned int degree)
+{
+ FE_Q_DG0<dim> fe(degree);
+ deallog << fe.get_name() << ' ';
+
+ std::vector<double> dofs(fe.dofs_per_cell);
+
+ std::vector<std::vector<double> >
+ values(1, std::vector<double>(fe.get_unit_support_points().size()));
+ f.value_list(fe.get_unit_support_points(), values[0]);
+ fe.interpolate(dofs, values[0]);
+ deallog << " value " << difference(fe,dofs,f);
+ fe.interpolate(dofs, values);
+ deallog << " vector " << difference(fe,dofs,f);
+
+ std::vector<Vector<double> >
+ vectors(fe.get_unit_support_points().size(), Vector<double>(1));
+ f.vector_value_list(fe.get_unit_support_points(), vectors);
+ fe.interpolate(dofs, vectors, 0);
+ deallog << " Vector " << difference(fe,dofs,f) << std::endl;
+}
+
+int main()
+{
+ std::ofstream logfile ("interpolate_q_dg0/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(2.e-15);
+
+ Q1WedgeFunction<1,1> w1;
+ check_q_dg0(w1,1);
+ check_q_dg0(w1,2);
+ check_q_dg0(w1,3);
+ Q1WedgeFunction<2,1> w2;
+ check_q_dg0(w2,1);
+ check_q_dg0(w2,2);
+ check_q_dg0(w2,3);
+ Q1WedgeFunction<2,2> w22;
+ check_q_dg0(w22,2);
+ check_q_dg0(w22,3);
+ Q1WedgeFunction<2,3> w23;
+ check_q_dg0(w23,3);
+ Q1WedgeFunction<3,1> w3;
+ check_q_dg0(w3,1);
+ check_q_dg0(w3,2);
+ check_q_dg0(w3,3);
+}
--- /dev/null
+
+DEAL::FE_Q_DG0<1>(QUnknownNodes(1)) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<1>(QUnknownNodes(2)) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<1>(QUnknownNodes(3)) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(1) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(2) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(3) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(2) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(3) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<2>(3) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<3>(1) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<3>(2) value 0 vector 0 Vector 0
+DEAL::FE_Q_DG0<3>(3) value 0 vector 0 Vector 0
--- /dev/null
+//---------------------------- abfprojection_01.cc ---------------------------
+// abfprojection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
+// Version:
+//
+// Copyright (C) 2003, 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- abfprojection_01.cc ---------------------------
+
+/*
+ * Project the function [1,1] onto a deformed grid and see whether the
+ * FESystem elements can represent it exactly. This shouldn't be a surprise,
+ * but it is nice to compare with the RT and ABF elements
+ */
+
+
+
+char logname[] = "q_dg0_projection_01/output";
+#include "deformed_projection.h"
+
+
+void test ()
+{
+ FESystem<2> fe (FE_Q_DG0<2>(3), 2);
+ check (fe);
+}
--- /dev/null
+
+DEAL::Dofs/cell 32Dofs/face 8
+DEAL::Dofs total 32
+DEAL::MM created
+DEAL::RHS created
+DEAL:cg::Starting value 0.16
+DEAL:cg::Convergence step 14 value 0
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <solution_0> <solution_1>
+0.00 0.00 1.00 1.00
+0.25 0.00 1.00 1.00
+0.50 0.00 1.00 1.00
+0.75 0.00 1.00 1.00
+1.00 0.00 1.00 1.00
+
+0.15 0.12 1.00 1.00
+0.36 0.12 1.00 1.00
+0.57 0.12 1.00 1.00
+0.79 0.12 1.00 1.00
+1.00 0.12 1.00 1.00
+
+0.30 0.25 1.00 1.00
+0.47 0.25 1.00 1.00
+0.65 0.25 1.00 1.00
+0.82 0.25 1.00 1.00
+1.00 0.25 1.00 1.00
+
+0.45 0.38 1.00 1.00
+0.59 0.38 1.00 1.00
+0.73 0.38 1.00 1.00
+0.86 0.38 1.00 1.00
+1.00 0.38 1.00 1.00
+
+0.60 0.50 1.00 1.00
+0.70 0.50 1.00 1.00
+0.80 0.50 1.00 1.00
+0.90 0.50 1.00 1.00
+1.00 0.50 1.00 1.00
+
+
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::Dofs total 32
+DEAL::MM created
+DEAL::RHS created
+DEAL:cg::Starting value 0.10
+DEAL:cg::Convergence step 14 value 0
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <solution_0> <solution_1>
+0.60 0.50 1.00 1.00
+0.70 0.50 1.00 1.00
+0.80 0.50 1.00 1.00
+0.90 0.50 1.00 1.00
+1.00 0.50 1.00 1.00
+
+0.57 0.62 1.00 1.00
+0.68 0.62 1.00 1.00
+0.79 0.62 1.00 1.00
+0.89 0.62 1.00 1.00
+1.00 0.62 1.00 1.00
+
+0.55 0.75 1.00 1.00
+0.66 0.75 1.00 1.00
+0.78 0.75 1.00 1.00
+0.89 0.75 1.00 1.00
+1.00 0.75 1.00 1.00
+
+0.53 0.88 1.00 1.00
+0.64 0.88 1.00 1.00
+0.76 0.88 1.00 1.00
+0.88 0.88 1.00 1.00
+1.00 0.88 1.00 1.00
+
+0.50 1.00 1.00 1.00
+0.62 1.00 1.00 1.00
+0.75 1.00 1.00 1.00
+0.88 1.00 1.00 1.00
+1.00 1.00 1.00 1.00
+
+
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::Dofs total 32
+DEAL::MM created
+DEAL::RHS created
+DEAL:cg::Starting value 0.19
+DEAL:cg::Convergence step 14 value 0
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <solution_0> <solution_1>
+0.00 0.00 1.00 1.00
+0.15 0.12 1.00 1.00
+0.30 0.25 1.00 1.00
+0.45 0.38 1.00 1.00
+0.60 0.50 1.00 1.00
+
+0.00 0.25 1.00 1.00
+0.14 0.34 1.00 1.00
+0.29 0.44 1.00 1.00
+0.43 0.53 1.00 1.00
+0.57 0.62 1.00 1.00
+
+0.00 0.50 1.00 1.00
+0.14 0.56 1.00 1.00
+0.28 0.62 1.00 1.00
+0.41 0.69 1.00 1.00
+0.55 0.75 1.00 1.00
+
+0.00 0.75 1.00 1.00
+0.13 0.78 1.00 1.00
+0.26 0.81 1.00 1.00
+0.39 0.84 1.00 1.00
+0.53 0.88 1.00 1.00
+
+0.00 1.00 1.00 1.00
+0.12 1.00 1.00 1.00
+0.25 1.00 1.00 1.00
+0.38 1.00 1.00 1.00
+0.50 1.00 1.00 1.00
+
+
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
--- /dev/null
+#!/usr/bin/gnuplot -persist
+#
+#
+# G N U P L O T
+# Version 4.0 patchlevel 0
+# last modified Thu Apr 15 14:44:22 CEST 2004
+# System: Linux 2.6.11.4-21.11-smp
+#
+# Copyright (C) 1986 - 1993, 1998, 2004, 2006
+# Thomas Williams, Colin Kelley and many others
+#
+# This is gnuplot version 4.0. Please refer to the documentation
+# for command syntax changes. The old syntax will be accepted
+# throughout the 4.0 series, but all save files use the new syntax.
+#
+# Type `help` to access the on-line reference manual.
+# The gnuplot FAQ is available from
+# http://www.gnuplot.info/faq/
+#
+# Send comments and requests for help to
+# <gnuplot-info@lists.sourceforge.net>
+# Send bugs, suggestions and mods to
+# <gnuplot-bugs@lists.sourceforge.net>
+#
+set terminal postscript eps color solid enhanced 20
+# set output
+unset clip points
+set clip one
+unset clip two
+set bar 1.000000
+set border 31 lt -1 lw 1.000
+set xdata
+set ydata
+set zdata
+set x2data
+set y2data
+set timefmt x "%d/%m/%y,%H:%M"
+set timefmt y "%d/%m/%y,%H:%M"
+set timefmt z "%d/%m/%y,%H:%M"
+set timefmt x2 "%d/%m/%y,%H:%M"
+set timefmt y2 "%d/%m/%y,%H:%M"
+set timefmt cb "%d/%m/%y,%H:%M"
+set boxwidth
+set style fill empty border
+set dummy x,y
+set format x "% g"
+set format y "% g"
+set format x2 "% g"
+set format y2 "% g"
+set format z "% g"
+set format cb "% g"
+set angles radians
+unset grid
+set key title ""
+set key right top Right noreverse enhanced box linetype -2 linewidth 1.000 samplen 4 spacing 1 width 0 height 0 autotitles
+unset label
+unset arrow
+unset style line
+unset style arrow
+unset logscale
+set offsets 0, 0, 0, 0
+set pointsize 1
+set encoding default
+unset polar
+unset parametric
+unset decimalsign
+set view 60, 20, 1, 1
+set samples 100, 100
+set isosamples 10, 10
+set surface
+unset contour
+set clabel '%8.3g'
+set mapping cartesian
+set datafile separator whitespace
+set hidden3d offset 1 trianglepattern 3 undefined 1 altdiagonal bentover
+set cntrparam order 4
+set cntrparam linear
+set cntrparam levels auto 5
+set cntrparam points 5
+set size ratio 0 1,1
+set origin 0,0
+set style data lines
+set style function lines
+set xzeroaxis lt -2 lw 1.000
+set yzeroaxis lt -2 lw 1.000
+set x2zeroaxis lt -2 lw 1.000
+set y2zeroaxis lt -2 lw 1.000
+set tics in
+set ticslevel 0.5
+set ticscale 1 0.5
+set mxtics default
+set mytics default
+set mztics default
+set mx2tics default
+set my2tics default
+set mcbtics default
+set xtics border mirror norotate autofreq
+set ytics border mirror norotate autofreq
+set ztics border nomirror norotate autofreq
+set nox2tics
+set noy2tics
+set cbtics border mirror norotate autofreq
+set title "" 0.000000,0.000000 font ""
+set timestamp "" bottom norotate 0.000000,0.000000 ""
+set rrange [ * : * ] noreverse nowriteback # (currently [0.00000:10.0000] )
+set trange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] )
+set urange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] )
+set vrange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] )
+set xlabel "x" 0.000000,0.000000 font ""
+set x2label "" 0.000000,0.000000 font ""
+set xrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set x2range [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set ylabel "y" 0.000000,0.000000 font ""
+set y2label "" 0.000000,0.000000 font ""
+set yrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set y2range [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set zlabel "" 0.000000,0.000000 font ""
+set zrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set cblabel "" 0.000000,0.000000 font ""
+set cbrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] )
+set zero 1e-08
+set lmargin -1
+set bmargin -1
+set rmargin -1
+set tmargin -1
+set locale "C"
+set pm3d scansautomatic flush begin noftriangles nohidden3d implicit corners2color mean
+unset pm3d
+set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB
+set palette rgbformulae 7, 5, 15
+set colorbox default
+set colorbox vertical origin 0.9,0.2 size 0.1,0.63 bdefault
+set loadpath
+set fontpath
+set fit noerrorvariables
+set output "RT1_00_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:3 title "v0_x"
+set output "RT1_00_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:4 title "v0_y"
+set output "RT1_01_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:5 title "v1_x"
+set output "RT1_01_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:6 title "v1_y"
+set output "RT1_02_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:7 title "v2_x"
+set output "RT1_02_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:8 title "v2_y"
+set output "RT1_03_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:9 title "v3_x"
+set output "RT1_03_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:10 title "v3_y"
+set output "RT1_04_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:11 title "v4_x"
+set output "RT1_04_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:12 title "v4_y"
+set output "RT1_05_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:13 title "v5_x"
+set output "RT1_05_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:14 title "v5_y"
+set output "RT1_06_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:15 title "v6_x"
+set output "RT1_06_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:16 title "v6_y"
+set output "RT1_07_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:17 title "v7_x"
+set output "RT1_07_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:18 title "v7_y"
+set output "RT1_08_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:19 title "v8_x"
+set output "RT1_08_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:20 title "v8_y"
+set output "RT1_09_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:21 title "v9_x"
+set output "RT1_09_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:22 title "v9_y"
+set output "RT1_10_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:23 title "v10_x"
+set output "RT1_10_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:24 title "v10_y"
+set output "RT1_11_x.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:25 title "v11_x"
+set output "RT1_11_y.eps"
+splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:26 title "v11_y"
+# EOF