]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement adaptivity. For this, need to add a few functions to FE_FaceQ like face...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 7 Aug 2013 13:13:39 +0000 (13:13 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 7 Aug 2013 13:13:39 +0000 (13:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@30246 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-51/step-51.cc
deal.II/include/deal.II/fe/fe_face.h
deal.II/include/deal.II/fe/fe_q_base.h
deal.II/source/fe/fe_face.cc
deal.II/source/fe/fe_q_base.cc
tests/fe/shapes_faceq/cmp/generic

index 281dbfefdd9330f90272104b8f854d9d2a696bdd..d67bbfb2483af7f978f98caa69304fed29beaa46 100644 (file)
 #include <deal.II/lac/solver_gmres.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
 #include <deal.II/numerics/matrix_tools.h>
 #include <deal.II/numerics/data_out.h>
 
@@ -270,7 +272,7 @@ private:
   void assemble_system (const bool reconstruct_trace = false);
   void solve ();
   void postprocess ();
-  void refine_mesh ();
+  void refine_grid (const unsigned int cylce);
   void output_results (const unsigned int cycle);
 
   Triangulation<dim>   triangulation;
@@ -652,8 +654,8 @@ Step51<dim>::postprocess()
   const unsigned int n_q_points = quadrature.size();
   std::vector<double> u_values(n_q_points);
   std::vector<Tensor<1,dim> > u_gradients(n_q_points);
-  FEValuesExtractors::Vector gradients(0);
-  FEValuesExtractors::Scalar values(dim);
+  FEValuesExtractors::Vector fluxes(0);
+  FEValuesExtractors::Scalar scalar(dim);
   FEValues<dim> fe_values_local(mapping, fe_local, quadrature, update_values);
   FullMatrix<double> cell_matrix(fe_u_post.dofs_per_cell,
                                  fe_u_post.dofs_per_cell);
@@ -669,8 +671,8 @@ Step51<dim>::postprocess()
       fe_values.reinit(cell);
       fe_values_local.reinit(cell_loc);
 
-      fe_values_local[values].get_function_values(solution_local, u_values);
-      fe_values_local[gradients].get_function_values(solution_local, u_gradients);
+      fe_values_local[scalar].get_function_values(solution_local, u_values);
+      fe_values_local[fluxes].get_function_values(solution_local, u_gradients);
       for (unsigned int i=1; i<fe_u_post.dofs_per_cell; ++i)
         {
           for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
@@ -757,31 +759,79 @@ void Step51<dim>::output_results (const unsigned int cycle)
 
 
 
-
 template <int dim>
-void Step51<dim>::run ()
+void Step51<dim>::refine_grid (const unsigned int cycle)
 {
   const bool do_cube = true;
-  if (!do_cube)
+  if (cycle == 0)
     {
-      GridGenerator::hyper_ball (triangulation);
-      static const HyperBallBoundary<dim> boundary;
-      triangulation.set_boundary(0, boundary);
-      triangulation.refine_global(6-2*dim);
+      if (!do_cube)
+        {
+          GridGenerator::hyper_ball (triangulation);
+          static const HyperBallBoundary<dim> boundary;
+          triangulation.set_boundary(0, boundary);
+          triangulation.refine_global(6-2*dim);
+        }
+      else
+        GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1);
     }
+  else
+    switch (refinement_mode)
+      {
+      case global_refinement:
+        {
+          if (do_cube)
+            {
+              triangulation.clear();
+              GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+              triangulation.refine_global(3-dim+cycle/2);
+            }
+          else
+            triangulation.refine_global (1);
+          break;
+        }
+
+      case adaptive_refinement:
+      {
+        Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
+        FEValuesExtractors::Scalar scalar(dim);
+        typename FunctionMap<dim>::type neumann_boundary;
+        KellyErrorEstimator<dim>::estimate (dof_handler_local,
+                                            QGauss<dim-1>(3),
+                                            neumann_boundary,
+                                            solution_local,
+                                            estimated_error_per_cell,
+                                            fe_local.component_mask(scalar));
+
+        GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                         estimated_error_per_cell,
+                                                         0.3, 0.);
+
+        triangulation.execute_coarsening_and_refinement ();
+
+        break;
+      }
+
+      default:
+      {
+        Assert (false, ExcNotImplemented());
+      }
+      }
+  }
+
+
+
+
+
+template <int dim>
+void Step51<dim>::run ()
+{
   for (unsigned int cycle=0; cycle<10; ++cycle)
     {
       std::cout << "Cycle " << cycle << ':' << std::endl;
-
-      if (do_cube)
-        {
-          triangulation.clear();
-          GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
-          triangulation.refine_global(3-dim+cycle/2);
-        }
-      else triangulation.refine_global(1);
-
+      
+      refine_grid (cycle);
       setup_system ();
       assemble_system (false);
       solve ();
index a091042adb7f7747f367d951ecec237756d8b567..3047f7a06afa4d2fc290183a12f3bd6daa7c862f 100644 (file)
@@ -25,9 +25,12 @@ DEAL_II_NAMESPACE_OPEN
 
 
 /**
- * A finite element, which is a tensor product polynomial on each face
- * and undefined in the interior of the cells. The basis functions on
- * the faces are from Polynomials::LagrangeEquidistant
+ * A finite element, which is a tensor product polynomial on each face and
+ * undefined in the interior of the cells. The basis functions on the faces
+ * are Lagrange polynomials based on the support points of the
+ * (dim-1)-dimensional Gauss--Lobatto quadrature rule. For element degree one
+ * and two, the polynomials hence correspond to the usual Lagrange polynomials
+ * on equidistant points.
  *
  * This finite element is the trace space of FE_RaviartThomas on the
  * faces and serves in hybridized methods.
@@ -39,13 +42,9 @@ DEAL_II_NAMESPACE_OPEN
  * element space, but all shape function values extracted will equal
  * to zero.
  *
- * @todo Polynomials::LagrangeEquidistant should be and will be
- * replaced by Polynomials::LagrangeGaussLobatto as soon as such a
- * polynomial set exists.
- *
  * @ingroup fe
- * @author Guido Kanschat
- * @date 2009, 2011
+ * @author Guido Kanschat, Martin Kronbichler
+ * @date 2009, 2011, 2013
  */
 template <int dim, int spacedim=dim>
 class FE_FaceQ : public FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim>
@@ -67,6 +66,31 @@ public:
    */
   virtual std::string get_name () const;
 
+  /**
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+   * element only provides interpolation matrices for elements of the same
+   * type and FE_Nothing. For all other elements, an exception of type
+   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+   */
+  virtual void
+  get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                                 FullMatrix<double>       &matrix) const;
+
+  /**
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+   * element only provides interpolation matrices for elements of the same
+   * type and FE_Nothing. For all other elements, an exception of type
+   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+   */
+  virtual void
+  get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                                    const unsigned int        subface,
+                                    FullMatrix<double>       &matrix) const;
+
   /**
    * Check for non-zero values on a face.
    *
@@ -78,6 +102,28 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
+  /**
+   * Return whether this element implements its hanging node constraints in
+   * the new way, which has to be used to make elements "hp compatible".
+   *
+   * For the FE_Q 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;
+
+  /**
+   * Return whether this element dominates the one given as argument when they
+   * meet at a common face, whether it is the other way around, whether
+   * neither dominates, or if either could dominate.
+   *
+   * For a definition of domination, see FiniteElementBase::Domination and in
+   * particular the @ref hp_paper "hp paper".
+   */
+  virtual
+  FiniteElementDomination::Domination
+  compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+
 private:
   /**
    * Return vector with dofs per vertex, line, quad, hex.
index 6f17bd38a4dfadf800549ffdd5c7c8a175aa0a87..0bbeb3bf738b5f19146b0bf16bc390429c1a1c02 100644 (file)
@@ -62,13 +62,11 @@ public:
   /**
    * Return the matrix interpolating from a face of of one element to the face
    * of the neighboring element.  The size of the matrix is then
-   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
-   *
-   * Derived elements will have to implement this function. They may only
-   * provide interpolation matrices for certain source finite elements, for
-   * example those from the same family. If they don't implement interpolation
-   * from a given element, then they must throw an exception of type
-   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
+   * FE_Q element family only provides interpolation matrices for elements of
+   * the same type and FE_Nothing. For all other elements, an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is
+   * thrown.
    */
   virtual void
   get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
@@ -77,13 +75,11 @@ public:
   /**
    * Return the matrix interpolating from a face of of one element to the face
    * of the neighboring element.  The size of the matrix is then
-   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
-   *
-   * Derived elements will have to implement this function. They may only
-   * provide interpolation matrices for certain source finite elements, for
-   * example those from the same family. If they don't implement interpolation
-   * from a given element, then they must throw an exception of type
-   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
+   * FE_Q element family only provides interpolation matrices for elements of
+   * the same type and FE_Nothing. For all other elements, an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is
+   * thrown.
    */
   virtual void
   get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
index b8637568fc543241cb4d50e318bcf71ce6d79d64..30c23b4c260c485878f91bf76374aebee823ec70 100644 (file)
 
 #include <deal.II/fe/fe_face.h>
 #include <deal.II/fe/fe_poly_face.templates.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/base/quadrature_lib.h>
 #include <sstream>
 
 DEAL_II_NAMESPACE_OPEN
 
+
+namespace
+{
+  std::vector<Point<1> >
+  get_QGaussLobatto_points (const unsigned int degree)
+  {
+    if (degree > 0)
+      {
+        QGaussLobatto<1> quad(degree+1);
+        return quad.get_points();
+      }
+    else
+      return std::vector<Point<1> >(1, Point<1>(0.5));
+  }
+}
+
 template <int dim, int spacedim>
 FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
   :
   FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim> (
-    TensorProductPolynomials<dim-1>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
+    TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
     FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
     std::vector<bool>(1,true))
 {
@@ -40,22 +58,42 @@ FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
     {
       const double step = 1./this->degree;
       Point<codim> p;
+      std::vector<Point<1> > points = get_QGaussLobatto_points(degree);
 
       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;
+              p(0) = points[ix][0];
               if (codim>1)
-                p(1) = iy * step;
+                p(1) = points[iy][0];
               if (codim>2)
-                p(2) = iz * step;
+                p(2) = points[iz][0];
 
               this->unit_face_support_points[k++] = p;
             }
       AssertDimension (k, this->unit_face_support_points.size());
     }
+
+  // initialize unit support points
+  this->unit_support_points.resize(GeometryInfo<dim>::faces_per_cell*
+                                   this->unit_face_support_points.size());
+  const unsigned int n_face_dofs = this->unit_face_support_points.size();
+  for (unsigned int i=0; i<n_face_dofs; ++i)
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        for (unsigned int e=0, c=0; e<dim; ++e)
+          if (d!=e)
+            {
+              this->unit_support_points[n_face_dofs*2*d+i][e] =
+                this->unit_face_support_points[i][c];
+              this->unit_support_points[n_face_dofs*(2*d+1)+i][e] =
+                this->unit_face_support_points[i][c];
+              this->unit_support_points[n_face_dofs*(2*d+1)+i][d] = 1;
+              ++c;
+            }
+      }
 }
 
 
@@ -86,6 +124,170 @@ FE_FaceQ<dim,spacedim>::get_name () const
 
 
 
+template <int dim, int spacedim>
+void
+FE_FaceQ<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+                               FullMatrix<double>       &interpolation_matrix) const
+{
+  // this function is similar to the respective method in FE_Q
+
+  // this is only implemented, if the source FE is also a FE_FaceQ element
+  AssertThrow ((dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe) != 0),
+               (typename FiniteElement<dim,spacedim>::
+                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 FaceQ element, so we will be able to do the work
+  const FE_FaceQ<dim,spacedim> &source_fe
+    = dynamic_cast<const FE_FaceQ<dim,spacedim>&>(x_source_fe);
+
+  // Make sure that the element for which the DoFs should be constrained is
+  // the one with the higher polynomial degree.  Actually the procedure will
+  // work also if this assertion is not satisfied. But the matrices produced
+  // in that case might lead to problems in the hp procedures, which use this
+  // method.
+  Assert (this->dofs_per_face <= source_fe.dofs_per_face,
+          (typename FiniteElement<dim,spacedim>::
+           ExcInterpolationNotImplemented ()));
+
+  // generate a quadrature with the unit face support points. 
+  Quadrature<dim-1> face_quadrature (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.
+  for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
+    {
+      const Point<dim-1> &p = face_quadrature.point (i);
+
+      for (unsigned int j=0; j<this->dofs_per_face; ++j)
+        {
+          double matrix_entry = this->poly_space.compute_value (j, p);
+
+          // Correct the interpolated value. I.e. if it is close to 1 or 0,
+          // make it exactly 1 or 0. Unfortunately, this is required to avoid
+          // problems with higher order elements.
+          if (std::fabs (matrix_entry - 1.0) < eps)
+            matrix_entry = 1.0;
+          if (std::fabs (matrix_entry) < eps)
+            matrix_entry = 0.0;
+
+          interpolation_matrix(i,j) = matrix_entry;
+        }
+    }
+
+  // make sure that the row sum of each of the matrices is 1 at this
+  // point. this must be so since the shape functions sum up to 1
+  for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
+    {
+      double sum = 0.;
+
+      for (unsigned int i=0; i<this->dofs_per_face; ++i)
+        sum += interpolation_matrix(j,i);
+
+      Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1),
+              ExcInternalError());
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceQ<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+                                  const unsigned int        subface,
+                                  FullMatrix<double>       &interpolation_matrix) const
+{
+  // this function is similar to the respective method in FE_Q
+
+  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));
+
+  // see if source is a FaceQ element
+  if (const FE_FaceQ<dim,spacedim> *source_fe
+      = dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe))
+    {
+
+      // Make sure that the element for which the DoFs should be constrained
+      // is the one with the higher polynomial degree.  Actually the procedure
+      // will work also if this assertion is not satisfied. But the matrices
+      // produced in that case might lead to problems in the hp procedures,
+      // which use this method.
+      Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+              (typename FiniteElement<dim,spacedim>::
+               ExcInterpolationNotImplemented ()));
+
+      // generate a quadrature with the unit face support points. 
+      Quadrature<dim-1> face_quadrature (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.
+      for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
+        {
+          const Point<dim-1> p =
+            GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(i),
+                                                            subface);
+
+          for (unsigned int j=0; j<this->dofs_per_face; ++j)
+            {
+              double matrix_entry = this->poly_space.compute_value (j, p);
+
+              // Correct the interpolated value. I.e. if it is close to 1 or 0,
+              // make it exactly 1 or 0. Unfortunately, this is required to avoid
+              // problems with higher order elements.
+              if (std::fabs (matrix_entry - 1.0) < eps)
+                matrix_entry = 1.0;
+              if (std::fabs (matrix_entry) < eps)
+                matrix_entry = 0.0;
+
+              interpolation_matrix(i,j) = matrix_entry;
+            }
+        }
+
+      // make sure that the row sum of each of the matrices is 1 at this
+      // point. this must be so since the shape functions sum up to 1
+      for (unsigned int j=0; j<source_fe->dofs_per_face; ++j)
+        {
+          double sum = 0.;
+
+          for (unsigned int i=0; i<this->dofs_per_face; ++i)
+            sum += interpolation_matrix(j,i);
+
+          Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1),
+                  ExcInternalError());
+        }
+    }
+  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
+    {
+      // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+    }
+  else
+    AssertThrow (false,(typename FiniteElement<dim,spacedim>::
+                        ExcInterpolationNotImplemented()));
+}
+
+
+
 template <int dim, int spacedim>
 bool
 FE_FaceQ<dim,spacedim>::has_support_on_face (
@@ -96,6 +298,7 @@ FE_FaceQ<dim,spacedim>::has_support_on_face (
 }
 
 
+
 template <int dim, int spacedim>
 std::vector<unsigned int>
 FE_FaceQ<dim,spacedim>::get_dpo_vector (const unsigned int deg)
@@ -109,6 +312,41 @@ FE_FaceQ<dim,spacedim>::get_dpo_vector (const unsigned int deg)
 
 
 
+template <int dim, int spacedim>
+bool
+FE_FaceQ<dim,spacedim>::hp_constraints_are_implemented () const
+{
+  return true;
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_FaceQ<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+  if (const FE_FaceQ<dim,spacedim> *fe_q_other
+      = dynamic_cast<const FE_FaceQ<dim,spacedim>*>(&fe_other))
+    {
+      if (this->degree < fe_q_other->degree)
+        return FiniteElementDomination::this_element_dominates;
+      else if (this->degree == fe_q_other->degree)
+        return FiniteElementDomination::either_element_can_dominate;
+      else
+        return FiniteElementDomination::other_element_dominates;
+    }
+  else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+    {
+      // the FE_Nothing has no degrees of freedom and it is typically used in
+      // a context where we don't require any continuity along the interface
+      return FiniteElementDomination::no_requirements;
+    }
+
+  Assert (false, ExcNotImplemented());
+  return FiniteElementDomination::neither_element_dominates;
+}
+
 // explicit instantiations
 #include "fe_face.inst"
 
index 9cae5f3dc44de3c00b43929e0c76895b4ddfb2a2..2540d1cfdc5c25d909f69b63b708d571433d6d78 100644 (file)
@@ -754,7 +754,7 @@ hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
     {
       // 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
+      // could simply 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
index 73baf864ac8a71d8c7bfab1c085604dfdea6492e..ecf589b23662d8b915d31fcee3feb4d353e1b2a3 100644 (file)
@@ -105,37 +105,37 @@ DEAL:Face2d-FaceQ2::0.50 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ2::
 DEAL:Face2d-FaceQ2::
 DEAL:Face2d-FaceQ3::0.00 0.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.12 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.25 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.38 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.12 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.25 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.38 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::0.00 0.50 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::0.50 0.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.06 1.00 1.00 1.00 1.00 1.44 1.80 0.69 1.06 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.12 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.19 1.00 1.00 1.00 1.00 0.97 1.92 1.13 0.98 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.06 1.00 1.00 1.00 1.00 1.40 1.73 0.81 1.06 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.12 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.19 1.00 1.00 1.00 1.00 0.89 1.91 1.26 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::
-DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.31 1.00 1.00 1.00 1.00 0.98 1.13 1.92 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.38 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.44 1.00 1.00 1.00 1.00 1.06 0.69 1.80 1.44 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.31 1.00 1.00 1.00 1.00 0.94 1.26 1.91 0.89 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.38 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.44 1.00 1.00 1.00 1.00 1.06 0.81 1.73 1.40 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::0.50 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::0.00 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.12 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.25 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.38 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.12 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.25 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.38 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::0.50 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::0.00 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.12 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04
-DEAL:Face2d-FaceQ3::0.25 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94
-DEAL:Face2d-FaceQ3::0.38 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12
+DEAL:Face2d-FaceQ3::0.12 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02
+DEAL:Face2d-FaceQ3::0.25 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88
+DEAL:Face2d-FaceQ3::0.38 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05
 DEAL:Face2d-FaceQ3::0.50 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00
 DEAL:Face2d-FaceQ3::
 DEAL:Face2d-FaceQ3::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.