From: Timo Heister Date: Tue, 8 Jan 2013 17:12:22 +0000 (+0000) Subject: reintegrate branch branch_feqdg0 X-Git-Tag: v8.0.0~1630 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4ad143ec219d2189ab48b6114e8bb1593553dc41;p=dealii.git reintegrate branch branch_feqdg0 git-svn-id: https://svn.dealii.org/trunk@27974 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/base/tensor_product_polynomials_const.h b/deal.II/include/deal.II/base/tensor_product_polynomials_const.h new file mode 100644 index 0000000000..570e1096ab --- /dev/null +++ b/deal.II/include/deal.II/base/tensor_product_polynomials_const.h @@ -0,0 +1,351 @@ +#ifndef __deal2__tensor_product_polynomials_const_h +#define __deal2__tensor_product_polynomials_const_h + + +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * @addtogroup Polynomials + * @{ + */ + +/** + * Tensor product of given polynomials. + * + * Given a vector of n one-dimensional polynomials + * P1 to Pn, this class generates + * ndim polynomials of the form + * Qijk(x,y,z) = + * Pi(x)Pj(y)Pk(z)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]dim or [0,1] + * dim, 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 P1(x)P1(y), + * P2(x)P1(y), P3(x)P1(y), + * ..., P1(x)P2(y), + * P2(x)P2(y), P3(x)P2(y), + * ... 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 +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. pols 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 + TensorProductPolynomialsConst (const std::vector &pols); + + /** + * Prints the list of the indices + * to out. + */ + void output_indices(std::ostream &out) const; + + /** + * Sets the ordering of the + * polynomials. Requires + * renumber.size()==n(). + * Stores a copy of + * renumber. + */ + void set_numbering(const std::vector &renumber); + + /** + * Gives read access to the + * renumber vector. + */ + const std::vector &get_numbering() const; + + /** + * Gives read access to the + * inverse renumber vector. + */ + const std::vector &get_numbering_inverse() const; + + /** + * Computes the value and the + * first and second derivatives + * of each tensor product + * polynomial at unit_point. + * + * 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 &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const; + + /** + * Computes the value of the + * ith tensor product + * polynomial at + * unit_point. Here i 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 + * values.size()==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 &p) const; + + /** + * Computes the grad of the + * ith tensor product + * polynomial at + * unit_point. Here i 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 + * grads.size()==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 &p) const; + + /** + * Computes the second + * derivative (grad_grad) of the + * ith tensor product + * polynomial at + * unit_point. Here i 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 + * grad_grads.size()==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 &p) const; + + /** + * Returns the number of tensor + * product polynomials plus the constant + * function. For n 1d polynomials + * this is ndim+1. + */ + unsigned int n () const; + + +private: + /** + * Copy of the vector pols of + * polynomials given to the + * constructor. + */ + std::vector > polynomials; + + /** + * Number of tensor product + * polynomials. See n(). + */ + unsigned int n_tensor_pols; + + /** + * Index map for reordering the + * polynomials. + */ + std::vector index_map; + + /** + * Index map for reordering the + * polynomials. + */ + std::vector index_map_inverse; + + /** + * Each tensor product polynomial + * 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. + */ + // 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 + * xdim for + * unsigned int x. Used in + * the constructor. + */ + static + unsigned int x_to_the_dim (const unsigned int x); +}; + +#ifndef DOXYGEN + +template +inline +const std::vector & +TensorProductPolynomialsConst::get_numbering() const +{ + return index_map; +} + + +template +inline +const std::vector & +TensorProductPolynomialsConst::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 +inline +unsigned int +TensorProductPolynomialsConst:: +x_to_the_dim (const unsigned int x) +{ + unsigned int y = 1; + for (int d=0; d +template +TensorProductPolynomialsConst:: +TensorProductPolynomialsConst(const std::vector &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 +#include +#include + +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 spacedim template parameter check + * the documentation of FiniteElement or the one of Triangulation. + * + *

Implementation

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

Numbering of the degrees of freedom (DoFs)

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

Q1 elements

+ *
    + *
  • 1D case: + * @verbatim + * 0---2---1 + * @endverbatim + * + *
  • 2D case: + * @verbatim + * 2-------3 + * | | + * | 5 | + * | | + * 0-------1 + * @endverbatim + * + *
  • 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: + *
      + *
    • Index 0: [ 0, 0, 0]; + *
    • Index 1: [ 1, 0, 0]; + *
    • Index 2: [ 0, 1, 0]; + *
    • Index 3: [ 1, 1, 0]; + *
    • Index 4: [ 0, 0, 1]; + *
    • Index 5: [ 1, 0, 1]; + *
    • Index 6: [ 0, 1, 1]; + *
    • Index 7: [ 1, 1, 1]; + *
    • Index 8: [1/2, 1/2, 1/2]; + *
    + *
+ *

Q2 elements

+ *
    + *
  • 1D case: + * @verbatim + * 0---2---1 + * @endverbatim + * Index 3 has the same coordinates as index 2 + * + *
  • 2D case: + * @verbatim + * 2---7---3 + * | | + * 4 8 5 + * | | + * 0---6---1 + * @endverbatim + * Index 9 has the same coordinates as index 2 + * + *
  • 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: + *
      + *
    • Index 0: [0, 0, 0]; + *
    • Index 1: [1, 0, 0]; + *
    • Index 2: [0, 1, 0]; + *
    • Index 3: [1, 1, 0]; + *
    • Index 4: [0, 0, 1]; + *
    • Index 5: [1, 0, 1]; + *
    • Index 6: [0, 1, 1]; + *
    • Index 7: [1, 1, 1]; + *
    • Index 8: [0, 1/2, 0]; + *
    • Index 9: [1, 1/2, 0]; + *
    • Index 10: [1/2, 0, 0]; + *
    • Index 11: [1/2, 1, 0]; + *
    • Index 12: [0, 1/2, 1]; + *
    • Index 13: [1, 1/2, 1]; + *
    • Index 14: [1/2, 0, 1]; + *
    • Index 15: [1/2, 1, 1]; + *
    • Index 16: [0, 0, 1/2]; + *
    • Index 17: [1, 0, 1/2]; + *
    • Index 18: [0, 1, 1/2]; + *
    • Index 19: [1, 1, 1/2]; + *
    • Index 20: [0, 1/2, 1/2]; + *
    • Index 21: [1, 1/2, 1/2]; + *
    • Index 22: [1/2, 0, 1/2]; + *
    • Index 23: [1/2, 1, 1/2]; + *
    • Index 24: [1/2, 1/2, 0]; + *
    • Index 25: [1/2, 1/2, 1]; + *
    • Index 26: [1/2, 1/2, 1/2]; + *
    • Index 27: [1/2, 1/2, 1/2]; + *
    + *
+ *

Q3 elements

+ *
    + *
  • 1D case: + * @verbatim + * 0--2-4-3--1 + * @endverbatim + * + *
  • 2D case: + * @verbatim + * 2--10-11-3 + * | | + * 5 14 15 7 + * | 16 | + * 4 12 13 6 + * | | + * 0--8--9--1 + * @endverbatim + *
+ *

Q4 elements

+ *
    + *
  • 1D case: + * @verbatim + * 0--2--3--4--1 + * @endverbatim + * Index 5 has the same coordinates as index 3 + * + *
  • 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 + *
+ * + */ +template +class FE_Q_DG0 : public FE_Poly,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 + * points.size()-1. + * 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 + * FE_Q_DG0(degree), 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 &local_dofs, + const std::vector &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, offset 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 &local_dofs, + const std::vector > &values, + unsigned int offset = 0) const; + + /** + * Interpolate a set of vector + * values, computed in the + * generalized support points. + */ + virtual void interpolate( + std::vector &local_dofs, + const VectorSlice > > &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 + * source.dofs_per_cell. + * + * These matrices are only available if the + * source element is also a @p FE_Q_DG0 + * element. Otherwise, an exception of type + * FiniteElement::ExcInterpolationNotImplemented + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + + + /** + * Return the matrix + * interpolating from a face of + * of one element to the face of + * the neighboring element. + * The size of the matrix is + * then source.dofs_per_face times + * this->dofs_per_face. + * + * 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::ExcInterpolationNotImplemented. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + + /** + * Return the matrix + * interpolating from a face of + * of one element to the face of + * the neighboring element. + * The size of the matrix is + * then source.dofs_per_face times + * this->dofs_per_face. + * + * 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::ExcInterpolationNotImplemented. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const; + + /** + * Check for non-zero values on a face. + * + * This function returns + * @p true, if the shape + * function @p shape_index has + * non-zero values on the face + * @p face_index. + * + * Implementation of the + * interface in + * FiniteElement + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + /** + * @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 > + hp_vertex_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as + * hp_vertex_dof_indices(), + * except that the function + * treats degrees of freedom on + * lines. + */ + virtual + std::vector > + hp_line_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as + * hp_vertex_dof_indices(), + * except that the function + * treats degrees of freedom on + * quads. + */ + virtual + std::vector > + hp_quad_dof_identities (const FiniteElement &fe_other) const; + + /** + * Return whether this element dominates + * the one given as argument when they + * meet at a common face, + * whether it is the other way around, + * whether neither dominates, or if + * either could dominate. + * + * For a definition of domination, see + * FiniteElementBase::Domination and in + * particular the @ref hp_paper "hp paper". + */ + virtual + FiniteElementDomination::Domination + compare_for_face_domination (const FiniteElement &fe_other) const; + //@} + + /** + * 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 *clone() const; + +private: + /** + * Returns the restriction_is_additive flags. + * Only the last component is true. + */ + static std::vector 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 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 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 + 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 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 friend class FE_Q_DG0; + + friend struct FE_Q_DG0::Implementation; +}; + + + +/*@}*/ + +/* -------------- declaration of explicit specializations ------------- */ + +template <> +void FE_Q_DG0<1>::initialize_unit_face_support_points (); + +template <> +std::vector +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 +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 +FE_Q_DG0<1,3>::face_lexicographic_to_hierarchic_numbering (const unsigned int); + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/source/base/tensor_product_polynomials_const.cc b/deal.II/source/base/tensor_product_polynomials_const.cc new file mode 100644 index 0000000000..27ca3417a3 --- /dev/null +++ b/deal.II/source/base/tensor_product_polynomials_const.cc @@ -0,0 +1,325 @@ +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +/* ------------------- TensorProductPolynomialsConst -------------- */ + +template <> +inline +void +TensorProductPolynomialsConst<1>:: +compute_index (const unsigned int i, + unsigned int (&indices)[1]) const +{ + Assert (i +inline +void +TensorProductPolynomialsConst<2>:: +compute_index (const unsigned int i, + unsigned int (&indices)[2]) const +{ + const unsigned int n_pols = polynomials.size(); + Assert (i +inline +void +TensorProductPolynomialsConst<3>:: +compute_index (const unsigned int i, + unsigned int (&indices)[3]) const +{ + const unsigned int n_pols = polynomials.size(); + Assert (i +void +TensorProductPolynomialsConst::output_indices(std::ostream &out) const +{ + unsigned int ix[dim]; + //ouput without constant locally constant function + for (unsigned int i=0; i +void +TensorProductPolynomialsConst::set_numbering( + const std::vector &renumber) +{ + Assert(renumber.size()==index_map.size(), + ExcDimensionMismatch(renumber.size(), index_map.size())); + + index_map=renumber; + for (unsigned int i=0; i +double +TensorProductPolynomialsConst::compute_value (const unsigned int i, + const Point &p) const +{ + unsigned int max_indices=pow(polynomials.size(),dim); + + Assert (i<=max_indices, ExcInternalError()); + if (i +double +TensorProductPolynomialsConst<0>::compute_value (const unsigned int , + const Point<0> &) const +{ + Assert (false, ExcNotImplemented()); + return 0.; +} + +template +Tensor<1,dim> +TensorProductPolynomialsConst::compute_grad (const unsigned int i, + const Point &p) const +{ + unsigned int max_indices=pow(polynomials.size(),dim); + + Assert (i<=max_indices, ExcInternalError()); + Tensor<1,dim> grad; + if (i tmp (2); + for (unsigned int d=0; d +Tensor<2,dim> +TensorProductPolynomialsConst::compute_grad_grad (const unsigned int i, + const Point &p) const +{ + unsigned int max_indices=pow(polynomials.size(),dim); + + Assert (i<=max_indices, ExcInternalError()); + Tensor<2,dim> grad_grad; + if (i tmp (3); + for (unsigned int d=0; d +void +TensorProductPolynomialsConst:: +compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &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 tmp (n_values_and_derivatives); + for (unsigned int d=0; d +unsigned int +TensorProductPolynomialsConst::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 diff --git a/deal.II/source/fe/fe_poly.cc b/deal.II/source/fe/fe_poly.cc index ac5052f76e..c2be11407a 100644 --- a/deal.II/source/fe/fe_poly.cc +++ b/deal.II/source/fe/fe_poly.cc @@ -13,6 +13,7 @@ #include #include +#include #include #include #include diff --git a/deal.II/source/fe/fe_poly.inst.in b/deal.II/source/fe/fe_poly.inst.in index 605b020ae4..061f74b81a 100644 --- a/deal.II/source/fe/fe_poly.inst.in +++ b/deal.II/source/fe/fe_poly.inst.in @@ -16,6 +16,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS { #if deal_II_dimension <= deal_II_space_dimension template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; + template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; #endif diff --git a/deal.II/source/fe/fe_q_dg0.cc b/deal.II/source/fe/fe_q_dg0.cc new file mode 100644 index 0000000000..ce21feb269 --- /dev/null +++ b/deal.II/source/fe/fe_q_dg0.cc @@ -0,0 +1,1811 @@ +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include + +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 struct int2type {}; +// +// namespace { +// static void SYMBOL (const int2type<1> & ) {} +// } +// +// template void g() { +// SYMBOL(int2type()); +// } +// +// 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 + invert_numbering (const std::vector &in) + { + std::vector out (in.size()); + for (unsigned int i=0; i +struct FE_Q_DG0::Implementation +{ + /** + * Initialize the hanging node + * constraints matrices. Called from the + * constructor in case the finite element + * is based on quadrature points. + */ + template + static + void initialize_constraints (const Quadrature<1> &, + FE_Q_DG0<1,spacedim> &) + { + // no constraints in 1d + } + + + template + 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 > constraint_points; + // Add midpoint + constraint_points.push_back (Point (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::child_to_cell_coordinates(Point(i*step),0)); + // subface 1 + for (unsigned int i=1; i<=n; ++i) + constraint_points.push_back ( + GeometryInfo::child_to_cell_coordinates(Point(i*step),1)); + } + + // Now construct relation between + // destination (child) and source (mother) + // dofs. + const std::vector > 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 + 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 > constraint_points; + + // Add midpoint + constraint_points.push_back (Point (0.5, 0.5)); + + // Add midpoints of lines of + // "mother-face" + constraint_points.push_back (Point (0, 0.5)); + constraint_points.push_back (Point (1, 0.5)); + constraint_points.push_back (Point (0.5, 0)); + constraint_points.push_back (Point (0.5, 1)); + + if (fe.degree>1) + { + const unsigned int n=fe.degree-1; + const double step=1./fe.degree; + std::vector > line_support_points(n); + for (unsigned int i=0; i qline(line_support_points); + + // auxiliary points in 2d + std::vector > p_line(n); + + // Add nodes of lines interior + // in the "mother-face" + + // line 5: use line 9 + QProjector::project_to_subface(qline, 0, 0, p_line); + for (unsigned int i=0; i (0.5, 0)); + // line 6: use line 10 + QProjector::project_to_subface(qline, 0, 1, p_line); + for (unsigned int i=0; i (0.5, 0)); + // line 7: use line 13 + QProjector::project_to_subface(qline, 2, 0, p_line); + for (unsigned int i=0; i (0, 0.5)); + // line 8: use line 14 + QProjector::project_to_subface(qline, 2, 1, p_line); + for (unsigned int i=0; i (0, 0.5)); + + // DoFs on bordering lines + // lines 9-16 + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int subface=0; + subface::max_children_per_face; ++subface) + { + QProjector::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 > 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 (ix*step, iy*step); + + // at the moment do this for + // isotropic face refinement only + for (unsigned int child=0; + child::max_children_per_cell; ++child) + for (unsigned int i=0; i::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 > polynomial_basis= + Polynomials::generate_complete_Lagrange_basis(points.get_points()); + + const TensorProductPolynomialsConst face_polynomials(polynomial_basis); + + fe.interface_constraints + .TableBase<2,double>::reinit (fe.interface_constraints_size()); + + for (unsigned int i=0; i 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 (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 +FE_Q_DG0::FE_Q_DG0 (const unsigned int degree) + : + FE_Poly, dim, spacedim> ( + TensorProductPolynomialsConst(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)), + FiniteElementData(get_dpo_vector(degree), + 1, degree, + FiniteElementData::L2), + get_riaf_vector(degree), + std::vector(1ll , std::vector(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 renumber (this->dofs_per_cell-1); + const FiniteElementData 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 tmp (degree); + this->prolongation = tmp.prolongation; + this->restriction = tmp.restriction; + }*/ + + initialize_quad_dof_index_permutation(); +} + + + +template +FE_Q_DG0::FE_Q_DG0 (const Quadrature<1> &points) + : + FE_Poly, dim, spacedim> ( + TensorProductPolynomialsConst(Polynomials::generate_complete_Lagrange_basis(points.get_points())), + FiniteElementData(get_dpo_vector(points.size()-1), + 1, points.size()-1, + FiniteElementData::L2), + get_riaf_vector(points.size()-1), + std::vector(1, std::vector(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 renumber (this->dofs_per_cell-1); + const FiniteElementData 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 +std::string +FE_Q_DG0::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 points(n_points); + const unsigned int dofs_per_cell = this->dofs_per_cell; + const std::vector > &unit_support_points = this->unit_support_points; + unsigned int index = 0; + + // Decode the support points + // in one coordinate direction. + for (unsigned int j=0; j1) ? (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; jdegree) > 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(QGaussLobatto(" << this->degree+1 << "))"; + else + namebuf << "FE_Q_DG0<" << dim << ">(QUnknownNodes(" << this->degree << "))"; + } + return namebuf.str(); +} + + + +template +FiniteElement * +FE_Q_DG0::clone() const +{ + return new FE_Q_DG0(*this); +} + +template +void +FE_Q_DG0::interpolate(std::vector &local_dofs, + const std::vector &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 +void +FE_Q_DG0::interpolate(std::vector &local_dofs, + const std::vector > &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; idofs_per_cell-1; ++i) + { + const std::pair 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 +void +FE_Q_DG0::interpolate( + std::vector &local_dofs, + const VectorSlice > > &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; idofs_per_cell-1; ++i) + { + const std::pair 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 +void +FE_Q_DG0:: +get_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q_DG0 element + typedef FE_Q_DG0 FEQDG0; + typedef FiniteElement FEL; + + AssertThrow ((x_source_fe.get_name().find ("FE_Q_DG0<") == 0) + || + (dynamic_cast(&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 &source_fe + = dynamic_cast&>(x_source_fe); + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix source_interpolation (this->dofs_per_cell, + source_fe.dofs_per_cell); + FullMatrix tmp (this->dofs_per_cell, + source_fe.dofs_per_cell); + for (unsigned int j=0; jdofs_per_cell-1; ++j) + { + // read in a point on this + // cell and evaluate the + // shape functions there + const Point p = this->unit_support_points[j]; + for (unsigned int i=0; idofs_per_cell-1; ++i) + cell_interpolation(j,i) = this->poly_space.compute_value (i, p); + + for (unsigned int i=0; idofs_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; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++i) + { + double sum = 0.; + for (unsigned int j=0; j +void +FE_Q_DG0<1>:: +get_face_interpolation_matrix (const FiniteElement<1,1> &/*x_source_fe*/, + FullMatrix &/*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 &/*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 &/*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 &/*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 &/*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 &/*interpolation_matrix*/) const +{ + typedef FiniteElement<1,3> FEL; + Assert (false, + FEL:: + ExcInterpolationNotImplemented ()); +} + + + +template +void +FE_Q_DG0:: +get_face_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q_DG0 element + typedef FE_Q_DG0 FEQDG0; + typedef FiniteElement FEL; + AssertThrow ((x_source_fe.get_name().find ("FE_Q_DG0<") == 0) + || + (dynamic_cast(&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 &source_fe + = dynamic_cast&>(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 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 face_quadrature + = QProjector::project_to_face (quad_face_support, 0); + for (unsigned int i=0; i &p = face_quadrature.point (i); + + for (unsigned int j=0; jdofs_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; jdofs_per_face; ++i) + sum += interpolation_matrix(j,i); + + Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1), + ExcInternalError()); + } +} + + + +template +void +FE_Q_DG0:: +get_subface_interpolation_matrix (const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix &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 *source_fe + = dynamic_cast *>(&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:: + ExcInterpolationNotImplemented ())); + + // generate a point on this + // cell and evaluate the + // shape functions there + const Quadrature + 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 subface_quadrature + = QProjector::project_to_subface (quad_face_support, 0, subface); + for (unsigned int i=0; idofs_per_face; ++i) + { + const Point &p = subface_quadrature.point (i); + + for (unsigned int j=0; jdofs_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; jdofs_per_face; ++j) + { + double sum = 0.; + + for (unsigned int i=0; idofs_per_face; ++i) + sum += interpolation_matrix(j,i); + + Assert (std::fabs(sum-1) < 2e-13*this->degree*dim, + ExcInternalError()); + } + } + else if (dynamic_cast *>(&x_source_fe) != 0) + { + // nothing to do here, the + // FE_Nothing has no degrees of + // freedom anyway + } + else + AssertThrow (false, + (typename FiniteElement:: + ExcInterpolationNotImplemented())); +} + + + +template +bool +FE_Q_DG0::hp_constraints_are_implemented () const +{ + return true; +} + + + + +template +std::vector > +FE_Q_DG0:: +hp_vertex_dof_identities (const FiniteElement &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*>(&fe_other) != 0) + { + return + std::vector > + (1, std::make_pair (0U, 0U)); + } + else if (dynamic_cast*>(&fe_other) != 0) + { + // the FE_Nothing has no + // degrees of freedom, so there + // are no equivalencies to be + // recorded + return std::vector > (); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector > (); + } +} + + + +template +std::vector > +FE_Q_DG0:: +hp_line_dof_identities (const FiniteElement &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 *fe_q_dg0_other = dynamic_cast*>(&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 > identities; + + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + const std::vector &index_map_inverse_other= + fe_q_dg0_other->poly_space.get_numbering_inverse(); + + for (unsigned int i=0; iunit_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*>(&fe_other) != 0) + { + // the FE_Nothing has no + // degrees of freedom, so there + // are no equivalencies to be + // recorded + return std::vector > (); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector > (); + } +} + + + +template +std::vector > +FE_Q_DG0:: +hp_quad_dof_identities (const FiniteElement &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 *fe_q_dg0_other = dynamic_cast*>(&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 > identities; + + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + const std::vector &index_map_inverse_other= + fe_q_dg0_other->poly_space.get_numbering_inverse(); + + for (unsigned int i1=0; i1unit_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*>(&fe_other) != 0) + { + // the FE_Nothing has no + // degrees of freedom, so there + // are no equivalencies to be + // recorded + return std::vector > (); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector > (); + } +} + + + +template +FiniteElementDomination::Domination +FE_Q_DG0:: +compare_for_face_domination (const FiniteElement &fe_other) const +{ + if (const FE_Q_DG0 *fe_q_dg0_other + = dynamic_cast*>(&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*>(&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 +void FE_Q_DG0::initialize_unit_support_points () +{ + // number of points: (degree+1)^dim+1 + unsigned int n = this->degree+1; + for (unsigned int i=1; idegree+1; + n++; + + this->unit_support_points.resize(n); + + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + + const double step = 1./this->degree; + Point 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 +void FE_Q_DG0::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; idegree+1; + n++; + + this->unit_support_points.resize(n); + + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + + Quadrature support_quadrature(points); + + Point p; + + for (unsigned int k=0; kunit_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 +void FE_Q_DG0::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; idegree+1; + + this->unit_face_support_points.resize(n); + + const std::vector &face_index_map_inverse= + FE_Q_DG0_Helper::invert_numbering(face_index_map); + + const double step = 1./this->degree; + Point 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 +void FE_Q_DG0::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; idegree+1; + + this->unit_face_support_points.resize(n); + + const std::vector< Point<1> > edge = points.get_points(); + + const std::vector &face_index_map_inverse= + FE_Q_DG0_Helper::invert_numbering(face_index_map); + + Point 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 +void +FE_Q_DG0::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; localdofs_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; idofs_per_line; ++i) + this->adjust_line_dof_index_for_line_orientation_table[i]=this->dofs_per_line-1-i - i; +} + +template +std::vector +FE_Q_DG0::get_riaf_vector(const unsigned int deg) +{ + std::vector riaf + (FiniteElementData (get_dpo_vector(deg),1,deg).dofs_per_cell, + false); + riaf[riaf.size()-1]=true; + return riaf; +} + + +template +std::vector +FE_Q_DG0::get_dpo_vector(const unsigned int deg) +{ + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +std::vector +FE_Q_DG0::get_dpo_vector_q(const unsigned int deg) +{ + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +std::vector +FE_Q_DG0::face_lexicographic_to_hierarchic_numbering (const unsigned int degree) +{ + const FiniteElementData face_data(FE_Q_DG0::get_dpo_vector_q(degree),1,degree); + std::vector face_renumber (face_data.dofs_per_cell); + FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber); + return face_renumber; +} + + + +template <> +std::vector +FE_Q_DG0<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +{ + return std::vector(); +} + +template <> +std::vector +FE_Q_DG0<1,2>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +{ + return std::vector(); +} + + +template <> +std::vector +FE_Q_DG0<1,3>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +{ + return std::vector(); +} + + +template +void +FE_Q_DG0::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 +bool +FE_Q_DG0::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::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::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::vertices_per_cell, + ExcInternalError()); + + for (unsigned int v=0; v::vertices_per_face; ++v) + if (GeometryInfo::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::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::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(quad_index) < + static_cast(GeometryInfo::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 +std::size_t +FE_Q_DG0::memory_consumption () const +{ + return (FiniteElement::memory_consumption() + + MemoryConsumption::memory_consumption (face_index_map)); +} + + + +// explicit instantiations +#include "fe_q_dg0.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_q_dg0.inst.in b/deal.II/source/fe/fe_q_dg0.inst.in new file mode 100644 index 0000000000..4f73f2c255 --- /dev/null +++ b/deal.II/source/fe/fe_q_dg0.inst.in @@ -0,0 +1,22 @@ +//--------------------------------------------------------------------------- +// $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; +#endif + } + + diff --git a/tests/fe/deformed_projection.h b/tests/fe/deformed_projection.h index fb0bfe0156..2179fa5ad3 100644 --- a/tests/fe/deformed_projection.h +++ b/tests/fe/deformed_projection.h @@ -52,6 +52,7 @@ #include #include #include +#include #include #include #include diff --git a/tests/fe/fe_prolongation_common.h b/tests/fe/fe_prolongation_common.h index e2d005ee96..3e21e582c9 100644 --- a/tests/fe/fe_prolongation_common.h +++ b/tests/fe/fe_prolongation_common.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include diff --git a/tests/fe/fe_prolongation_q_dg0.cc b/tests/fe/fe_prolongation_q_dg0.cc new file mode 100644 index 0000000000..e3dbcbf699 --- /dev/null +++ b/tests/fe/fe_prolongation_q_dg0.cc @@ -0,0 +1,34 @@ +//---------------------------------------------------------------------- +// $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); +} diff --git a/tests/fe/fe_prolongation_q_dg0/cmp/generic b/tests/fe/fe_prolongation_q_dg0/cmp/generic new file mode 100644 index 0000000000..5aa06d32dd --- /dev/null +++ b/tests/fe/fe_prolongation_q_dg0/cmp/generic @@ -0,0 +1,465 @@ + +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 diff --git a/tests/fe/fe_q_dg0.cc b/tests/fe/fe_q_dg0.cc new file mode 100644 index 0000000000..a0d26522cc --- /dev/null +++ b/tests/fe/fe_q_dg0.cc @@ -0,0 +1,969 @@ +// 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 +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace Step22 +{ + using namespace dealii; + + template + struct InnerPreconditioner; + + template <> + struct InnerPreconditioner<2> + { + typedef SparseDirectUMFPACK type; + }; + + template <> + struct InnerPreconditioner<3> + { + typedef SparseILU type; + }; + + template + class StokesProblem + { + public: + StokesProblem (const unsigned int degree, FESystem & 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 &calc_solution, + Vector &output_vector, + const Quadrature &quadrature, bool norm); + + const unsigned int degree; + + Triangulation triangulation; + FESystem &fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + + std_cxx1x::shared_ptr::type> + A_preconditioner; + + ConvergenceTable convergence_table; + }; + + template + class ExactSolution : public Function + { + public: + ExactSolution () : Function(dim+1) {} + + /*virtual*/ double value (const Point &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 &p, + const unsigned int component) const; + }; + + + template + double + ExactSolution::value (const Point &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 + Tensor<1,dim> + ExactSolution::gradient (const Point &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 + double + ExactSolution::laplacian (const Point &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 + class JumpFunction : public Function + { + public: + JumpFunction () : Function(1) {} + + double jump (const Point< dim> &p, + const Point &normal) const; + }; + + template + double + JumpFunction::jump (const Point &p, + const Point &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 + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + + /*virtual*/ double value (const Point &p, + const unsigned int component) const; + + const ExactSolution solution; + }; + + template + double + RightHandSide::value (const Point &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 InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + const SmartPointer preconditioner; + }; + + + template + InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (&preconditioner) + {} + + template + void InverseMatrix::vmult (Vector &dst, + const Vector &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 SchurComplement : public Subscriptor + { + public: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix, Preconditioner> + &A_inverse); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer, + Preconditioner> > A_inverse; + + mutable Vector tmp1, tmp2; + }; + + + + template + SchurComplement:: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix,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 + void SchurComplement::vmult (Vector &dst, + const Vector &src) const + { + system_matrix->block(0,1).vmult (tmp1, src); + A_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); + } + + template + StokesProblem::StokesProblem (const unsigned int degree, FESystem & fe_) + : + degree (degree), + triangulation (Triangulation::maximum_smoothing), + fe (fe_), + dof_handler (triangulation) + {} + + template + void StokesProblem::setup_dofs () + { + A_preconditioner.reset (); + system_matrix.clear (); + + dof_handler.distribute_dofs (fe); + DoFRenumbering::Cuthill_McKee (dof_handler); + + std::vector block_component (dim+1,0); + block_component[dim] = 1; + DoFRenumbering::component_wise (dof_handler, block_component); + + { + constraints.clear (); + std::vector 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(), + constraints, + component_mask); + + /*std::vector boundary_dofs (dof_handler.n_dofs(), false); + + std::vectorboundary_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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + std::vector 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 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 + void StokesProblem::assemble_system () + { + system_matrix=0; + system_rhs=0; + + QGauss quadrature_formula(degree+2); + QGauss quadrature_face(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + + FEFaceValues 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + const JumpFunction jumpfunction; + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + 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::faces_per_cell; ++face_no) + { + typename DoFHandler::face_iterator face= cell->face(face_no); + if (face->at_boundary()==false) + { + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no); + if (neighbor->index() > cell->index()) + { + fe_v_face.reinit (cell, face_no); + + const std::vector > &normals = + fe_v_face.get_normal_vectors (); + const std::vector > &quad_points= + fe_v_face.get_quadrature_points(); + + for (unsigned int q=0; qget_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::type>(new + typename InnerPreconditioner::type()); + A_preconditioner->initialize (system_matrix.block(0,0), + typename + InnerPreconditioner::type::AdditionalData()); + + } + + template + void StokesProblem::solve () + { + const InverseMatrix, + typename InnerPreconditioner::type> + A_inverse (system_matrix.block(0,0), *A_preconditioner); + Vector tmp (solution.block(0).size()); + + { + Vector 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::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 preconditioner; + preconditioner.initialize (system_matrix.block(1,1), + SparseILU::AdditionalData()); + + InverseMatrix,SparseILU > + 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 + void + StokesProblem::output_results (const unsigned int refinement_cycle) + { + const ComponentSelectFunction pressure_mask (dim, dim+1); + const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), + dim+1); + + ExactSolution exactsolution; + + const unsigned int n_active_cells=triangulation.n_active_cells(); + + Vector difference_per_cell (n_active_cells); + + QGauss 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: " + <0)?sqrt(L2_error_pressure):0; + + /*Vector difference(dim+1); + const Point point(.25,0); + VectorTools::point_difference (dof_handler, solution, exactsolution, + difference, point); + std::cout<<"difference at 0.25,0 "< solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (); + + } + + template + void StokesProblem::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 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 + void StokesProblem::divergence_velocity + (const BlockVector &calc_solution, + Vector &output_vector, + const Quadrature &quadrature, bool norm) + { + output_vector = 0; + + FEValues 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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + std::vector 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 fe(FE_Q(degree+1), dim, FE_Q(degree), 1); + deallog << fe.get_name() << ":" << std::endl; + StokesProblem<2> flow_problem(degree, fe); + flow_problem.run (); + } + { + degree=1; + FESystem<2> fe(FE_Q(degree+1), dim, FE_Q_DG0(degree), 1); + deallog << fe.get_name() << ":" << std::endl; + StokesProblem<2> flow_problem(degree, fe); + flow_problem.run (); + } + { + degree=2; + FESystem<2> fe(FE_Q(degree+1), dim, FE_Q_DG0(degree), 1); + deallog << fe.get_name() << ":" << std::endl; + StokesProblem<2> flow_problem(degree, fe); + flow_problem.run (); + } + + return 0; +} diff --git a/tests/fe/fe_q_dg0/cmp/generic b/tests/fe/fe_q_dg0/cmp/generic new file mode 100644 index 0000000000..4d730cf84b --- /dev/null +++ b/tests/fe/fe_q_dg0/cmp/generic @@ -0,0 +1,22 @@ + +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 diff --git a/tests/fe/fe_restriction_common.h b/tests/fe/fe_restriction_common.h index 1c92e184a1..2ee7691c1c 100644 --- a/tests/fe/fe_restriction_common.h +++ b/tests/fe/fe_restriction_common.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include diff --git a/tests/fe/fe_restriction_q_dg0.cc b/tests/fe/fe_restriction_q_dg0.cc new file mode 100644 index 0000000000..03abeeaaa0 --- /dev/null +++ b/tests/fe/fe_restriction_q_dg0.cc @@ -0,0 +1,34 @@ +//---------------------------------------------------------------------- +// $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); +} diff --git a/tests/fe/fe_restriction_q_dg0/cmp/generic b/tests/fe/fe_restriction_q_dg0/cmp/generic new file mode 100644 index 0000000000..54abe855c6 --- /dev/null +++ b/tests/fe/fe_restriction_q_dg0/cmp/generic @@ -0,0 +1,465 @@ + +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 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ diff --git a/tests/fe/fe_support_points_common.h b/tests/fe/fe_support_points_common.h index e5ff06123b..f2cc34af41 100644 --- a/tests/fe/fe_support_points_common.h +++ b/tests/fe/fe_support_points_common.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include diff --git a/tests/fe/fe_support_points_q_dg0.cc b/tests/fe/fe_support_points_q_dg0.cc new file mode 100644 index 0000000000..2e6c6d4f44 --- /dev/null +++ b/tests/fe/fe_support_points_q_dg0.cc @@ -0,0 +1,32 @@ +//---------------------------------------------------------------------- +// $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); +} diff --git a/tests/fe/fe_support_points_q_dg0/cmp/generic b/tests/fe/fe_support_points_q_dg0/cmp/generic new file mode 100644 index 0000000000..ef1829c25b --- /dev/null +++ b/tests/fe/fe_support_points_q_dg0/cmp/generic @@ -0,0 +1,208 @@ + +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 diff --git a/tests/fe/fe_support_points_q_qg0.cc b/tests/fe/fe_support_points_q_qg0.cc new file mode 100644 index 0000000000..b1b3c3719c --- /dev/null +++ b/tests/fe/fe_support_points_q_qg0.cc @@ -0,0 +1,33 @@ +//---------------------------------------------------------------------- +// $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); +} diff --git a/tests/fe/fe_support_points_q_qg0/cmp/generic b/tests/fe/fe_support_points_q_qg0/cmp/generic new file mode 100644 index 0000000000..d2c21ba1a6 --- /dev/null +++ b/tests/fe/fe_support_points_q_qg0/cmp/generic @@ -0,0 +1,213 @@ + +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 diff --git a/tests/fe/injection_common.h b/tests/fe/injection_common.h index 39b0ade345..004818069b 100644 --- a/tests/fe/injection_common.h +++ b/tests/fe/injection_common.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include diff --git a/tests/fe/injection_q_dg0.cc b/tests/fe/injection_q_dg0.cc new file mode 100644 index 0000000000..82b9757f57 --- /dev/null +++ b/tests/fe/injection_q_dg0.cc @@ -0,0 +1,26 @@ +//---------------------------- 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 +void test () +{ + for (unsigned int i=1; i<4; ++i) + for (unsigned int j=i; j<4; ++j) + do_check (FE_Q_DG0(i), FE_Q_DG0(j)); +} diff --git a/tests/fe/interpolate_common.h b/tests/fe/interpolate_common.h index 21c0f9b0a6..bb030ae8d8 100644 --- a/tests/fe/interpolate_common.h +++ b/tests/fe/interpolate_common.h @@ -16,6 +16,7 @@ #include #include #include +#include #include diff --git a/tests/fe/interpolate_q_dg0.cc b/tests/fe/interpolate_q_dg0.cc new file mode 100644 index 0000000000..efd3275049 --- /dev/null +++ b/tests/fe/interpolate_q_dg0.cc @@ -0,0 +1,70 @@ +//---------------------------------------------------------------------- +// $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 +#include + +#include + +// FE_Q_DG0::interpolate(...) + +template +void check_q_dg0(const Function& f, + const unsigned int degree) +{ + FE_Q_DG0 fe(degree); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(1, std::vector(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 > + vectors(fe.get_unit_support_points().size(), Vector(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); +} diff --git a/tests/fe/interpolate_q_dg0/cmp/generic b/tests/fe/interpolate_q_dg0/cmp/generic new file mode 100644 index 0000000000..697d43ac5c --- /dev/null +++ b/tests/fe/interpolate_q_dg0/cmp/generic @@ -0,0 +1,13 @@ + +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 diff --git a/tests/fe/q_dg0_projection_01.cc b/tests/fe/q_dg0_projection_01.cc new file mode 100644 index 0000000000..22d6987350 --- /dev/null +++ b/tests/fe/q_dg0_projection_01.cc @@ -0,0 +1,30 @@ +//---------------------------- 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); +} diff --git a/tests/fe/q_dg0_projection_01/cmp/generic b/tests/fe/q_dg0_projection_01/cmp/generic new file mode 100644 index 0000000000..0b529a682b --- /dev/null +++ b/tests/fe/q_dg0_projection_01/cmp/generic @@ -0,0 +1,143 @@ + +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. +# +# +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. +# +# +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. +# +# +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 diff --git a/tests/fe/q_dg0_projection_01/rt.gpl b/tests/fe/q_dg0_projection_01/rt.gpl new file mode 100644 index 0000000000..6a47c2f34e --- /dev/null +++ b/tests/fe/q_dg0_projection_01/rt.gpl @@ -0,0 +1,184 @@ +#!/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 +# +# Send bugs, suggestions and mods to +# +# +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