]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Test nedelec_2 adapted and a bug in FE_Nedelec::get_subface_interpolation_matrix...
authorMarkus Buerg <buerg@math.tamu.edu>
Mon, 13 Sep 2010 16:00:22 +0000 (16:00 +0000)
committerMarkus Buerg <buerg@math.tamu.edu>
Mon, 13 Sep 2010 16:00:22 +0000 (16:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@21947 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_nedelec.cc
tests/fe/nedelec_2.cc
tests/fe/nedelec_2/cmp/generic

index d0bb20cd7861e4943f8ea3180afca095ece77b91..0a07eb57eae3e5084a57571404f110d993ce5255 100644 (file)
@@ -1,16 +1,3 @@
-//---------------------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 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.
-//
-//---------------------------------------------------------------------------
-
 #include <base/logstream.h>
 #include <base/utilities.h>
 #include <base/quadrature.h>
@@ -1767,8 +1754,10 @@ const
 template <int dim>
 void
 FE_Nedelec<dim>::get_subface_interpolation_matrix(
-                                                 const FiniteElement<dim>&, const unsigned int,
-   FullMatrix<double>&) const
+                                                  const FiniteElement<dim>&,
+                                                  const unsigned int,
+                                                  FullMatrix<double>&)
+const
 {
   Assert (false, ExcNotImplemented ());
 }
@@ -1850,7 +1839,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // interpolation part.
           for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
             {
-              const double weight = 2.0 * edge_quadrature.weight (q_point);
+              const double weight = edge_quadrature.weight (q_point);
 
               for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
                 interpolation_matrix (0, dof)
@@ -1905,7 +1894,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                          const double tmp = 2.0 * quadrature_points[q_point] (0)
                                             + shift[subface];
                          const double weight
-                           = std::sqrt (2.0 * edge_quadrature.weight (q_point));
+                           = std::sqrt (edge_quadrature.weight (q_point));
 
                   for (unsigned int i = 0; i < deg; ++i)
                     assembling_matrix (i, q_point)
@@ -1934,7 +1923,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                         = 2.0 * quadrature_points[q_point] (0)
                           + shift[subface];
                       const double weight
-                        = 2.0 * edge_quadrature.weight (q_point)
+                        = edge_quadrature.weight (q_point)
                           * (this->shape_value_component
                              (dof, Point<dim> (0.0,
                                                quadrature_points[q_point] (0)),
@@ -1978,7 +1967,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                        ++q_point)
                     {
                       const double
-                        weight = 2.0 * edge_quadrature.weight (q_point);
+                        weight = edge_quadrature.weight (q_point);
 
                       for (unsigned int i = 0; i < 2; ++i)
                         for (unsigned int dof = 0; dof < this->dofs_per_face;
@@ -2072,12 +2061,9 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // face.
                       std::vector<std::vector<Point<dim> > >
                         edge_quadrature_points_full_dim
-                        (GeometryInfo<dim>::lines_per_face);
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<Point<dim> > (n_edge_points));
                   
-                      for (unsigned int line = 0;
-                           line < GeometryInfo<dim>::lines_per_face; ++line)
-                        edge_quadrature_points_full_dim.resize (n_edge_points);
-                                  
                       for (unsigned int q_point = 0; q_point < n_edge_points;
                            ++q_point)
                         {
@@ -2097,35 +2083,32 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           0.5, 0.0);
                         }
                   
-                      Vector<double> solution (deg);
-                      Vector<double> system_rhs (deg);
-
-                      for (unsigned int dof = 0; dof < this->dofs_per_face;
-                           ++dof)
-                        {
                                   // Set up the system matrix.
                                   // This can be used for all
                                   // edges.
-                             for (unsigned int q_point = 0;
-                                  q_point < n_edge_points; ++q_point)
-                               {
-                                     const double tmp
-                                       = 2.0 * edge_quadrature_points[q_point] (0);
-                                     const double weight
-                                       = std::sqrt (2.0 * edge_quadrature.weight
-                                                          (q_point));
-
-                              for (unsigned int i = 0; i < deg; ++i)
-                                assembling_matrix (i, q_point)
-                                  = weight
-                                    * lobatto_polynomials_grad[i + 1].value
-                                      (tmp);
-                               }
+                         for (unsigned int q_point = 0;
+                              q_point < n_edge_points; ++q_point)
+                           {
+                                 const double tmp
+                                   = 2.0 * edge_quadrature_points[q_point] (0);
+                                 const double weight
+                                   = std::sqrt (edge_quadrature.weight (q_point));
+
+                          for (unsigned int i = 0; i < deg; ++i)
+                            assembling_matrix (i, q_point)
+                               = weight * lobatto_polynomials_grad[i + 1].value
+                                          (tmp);
+                           }
                      
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.invert (system_matrix);
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.invert (system_matrix);
                      
+                      Vector<double> solution (deg);
+                      Vector<double> system_rhs (deg);
+                      
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
                           for (unsigned int line = 0;
                                line < GeometryInfo<dim>::lines_per_face;
                                ++line)
@@ -2137,7 +2120,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                    q_point < n_edge_points; ++q_point)
                                 {
                                   const double right_hand_side_value
-                                    = std::sqrt (2.0 * edge_quadrature.weight
+                                    = std::sqrt (edge_quadrature.weight
                                                        (q_point))
                                       * (this->shape_value_component
                                          (this->face_to_cell_index (dof, 4),
@@ -2168,9 +2151,44 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                     = solution (i);
                             }
 
-                          assembling_matrix.reinit (deg * this->degree,
-                                                    n_face_points);
-                          system_rhs.reinit (assembling_matrix.m ());
+                      assembling_matrix.reinit (deg * this->degree,
+                                                n_face_points);
+                     
+                      for (unsigned int q_point = 0; q_point < n_face_points;
+                           ++q_point)
+                        {
+                          const Point<dim - 1> quadrature_point
+                                (2.0 * face_quadrature_points[q_point] (0),
+                                 2.0 * face_quadrature_points[q_point] (1));
+                          const double weight
+                            = std::sqrt (face_quadrature.weight (q_point));
+                        
+                          for (unsigned int i = 0; i <= deg; ++i)
+                            {
+                              const double tmp
+                                = weight * legendre_polynomials[i].value
+                                           (quadrature_point (0));
+                              
+                              for (unsigned int j = 0; j < deg; ++j)
+                                   assembling_matrix (i * deg + j, q_point)
+                                     = tmp * lobatto_polynomials[j + 2].value
+                                             (quadrature_point (1));
+                            }
+                        }
+                               
+                      system_matrix.reinit (assembling_matrix.m (),
+                                            assembling_matrix.m ());
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.reinit (system_matrix.m (),
+                                                system_matrix.m ());
+                      system_matrix_inv.invert (system_matrix);
+                      solution.reinit (system_matrix_inv.m ());
+                      system_rhs.reinit (system_matrix_inv.m ());
+                      
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
+                        {
                           system_rhs = 0;
 
                                   // Now we project the remaining
@@ -2204,42 +2222,23 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 1);
                         
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
+                                *= face_quadrature.weight (q_point);
                               
-                              const double weight
-                                = std::sqrt (4.0 * face_quadrature.weight
-                                                   (q_point));
-                        
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp1 = weight * L_i;
-                                  const double tmp2
-                                    = right_hand_side_value * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                      * legendre_polynomials[i].value
+                                        (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
-                                    {
-                                         const double l_j
-                                           = lobatto_polynomials[j + 2].value
-                                             (quadrature_point (1));
-                                         
-                                         assembling_matrix (i * deg + j, q_point)
-                                           = tmp1 * l_j;
-                                      system_rhs (i * deg + j) +=  tmp2 * l_j;
-                                    }
+                                    system_rhs (i * deg + j)
+                                      += tmp
+                                         * lobatto_polynomials[j + 2].value
+                                              (quadrature_point (1));
                                 }
                             }
 
-                          system_matrix.reinit (assembling_matrix.m (),
-                                                assembling_matrix.m ());
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.reinit (system_matrix.m (),
-                                                    system_matrix.m ());
-                          system_matrix_inv.invert (system_matrix);
-                          solution.reinit (system_matrix_inv.m ());
                           system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
@@ -2284,15 +2283,14 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 0);
 
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
+                                *= face_quadrature.weight (q_point);
                         
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp = right_hand_side_value
-                                                     * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                      * legendre_polynomials[i].value
+                                        (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
                                     system_rhs (i * deg + j)
@@ -2341,27 +2339,29 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                   for (unsigned int q_point = 0; q_point < n_edge_points;
                        ++q_point)
                     {
-//                       const double weight
-//                         = 2.0 * edge_quadrature_x.weight (q_point);
+                      const double weight
+                        = edge_quadrature_x.weight (q_point);
 
                       for (unsigned int i = 0; i < 2; ++i)
                         for (unsigned int dof = 0; dof < this->dofs_per_face;
                              ++dof)
                           {
                             interpolation_matrix (i * source_fe.degree, dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (0.5 * (i + 1),
-                                   edge_quadrature_y_points[q_point] (0), 0.0),
-                                  1);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (0.5 * (i + 1),
+                                     edge_quadrature_y_points[q_point] (0), 0.0),
+                                    1);
                             interpolation_matrix
                             ((i + 2) * source_fe.degree, dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (edge_quadrature_x_points[q_point] (0),
-                                   0.5 * i, 0.0), 0);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (edge_quadrature_x_points[q_point] (0),
+                                     0.5 * i, 0.0), 0);
                           }
                     }
                
@@ -2434,22 +2434,16 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // weights on the four
                                   // edges of a face.
                       std::vector<std::vector<double> > edge_quadrature_points
-                        (GeometryInfo<dim>::lines_per_face);
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<double> (n_edge_points));
                       std::vector<std::vector<double> >
                         edge_quadrature_weights
-                          (GeometryInfo<dim>::lines_per_face);
+                          (GeometryInfo<dim>::lines_per_face,
+                           std::vector<double> (n_edge_points));
                       std::vector<std::vector<Point<dim> > >
                         edge_quadrature_points_full_dim
-                          (GeometryInfo<dim>::lines_per_face);
-                  
-                      for (unsigned int line = 0;
-                           line < GeometryInfo<dim>::lines_per_face; ++line)
-                        {
-                             edge_quadrature_points.resize (n_edge_points);
-                          edge_quadrature_points_full_dim.resize
-                            (n_edge_points);
-                          edge_quadrature_weights.resize (n_edge_points);
-                        }
+                          (GeometryInfo<dim>::lines_per_face,
+                           std::vector<Point<dim> > (n_edge_points));
                   
                       for (unsigned int q_point = 0; q_point < n_edge_points;
                            ++q_point)
@@ -2480,46 +2474,41 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                               (edge_quadrature_x_points[q_point] (0), 0.5,
                                0.0);
                           edge_quadrature_weights[0][q_point]
-                            = std::sqrt (2.0 * edge_quadrature_y.weight
-                                               (q_point));
+                            = std::sqrt (edge_quadrature_y.weight (q_point));
                           edge_quadrature_weights[1][q_point]
                             = edge_quadrature_weights[0][q_point];
                           edge_quadrature_weights[2][q_point]
-                            = std::sqrt (2.0 * edge_quadrature_x.weight
-                                               (q_point));
+                            = std::sqrt (edge_quadrature_x.weight (q_point));
                           edge_quadrature_weights[3][q_point]
                             = edge_quadrature_weights[2][q_point];
                         }
                   
-                      Vector<double> system_rhs (system_matrix.m ());
-                      Vector<double> solution (system_rhs.size ());
-
-                      for (unsigned int dof = 0; dof < this->dofs_per_face;
-                          ++dof)
-                        {
                                   // Set up the system matrix.
                                   // This can be used for all
                                   // edges.
-                             for (unsigned int q_point = 0;
-                                  q_point < n_edge_points; ++q_point)
-                               {
-                                      const double tmp
-                                        = 2.0 * edge_quadrature_y_points[q_point] (0);
-                                      const double weight
-                                        = std::sqrt (2.0 * edge_quadrature_y.weight
-                                                           (q_point));
-
-                               for (unsigned int i = 0; i < deg; ++i)
-                                 assembling_matrix (i, q_point)
-                                   = weight
-                                     * lobatto_polynomials_grad[i + 1].value
-                                       (tmp);
-                               }
+                         for (unsigned int q_point = 0;
+                              q_point < n_edge_points; ++q_point)
+                           {
+                                 const double tmp
+                                   = 2.0 * edge_quadrature_y_points[q_point] (0);
+                                 const double weight
+                                   = std::sqrt (edge_quadrature_y.weight (q_point));
+
+                          for (unsigned int i = 0; i < deg; ++i)
+                            assembling_matrix (i, q_point)
+                              = weight
+                                * lobatto_polynomials_grad[i + 1].value (tmp);
+                           }
                      
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.invert (system_matrix);
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.invert (system_matrix);
                      
+                      Vector<double> system_rhs (system_matrix.m ());
+                      Vector<double> solution (system_rhs.size ());
+
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                          ++dof)
                           for (unsigned int line = 0;
                                line < GeometryInfo<dim - 1>::lines_per_cell;
                                ++line)
@@ -2565,9 +2554,45 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // functions. First on the
                                   // horizontal ones, then on
                                   // the vertical ones. 
-                          assembling_matrix.reinit (deg * this->degree,
-                                                    n_face_points);
-                          system_rhs.reinit (assembling_matrix.m ());
+                      assembling_matrix.reinit (deg * this->degree,
+                                                n_face_points);
+                      
+                      for (unsigned int q_point = 0;
+                           q_point < n_face_points; ++q_point)
+                        {
+                          const Point<dim - 1> quadrature_point
+                                (2.0 * face_quadrature_points[q_point] (0)
+                                 - 1.0,
+                                 2.0 * face_quadrature_points[q_point] (1));
+                          const double weight
+                            = std::sqrt (face_quadrature.weight (q_point));
+                        
+                          for (unsigned int i = 0; i <= deg; ++i)
+                            {
+                              const double tmp
+                                = weight * legendre_polynomials[i].value
+                                           (quadrature_point (0));
+                           
+                              for (unsigned int j = 0; j < deg; ++j)
+                                   assembling_matrix (i * deg + j, q_point)
+                                     = tmp * lobatto_polynomials[j + 2].value
+                                             (quadrature_point (1));
+                            }
+                        }
+                            
+                      system_matrix.reinit (assembling_matrix.m (),
+                                            assembling_matrix.m ());
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.reinit (system_matrix.m (),
+                                                system_matrix.m ());
+                      system_matrix_inv.invert (system_matrix);
+                      solution.reinit (system_matrix_inv.m ());
+                      system_rhs.reinit (assembling_matrix.m ());
+                      
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                          ++dof)
+                        {
                           system_rhs = 0;
 
                           for (unsigned int q_point = 0;
@@ -2597,42 +2622,23 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 1);
                         
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
-                                
-                              const double weight
-                                = std::sqrt (4.0 * face_quadrature.weight
-                                                   (q_point));
+                                *= face_quadrature.weight (q_point);
                         
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp1 = weight * L_i;
-                                  const double tmp2 = right_hand_side_value
-                                                      * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                      * legendre_polynomials[i].value
+                                        (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
-                                    {
-                                         const double l_j
-                                           = lobatto_polynomials[j + 2].value
+                                    system_rhs (i * deg + j)
+                                      = tmp
+                                        * lobatto_polynomials[j + 2].value
                                              (quadrature_point (1));
-                                         
-                                         assembling_matrix (i * deg + j, q_point)
-                                           = tmp1 * l_j;
-                                      system_rhs (i * deg + j) +=  tmp2 * l_j;
-                                    }
                                 }
                             }
 
-                          system_matrix.reinit (assembling_matrix.m (),
-                                                assembling_matrix.m ());
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.reinit (system_matrix.m (),
-                                                    system_matrix.m ());
-                          system_matrix_inv.invert (system_matrix);
-                          solution.reinit (system_matrix_inv.m ());
                           system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
@@ -2678,15 +2684,14 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 0);
 
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
+                                *= face_quadrature.weight (q_point);
                         
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
                                   const double tmp
-                                    = right_hand_side_value * L_i;
+                                    = right_hand_side_value
+                                      * legendre_polynomials[i].value
+                                        (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
                                     system_rhs (i * deg + j)
@@ -2735,27 +2740,29 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                   for (unsigned int q_point = 0; q_point < n_edge_points;
                        ++q_point)
                     {
-//                       const double weight
-//                         = 2.0 * edge_quadrature_x.weight (q_point);
+                      const double weight
+                        = edge_quadrature_x.weight (q_point);
 
                       for (unsigned int i = 0; i < 2; ++i)
                         for (unsigned int dof = 0; dof < this->dofs_per_face;
                              ++dof)
                           {
                             interpolation_matrix (i * source_fe.degree, dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (0.5 * i,
-                                   edge_quadrature_y_points[q_point] (0), 0.0),
-                                  1);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (0.5 * i,
+                                     edge_quadrature_y_points[q_point] (0), 0.0),
+                                    1);
                             interpolation_matrix ((i + 2) * source_fe.degree,
                                                   dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (edge_quadrature_x_points[q_point] (0),
-                                   0.5 * (i + 1), 0.0), 0);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (edge_quadrature_x_points[q_point] (0),
+                                     0.5 * (i + 1), 0.0), 0);
                           }
                     }
                
@@ -2824,22 +2831,16 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // edges of a face.
                       std::vector<std::vector<double> >
                         edge_quadrature_points
-                        (GeometryInfo<dim>::lines_per_face);
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<double> (n_edge_points));
                       std::vector<std::vector<double> >
                         edge_quadrature_weights
-                        (GeometryInfo<dim>::lines_per_face);
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<double> (n_edge_points));
                       std::vector<std::vector<Point<dim> > >
                         edge_quadrature_points_full_dim
-                        (GeometryInfo<dim>::lines_per_face);
-                  
-                      for (unsigned int line = 0;
-                           line < GeometryInfo<dim>::lines_per_face; ++line)
-                        {
-                             edge_quadrature_points.resize (n_edge_points);
-                          edge_quadrature_points_full_dim.resize
-                            (n_edge_points);
-                          edge_quadrature_weights.resize (n_edge_points);
-                        }
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<Point<dim> > (n_edge_points));
                   
                       for (unsigned int q_point = 0; q_point < n_edge_points;
                            ++q_point)
@@ -2870,47 +2871,44 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                               (edge_quadrature_x_points[q_point] (0), 1.0,
                                0.0);
                           edge_quadrature_weights[0][q_point]
-                            = std::sqrt (2.0 * edge_quadrature_y.weight
-                                               (q_point));
+                            = std::sqrt (edge_quadrature_y.weight (q_point));
                           edge_quadrature_weights[1][q_point]
                             = edge_quadrature_weights[0][q_point];
                           edge_quadrature_weights[2][q_point]
-                            = std::sqrt (2.0 * edge_quadrature_x.weight
-                                               (q_point));
+                            = std::sqrt (edge_quadrature_x.weight (q_point));
                           edge_quadrature_weights[3][q_point]
                             = edge_quadrature_weights[2][q_point];
                         }
                       
-                      Vector<double> system_rhs (system_matrix.m ());
-                      Vector<double> solution (system_rhs.size ());
-
-                      for (unsigned int dof = 0; dof < this->dofs_per_face;
-                           ++dof)
-                        {
                                   // Set up the system matrix.
                                   // This can be used for all
                                   // edges.
-                             for (unsigned int q_point = 0;
-                                  q_point < n_edge_points; ++q_point)
-                               {
-                                     const double weight
-                                       = std::sqrt (2.0 * edge_quadrature_y.weight
-                                                          (q_point));
-                                     const double tmp
-                                       = 2.0 * edge_quadrature_y_points[q_point] (0)
-                                         - 1.0;
-
-                              for (unsigned int i = 0; i < deg; ++i)
-                                assembling_matrix (i, q_point)
-                                  = weight
-                                    * lobatto_polynomials_grad[i + 1].value
-                                      (tmp);
-                               }
+                         for (unsigned int q_point = 0;
+                              q_point < n_edge_points; ++q_point)
+                           {
+                                 const double weight
+                               = std::sqrt (edge_quadrature_y.weight
+                                            (q_point));
+                                     const double tmp
+                                       = 2.0 * edge_quadrature_y_points[q_point] (0)
+                                         - 1.0;
+
+                          for (unsigned int i = 0; i < deg; ++i)
+                            assembling_matrix (i, q_point)
+                              = weight
+                                * lobatto_polynomials_grad[i + 1].value
+                                  (tmp);
+                        }
                      
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.invert (system_matrix);
+                      assembling_matrix.mTmult (system_matrix,
+                                               assembling_matrix);
+                      system_matrix_inv.invert (system_matrix);
                      
+                      Vector<double> system_rhs (system_matrix.m ());
+                      Vector<double> solution (system_rhs.size ());
+
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
                           for (unsigned int line = 0;
                                line < GeometryInfo<dim - 1>::lines_per_cell;
                                ++line)
@@ -2950,16 +2948,50 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                     = solution (i);
                             }
 
-                          assembling_matrix.reinit (deg * this->degree,
-                                                    n_face_points);
-                          system_rhs.reinit (assembling_matrix.m ());
-                          system_rhs = 0;
+                      assembling_matrix.reinit (deg * this->degree,
+                                                n_face_points);
+                      
+                      for (unsigned int q_point = 0;
+                           q_point < n_face_points; ++q_point)
+                        {
+                           const double weight
+                             = std::sqrt (face_quadrature.weight (q_point));
+                           const Point<dim - 1> quadrature_point
+                             (2.0 * face_quadrature_points[q_point] (0),
+                              2.0 * face_quadrature_points[q_point] (1) - 1.0);
+                           
+                           for (unsigned int i = 0; i <= deg; ++i)
+                             {
+                               const double tmp
+                                 = weight * legendre_polynomials[i].value
+                                            (quadrature_point (0));
+                           
+                               for (unsigned int j = 0; j < deg; ++j)
+                                    assembling_matrix (i * deg + j, q_point)
+                                      = tmp * lobatto_polynomials[j + 2].value
+                                              (quadrature_point (1));
+                             }
+                        }
+                      
+                      system_matrix.reinit (assembling_matrix.m (),
+                                            assembling_matrix.m ());
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.reinit (system_matrix.m (),
+                                                system_matrix.m ());
+                      system_matrix_inv.invert (system_matrix);
+                      solution.reinit (system_matrix_inv.m ());
+                      system_rhs.reinit (assembling_matrix.m ());
+                      system_rhs = 0;
 
                                   // Now we project the remaining
                                   // part on the face shape
                                   // functions. First on the
                                   // horizontal ones, then on
                                   // the vertical ones. 
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
+                        {
                           for (unsigned int q_point = 0;
                                q_point < n_face_points; ++q_point)
                             {
@@ -2986,42 +3018,23 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 1);
                         
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
+                                *= face_quadrature.weight (q_point);
                               
-                              const double weight
-                                = std::sqrt (4.0 * reference_face_quadrature.weight
-                                                   (q_point));
-                        
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp1 = weight * L_i;
-                                  const double tmp2
-                                    = right_hand_side_value * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                      * legendre_polynomials[i].value
+                                        (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
-                                    {
-                                         const double l_j
-                                           = lobatto_polynomials[j + 2].value
-                                             (quadrature_point (1));
-                                         
-                                         assembling_matrix (i * deg + j, q_point)
-                                           = tmp1 * l_j;
-                                      system_rhs (i * deg + j) +=  tmp2 * l_j;
-                                    }
+                                      system_rhs (i * deg + j)
+                                        +=  tmp
+                                            * lobatto_polynomials[j + 2].value
+                                                 (quadrature_point (1));
                                 }
                             }
 
-                          system_matrix.reinit (assembling_matrix.m (),
-                                                assembling_matrix.m ());
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.reinit (system_matrix.m (),
-                                                    system_matrix.m ());
-                          system_matrix_inv.invert (system_matrix);
-                          solution.reinit (system_matrix_inv.m ());
                           system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
@@ -3065,14 +3078,15 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                          (i * source_fe.degree + j,
                                           quadrature_point, 0);
 
-                              right_hand_side_value *= 4.0 * face_quadrature.weight (q_point);
+                              right_hand_side_value *= face_quadrature.weight
+                                                       (q_point);
                        
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp = right_hand_side_value * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                       * legendre_polynomials[i].value
+                                         (quadrature_point (0));
                           
                                   for (unsigned int j = 0; j < deg; ++j)
                                     system_rhs (i * deg + j)
@@ -3115,27 +3129,28 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                   for (unsigned int q_point = 0; q_point < n_edge_points;
                        ++q_point)
                     {
-//                       const double weight
-//                         = 2.0 * edge_quadrature.weight (q_point);
+                      const double weight = edge_quadrature.weight (q_point);
 
                       for (unsigned int i = 0; i < 2; ++i)
                         for (unsigned int dof = 0; dof < this->dofs_per_face;
                              ++dof)
                           {
                             interpolation_matrix (i * source_fe.degree, dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (0.5 * (i + 1),
-                                   edge_quadrature_points[q_point] (0), 0.0),
-                                  1);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (0.5 * (i + 1),
+                                     edge_quadrature_points[q_point] (0), 0.0),
+                                    1);
                             interpolation_matrix ((i + 2) * source_fe.degree,
                                                   dof)
-                              += this->shape_value_component
-                                 (this->face_to_cell_index (dof, 4),
-                                  Point<dim>
-                                  (edge_quadrature_points[q_point] (0),
-                                   0.5 * (i + 1), 0.0), 0);
+                              += weight
+                                 * this->shape_value_component
+                                   (this->face_to_cell_index (dof, 4),
+                                    Point<dim>
+                                    (edge_quadrature_points[q_point] (0),
+                                     0.5 * (i + 1), 0.0), 0);
                           }
                     }
                
@@ -3209,12 +3224,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                   // face.
                       std::vector<std::vector<Point<dim> > >
                         edge_quadrature_points_full_dim
-                        (GeometryInfo<dim>::lines_per_face);
-                  
-                      for (unsigned int line = 0;
-                           line < GeometryInfo<dim>::lines_per_face; ++line)
-                        edge_quadrature_points_full_dim.resize
-                          (n_edge_points);
+                        (GeometryInfo<dim>::lines_per_face,
+                         std::vector<Point<dim> > (n_edge_points));
                   
                       for (unsigned int q_point = 0; q_point < n_edge_points;
                            ++q_point)
@@ -3233,36 +3244,34 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           1.0, 0.0);
                         }
                   
-                      Vector<double> system_rhs (system_matrix.m ());
-                      Vector<double> solution (system_rhs.size ());
-
-                      for (unsigned int dof = 0; dof < this->dofs_per_face;
-                           ++dof)
-                        {
                                   // Set up the system matrix.
                                   // This can be used for all
                                   // edges.
-                             for (unsigned int q_point = 0;
-                                  q_point < n_edge_points; ++q_point)
-                               {
-                                     const double tmp
-                                       = 2.0 * edge_quadrature_points[q_point] (0)
-                                         - 1.0;
-                                     const double weight
-                                       = std::sqrt (2.0 * edge_quadrature.weight
-                                                          (q_point));
-
-                              for (unsigned int i = 0; i < deg; ++i)
-                                assembling_matrix (i, q_point)
-                                  = weight
-                                    * lobatto_polynomials_grad[i + 1].value
-                                      (tmp);
-                               }
+                         for (unsigned int q_point = 0;
+                           q_point < n_edge_points; ++q_point)
+                       {
+                                     const double tmp
+                                       = 2.0 * edge_quadrature_points[q_point] (0)
+                                         - 1.0;
+                                     const double weight
+                                   = std::sqrt (edge_quadrature.weight (q_point));
+
+                          for (unsigned int i = 0; i < deg; ++i)
+                            assembling_matrix (i, q_point)
+                              = weight
+                                * lobatto_polynomials_grad[i + 1].value
+                                  (tmp);
+                               }
                      
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.invert (system_matrix);
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.invert (system_matrix);
                      
+                      Vector<double> system_rhs (system_matrix.m ());
+                      Vector<double> solution (system_rhs.size ());
+
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
                           for (unsigned int line = 0;
                                line < GeometryInfo<dim - 1>::lines_per_cell;
                                ++line)
@@ -3274,8 +3283,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                    q_point < n_edge_points; ++q_point)
                                 {
                                   const double right_hand_side_value
-                                    = std::sqrt (2.0 * edge_quadrature.weight
-                                                       (q_point))
+                                    = std::sqrt (edge_quadrature.weight
+                                                 (q_point))
                                       * (this->shape_value_component
                                          (this->face_to_cell_index (dof, 4),
                                           edge_quadrature_points_full_dim[line][q_point],
@@ -3306,11 +3315,46 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                     = solution (i);
                             }
 
-                          assembling_matrix.reinit (deg * this->degree,
-                                                    n_face_points);
-                          system_rhs.reinit (assembling_matrix.m ());
-                          system_rhs = 0;
+                      assembling_matrix.reinit (deg * this->degree,
+                                                n_face_points);
+                      
+                      for (unsigned int q_point = 0;
+                           q_point < n_face_points; ++q_point)
+                        {
+                          const double weight
+                            = std::sqrt (face_quadrature.weight
+                                         (q_point));
+                          const Point<dim - 1> quadrature_point
+                             (2.0 * face_quadrature_points[q_point] (0) - 1.0,
+                              2.0 * face_quadrature_points[q_point] (1) - 1.0);
+                        
+                          for (unsigned int i = 0; i <= deg; ++i)
+                            {
+                              const double tmp
+                                = weight * legendre_polynomials[i].value
+                                           (quadrature_point (0));
+                           
+                              for (unsigned int j = 0; j < deg; ++j)
+                                   assembling_matrix (i * deg + j, q_point)
+                                     = tmp * lobatto_polynomials[j + 2].value
+                                             (quadrature_point (1));
+                            }
+                        }
+
+                      system_matrix.reinit (assembling_matrix.m (),
+                                            assembling_matrix.m ());
+                      assembling_matrix.mTmult (system_matrix,
+                                                assembling_matrix);
+                      system_matrix_inv.reinit (system_matrix.m (),
+                                                system_matrix.m ());
+                      system_matrix_inv.invert (system_matrix);
+                      solution.reinit (system_matrix.m ());
+                      system_rhs.reinit (assembling_matrix.m ());
+                      system_rhs = 0;
 
+                      for (unsigned int dof = 0; dof < this->dofs_per_face;
+                           ++dof)
+                        {
                                   // Now we project the remaining
                                   // part on the face shape
                                   // functions. First on the
@@ -3343,42 +3387,23 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 1);
                         
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
-                              
-                              const double weight
-                                = std::sqrt (4.0 * face_quadrature.weight
-                                                   (q_point));
+                                *= face_quadrature.weight (q_point);
                         
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
-                                  const double L_i
-                                    = legendre_polynomials[i].value
-                                      (quadrature_point (0));
-                                  const double tmp1 = weight * L_i;
-                                  const double tmp2 = right_hand_side_value
-                                                      * L_i;
+                                  const double tmp
+                                    = right_hand_side_value
+                                       * legendre_polynomials[i].value
+                                         (quadrature_point (0));
                            
                                   for (unsigned int j = 0; j < deg; ++j)
-                                    {
-                                         const double l_j
-                                           = lobatto_polynomials[j + 2].value
-                                             (quadrature_point (1));
-                                         
-                                         assembling_matrix (i * deg + j, q_point)
-                                           = tmp1 * l_j;
-                                      system_rhs (i * deg + j) +=  tmp2 * l_j;
-                                    }
+                                    system_rhs (i * deg + j)
+                                      +=  tmp
+                                          * lobatto_polynomials[j + 2].value
+                                               (quadrature_point (1));
                                 }
                             }
 
-                          system_matrix.reinit (assembling_matrix.m (),
-                                                assembling_matrix.m ());
-                          assembling_matrix.mTmult (system_matrix,
-                                                    assembling_matrix);
-                          system_matrix_inv.reinit (system_matrix.m (),
-                                                    system_matrix.m ());
-                          system_matrix_inv.invert (system_matrix);
-                          solution.reinit (system_matrix.m ());
                           system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
@@ -3424,7 +3449,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix
                                           quadrature_point, 0);
 
                               right_hand_side_value
-                                *= 4.0 * face_quadrature.weight (q_point);
+                                *= face_quadrature.weight (q_point);
                         
                               for (unsigned int i = 0; i <= deg; ++i)
                                 {
@@ -3825,7 +3850,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                   local_dofs[(i + 8) * this->degree] = 0.0;
               }
 
-            if (offset+1 < dim)
+            if (offset + 1 < dim)
               {
                 for (unsigned int i = 0; i < 2; ++i)
                   for (unsigned int j = 0; j < 2; ++j)
@@ -3923,7 +3948,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                     system_rhs = 0;
 
                     if ((((line == 0) || (line == 1) || (line == 4) ||
-                          (line == 5)) && (offset+1 < dim)) ||
+                          (line == 5)) && (offset + 1 < dim)) ||
                         (((line == 2) || (line == 3) || (line == 6) ||
                           (line == 7)) && (offset == 0)) || (line > 7))
                       {
@@ -4010,7 +4035,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                       {
                         case 0:
                           {
-                            if (offset+1 < dim)
+                            if (offset + 1 < dim)
                               {
                                   // Set up the right hand side
                                   // for the horizontal shape
@@ -4114,7 +4139,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
 
                         case 1:
                           {
-                            if (offset+1 < dim)
+                            if (offset + 1 < dim)
                               {
                                   // Set up the right hand side
                                   // for the horizontal shape
@@ -4437,7 +4462,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
 
                         case 4:
                           {
-                            if (offset+1 < dim)
+                            if (offset + 1 < dim)
                               {
                                   // Set up the right hand side
                                   // for the horizontal shape
@@ -4550,7 +4575,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                           }
 
                         default:
-                          if (offset+1 < dim)
+                          if (offset + 1 < dim)
                             {
                                   // Set up the right hand side
                                   // for the horizontal shape
@@ -4733,7 +4758,7 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                 system_rhs.reinit (system_matrix_inv.m ());
                 solution.reinit (system_matrix.m ());
 
-                if (offset+1 < dim)
+                if (offset + 1 < dim)
                   {
                     if (offset == 0)
                       {
@@ -5569,7 +5594,8 @@ const
                                                             + GeometryInfo<dim>::lines_per_cell
                                                             * n_edge_points];
 
-                      for (int i = 2;i < (int) GeometryInfo<dim>::lines_per_face; ++i)
+                      for (unsigned int i = 2;
+                           i < GeometryInfo<dim>::lines_per_face; ++i)
                         for (unsigned int j = 0; j <= deg; ++j)
                           tmp -= local_dofs[edge_indices[face][i]
                                             * this->degree + j]
index eb089f25deceb12c6ce59199794f4b6cce92e89e..e6abc991d3861e7e01ef4b5d73251c9c3c18a866 100644 (file)
 
 template <int dim>
 void
-plot (const Triangulation<dim> &tr)
+plot (const Triangulation<dim> &tr, const unsigned int p)
 {
   deallog << dim << "d, "
           << tr.n_active_cells() << " CELLS" << std::endl;
 
-  FE_Nedelec<dim> fe_ned(1);
+  FE_Nedelec<dim> fe_ned (p);
   
   DoFHandler<dim> dof(tr);
   dof.distribute_dofs(fe_ned);
@@ -87,7 +87,7 @@ plot (const Triangulation<dim> &tr)
 
 template<int dim>
 inline void
-check ()
+check (const unsigned int p)
 {
   Triangulation<dim> tr;
   GridGenerator::hyper_cube(tr, 0., 1.);
@@ -95,12 +95,12 @@ check ()
                                    // first everything on a
                                    // once-refined grid
   tr.refine_global (1);
-  plot (tr);
+  plot (tr, p);
 
                                    // then same with one cell refined
   tr.begin_active()->set_refine_flag ();
   tr.execute_coarsening_and_refinement ();
-  plot (tr);
+  plot (tr, p);
 }
 
 
@@ -113,9 +113,12 @@ main()
   deallog.attach(logfile);
   deallog.depth_console(0);
   deallog.threshold_double(1.e-10);
-  
-  check<2>();
-  check<3>();
+  deallog << "Degree 0:" << std::endl;
+  check<2> (0);
+  check<3> (0);
+  deallog << "Degree 1:" << std::endl;
+  check<2> (1);
+  check<3> (1);
   
   return 0;
 }
index 516283be33f481da46dd21938514fea4fc75869e..8097c289be0fb0750e4e051632f17ff4a592fe77 100644 (file)
@@ -1,4 +1,5 @@
 
+DEAL::Degree 0:
 DEAL::2d, 4 CELLS
 DEAL::  CELL 1.0
 DEAL::, xq=0.00 0.00, f=[4.00 0]
@@ -299,3 +300,304 @@ DEAL::, xq=0.50 0.25 0.50, f=[47.00 8.00 18.00]
 DEAL::, xq=0.25 0.50 0.50, f=[36.00 32.00 30.00]
 DEAL::, xq=0.50 0.50 0.50, f=[36.00 8.00 20.00]
 DEAL::
+DEAL::Degree 1:
+DEAL::2d, 4 CELLS
+DEAL::  CELL 1.0
+DEAL::, xq=0.00 0.00, f=[-9.32 -3.46]
+DEAL::, xq=0.50 0.00, f=[25.32 -6.39]
+DEAL::, xq=0.00 0.50, f=[-12.25 3.46]
+DEAL::, xq=0.50 0.50, f=[36.25 14.39]
+DEAL::
+DEAL::  CELL 1.1
+DEAL::, xq=0.50 0.00, f=[-23.96 -6.39]
+DEAL::, xq=1.00 0.00, f=[79.96 -21.03]
+DEAL::, xq=0.50 0.50, f=[-26.89 14.39]
+DEAL::, xq=1.00 0.50, f=[90.89 69.03]
+DEAL::
+DEAL::  CELL 1.2
+DEAL::, xq=0.00 0.50, f=[-12.25 -35.67]
+DEAL::, xq=0.50 0.50, f=[36.25 -38.60]
+DEAL::, xq=0.00 1.00, f=[-41.53 123.67]
+DEAL::, xq=0.50 1.00, f=[145.53 134.60]
+DEAL::
+DEAL::  CELL 1.3
+DEAL::, xq=0.50 0.50, f=[-26.89 -38.60]
+DEAL::, xq=1.00 0.50, f=[90.89 -50.32]
+DEAL::, xq=0.50 1.00, f=[-53.24 134.60]
+DEAL::, xq=1.00 1.00, f=[189.24 178.32]
+DEAL::
+DEAL::2d, 7 CELLS
+DEAL::  CELL 1.1
+DEAL::, xq=0.50 0.00, f=[-9.32 -3.46]
+DEAL::, xq=1.00 0.00, f=[25.32 -6.39]
+DEAL::, xq=0.50 0.50, f=[-12.25 3.46]
+DEAL::, xq=1.00 0.50, f=[36.25 14.39]
+DEAL::
+DEAL::  CELL 1.2
+DEAL::, xq=0.00 0.50, f=[-26.89 -21.03]
+DEAL::, xq=0.50 0.50, f=[90.89 -23.96]
+DEAL::, xq=0.00 1.00, f=[-29.82 69.03]
+DEAL::, xq=0.50 1.00, f=[101.82 79.96]
+DEAL::
+DEAL::  CELL 1.3
+DEAL::, xq=0.50 0.50, f=[-12.25 -23.96]
+DEAL::, xq=1.00 0.50, f=[36.25 -38.60]
+DEAL::, xq=0.50 1.00, f=[-41.53 79.96]
+DEAL::, xq=1.00 1.00, f=[145.53 134.60]
+DEAL::
+DEAL::  CELL 2.0
+DEAL::, xq=0.00 0.00, f=[-112.34 -100.63]
+DEAL::, xq=0.25 0.00, f=[400.34 -106.49]
+DEAL::, xq=0.00 0.25, f=[-118.20 356.63]
+DEAL::, xq=0.25 0.25, f=[422.20 378.49]
+DEAL::
+DEAL::  CELL 2.1
+DEAL::, xq=0.25 0.00, f=[-141.63 -106.49]
+DEAL::, xq=0.50 0.00, f=[509.63 -5.20]
+DEAL::, xq=0.25 0.25, f=[-147.48 378.49]
+DEAL::, xq=0.50 0.25, f=[531.48 1.73]
+DEAL::
+DEAL::  CELL 2.2
+DEAL::, xq=0.00 0.25, f=[-118.20 -165.05]
+DEAL::, xq=0.25 0.25, f=[422.20 -170.91]
+DEAL::, xq=0.00 0.50, f=[-56.33 597.05]
+DEAL::, xq=0.25 0.50, f=[61.44 618.91]
+DEAL::
+DEAL::  CELL 2.3
+DEAL::, xq=0.25 0.25, f=[-147.48 -170.91]
+DEAL::, xq=0.50 0.25, f=[531.48 -1.73]
+DEAL::, xq=0.25 0.50, f=[2.56 618.91]
+DEAL::, xq=0.50 0.50, f=[120.33 5.20]
+DEAL::
+DEAL::3d, 8 CELLS
+DEAL::  CELL 1.0
+DEAL::, xq=0.00 0.00 0.00, f=[-9.32 -3.46 -26.89]
+DEAL::, xq=0.50 0.00 0.00, f=[25.32 -6.39 -29.82]
+DEAL::, xq=0.00 0.50 0.00, f=[-12.25 3.46 -32.75]
+DEAL::, xq=0.50 0.50 0.00, f=[36.25 14.39 -35.67]
+DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 90.89]
+DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 101.82]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 112.75]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 123.67]
+DEAL::
+DEAL::  CELL 1.1
+DEAL::, xq=0.50 0.00 0.00, f=[-85.45 -6.39 -29.82]
+DEAL::, xq=1.00 0.00 0.00, f=[309.45 -82.53 -100.09]
+DEAL::, xq=0.50 0.50 0.00, f=[-88.38 14.39 -35.67]
+DEAL::, xq=1.00 0.50 0.00, f=[320.38 298.53 -103.02]
+DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 101.82]
+DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 364.09]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 123.67]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 375.02]
+DEAL::
+DEAL::  CELL 1.2
+DEAL::, xq=0.00 0.50 0.00, f=[-12.25 -144.02 -32.75]
+DEAL::, xq=0.50 0.50 0.00, f=[36.25 -146.95 -35.67]
+DEAL::, xq=0.00 1.00 0.00, f=[-149.87 528.02 -161.59]
+DEAL::, xq=0.50 1.00 0.00, f=[549.87 538.95 -164.52]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 112.75]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 123.67]
+DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 593.59]
+DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 604.52]
+DEAL::
+DEAL::  CELL 1.3
+DEAL::, xq=0.50 0.50 0.00, f=[-88.38 -146.95 -35.67]
+DEAL::, xq=1.00 0.50 0.00, f=[320.38 -205.51 -103.02]
+DEAL::, xq=0.50 1.00 0.00, f=[-208.44 538.95 -164.52]
+DEAL::, xq=1.00 1.00 0.00, f=[768.44 757.51 -217.22]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 123.67]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 375.02]
+DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 604.52]
+DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 801.22]
+DEAL::
+DEAL::  CELL 1.4
+DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 -264.07]
+DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 -267.00]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 -269.93]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 -272.86]
+DEAL::, xq=0.00 0.00 1.00, f=[-258.22 -252.36 976.07]
+DEAL::, xq=0.50 0.00 1.00, f=[954.22 -255.29 987.00]
+DEAL::, xq=0.00 0.50 1.00, f=[-261.15 932.36 997.93]
+DEAL::, xq=0.50 0.50 1.00, f=[965.15 943.29 1008.86]
+DEAL::
+DEAL::  CELL 1.5
+DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 -267.00]
+DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 -322.64]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 -272.86]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 -325.57]
+DEAL::, xq=0.50 0.00 1.00, f=[-316.78 -255.29 987.00]
+DEAL::, xq=1.00 0.00 1.00, f=[1172.78 -313.85 1194.64]
+DEAL::, xq=0.50 0.50 1.00, f=[-319.71 943.29 1008.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1183.71 1161.85 1205.57]
+DEAL::
+DEAL::  CELL 1.6
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 -269.93]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 -272.86]
+DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 -369.49]
+DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 -372.42]
+DEAL::, xq=0.00 0.50 1.00, f=[-261.15 -360.70 997.93]
+DEAL::, xq=0.50 0.50 1.00, f=[965.15 -363.63 1008.86]
+DEAL::, xq=0.00 1.00 1.00, f=[-366.56 1336.70 1369.49]
+DEAL::, xq=0.50 1.00 1.00, f=[1358.56 1347.63 1380.42]
+DEAL::
+DEAL::  CELL 1.7
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 -272.86]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 -325.57]
+DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 -372.42]
+DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 -413.41]
+DEAL::, xq=0.50 0.50 1.00, f=[-319.71 -363.63 1008.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1183.71 -407.56 1205.57]
+DEAL::, xq=0.50 1.00 1.00, f=[-410.48 1347.63 1380.42]
+DEAL::, xq=1.00 1.00 1.00, f=[1522.48 1511.56 1533.41]
+DEAL::
+DEAL::3d, 15 CELLS
+DEAL::  CELL 1.1
+DEAL::, xq=0.50 0.00 0.00, f=[-9.32 -3.46 -26.89]
+DEAL::, xq=1.00 0.00 0.00, f=[25.32 -6.39 -29.82]
+DEAL::, xq=0.50 0.50 0.00, f=[-12.25 3.46 -32.75]
+DEAL::, xq=1.00 0.50 0.00, f=[36.25 14.39 -35.67]
+DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 90.89]
+DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 101.82]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 112.75]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 123.67]
+DEAL::
+DEAL::  CELL 1.2
+DEAL::, xq=0.00 0.50 0.00, f=[-88.38 -82.53 -105.95]
+DEAL::, xq=0.50 0.50 0.00, f=[320.38 -85.45 -32.75]
+DEAL::, xq=0.00 1.00 0.00, f=[-91.31 298.53 -108.88]
+DEAL::, xq=0.50 1.00 0.00, f=[331.31 309.45 -111.81]
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 385.95]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 112.75]
+DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 396.88]
+DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 407.81]
+DEAL::
+DEAL::  CELL 1.3
+DEAL::, xq=0.50 0.50 0.00, f=[-12.25 -85.45 -32.75]
+DEAL::, xq=1.00 0.50 0.00, f=[36.25 -158.66 -35.67]
+DEAL::, xq=0.50 1.00 0.00, f=[-161.59 309.45 -111.81]
+DEAL::, xq=1.00 1.00 0.00, f=[593.59 582.66 -170.37]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 112.75]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 123.67]
+DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 407.81]
+DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 626.37]
+DEAL::
+DEAL::  CELL 1.4
+DEAL::, xq=0.00 0.00 0.50, f=[-208.44 -205.51 -223.08]
+DEAL::, xq=0.50 0.00 0.50, f=[768.44 -15.18 -226.01]
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 757.51 -228.94]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 47.18 -231.86]
+DEAL::, xq=0.00 0.00 1.00, f=[-217.22 -211.37 823.08]
+DEAL::, xq=0.50 0.00 1.00, f=[801.22 -214.29 834.01]
+DEAL::, xq=0.00 0.50 1.00, f=[-220.15 779.37 844.94]
+DEAL::, xq=0.50 0.50 1.00, f=[812.15 790.29 855.86]
+DEAL::
+DEAL::  CELL 1.5
+DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 -226.01]
+DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 -287.50]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 -231.86]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 -290.43]
+DEAL::, xq=0.50 0.00 1.00, f=[-281.64 -214.29 834.01]
+DEAL::, xq=1.00 0.00 1.00, f=[1041.64 -278.72 1063.50]
+DEAL::, xq=0.50 0.50 1.00, f=[-284.57 790.29 855.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1052.57 1030.72 1074.43]
+DEAL::
+DEAL::  CELL 1.6
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 -228.94]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 -231.86]
+DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 -334.35]
+DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 -337.28]
+DEAL::, xq=0.00 0.50 1.00, f=[-220.15 -325.57 844.94]
+DEAL::, xq=0.50 0.50 1.00, f=[812.15 -328.49 855.86]
+DEAL::, xq=0.00 1.00 1.00, f=[-331.42 1205.57 1238.35]
+DEAL::, xq=0.50 1.00 1.00, f=[1227.42 1216.49 1249.28]
+DEAL::
+DEAL::  CELL 1.7
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 -231.86]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 -290.43]
+DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 -337.28]
+DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 -378.27]
+DEAL::, xq=0.50 0.50 1.00, f=[-284.57 -328.49 855.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1052.57 -372.42 1074.43]
+DEAL::, xq=0.50 1.00 1.00, f=[-375.35 1216.49 1249.28]
+DEAL::, xq=1.00 1.00 1.00, f=[1391.35 1380.42 1402.27]
+DEAL::
+DEAL::  CELL 2.0
+DEAL::, xq=0.00 0.00 0.00, f=[-826.83 -815.11 -861.96]
+DEAL::, xq=0.25 0.00 0.00, f=[3066.83 -820.97 -867.82]
+DEAL::, xq=0.00 0.25 0.00, f=[-832.68 3023.11 -873.68]
+DEAL::, xq=0.25 0.25 0.00, f=[3088.68 3044.97 -879.53]
+DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 3197.96]
+DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 3219.82]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 3241.68]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 3263.53]
+DEAL::
+DEAL::  CELL 2.1
+DEAL::, xq=0.25 0.00 0.00, f=[-979.09 -820.97 -867.82]
+DEAL::, xq=0.50 0.00 0.00, f=[3635.09 -31.84 -115.22]
+DEAL::, xq=0.25 0.25 0.00, f=[-984.95 3044.97 -879.53]
+DEAL::, xq=0.50 0.25 0.00, f=[3656.95 28.38 -49.81]
+DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 3219.82]
+DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 120.33]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 3263.53]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 51.46]
+DEAL::
+DEAL::  CELL 2.2
+DEAL::, xq=0.00 0.25 0.00, f=[-832.68 -1096.22 -873.68]
+DEAL::, xq=0.25 0.25 0.00, f=[3088.68 -1102.08 -879.53]
+DEAL::, xq=0.00 0.50 0.00, f=[-394.95 4072.22 -83.31]
+DEAL::, xq=0.25 0.50 0.00, f=[422.57 4094.08 134.29]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 3241.68]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 3263.53]
+DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 117.36]
+DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.63]
+DEAL::
+DEAL::  CELL 2.3
+DEAL::, xq=0.25 0.25 0.00, f=[-984.95 -1102.08 -879.53]
+DEAL::, xq=0.50 0.25 0.00, f=[3656.95 43.56 -49.81]
+DEAL::, xq=0.25 0.50 0.00, f=[-190.57 4094.08 134.29]
+DEAL::, xq=0.50 0.50 0.00, f=[626.95 -40.09 -141.87]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 3263.53]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 51.46]
+DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.63]
+DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 149.12]
+DEAL::
+DEAL::  CELL 2.4
+DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 -1336.33]
+DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 -1342.19]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 -1348.05]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 -1353.90]
+DEAL::, xq=0.00 0.00 0.50, f=[-380.71 -927.78 4968.33]
+DEAL::, xq=0.25 0.00 0.50, f=[452.27 16.47 4990.19]
+DEAL::, xq=0.00 0.25 0.50, f=[571.27 998.27 5012.05]
+DEAL::, xq=0.25 0.25 0.50, f=[-558.41 -20.16 5033.90]
+DEAL::
+DEAL::  CELL 2.5
+DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 -1342.19]
+DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 -56.33]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 -1353.90]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 -24.49]
+DEAL::, xq=0.25 0.00 0.50, f=[560.05 16.47 4990.19]
+DEAL::, xq=0.50 0.00 0.50, f=[488.39 104.59 179.22]
+DEAL::, xq=0.25 0.25 0.50, f=[59.35 -20.16 5033.90]
+DEAL::, xq=0.50 0.25 0.50, f=[149.04 -103.77 76.78]
+DEAL::
+DEAL::  CELL 2.6
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 -1348.05]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 -1353.90]
+DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 386.73]
+DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.78]
+DEAL::, xq=0.00 0.25 0.50, f=[571.27 -446.27 5012.05]
+DEAL::, xq=0.25 0.25 0.50, f=[-558.41 7.31 5033.90]
+DEAL::, xq=0.00 0.50 0.50, f=[-448.24 1479.78 139.22]
+DEAL::, xq=0.25 0.50 0.50, f=[480.14 -29.31 198.13]
+DEAL::
+DEAL::  CELL 2.7
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 -1353.90]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 -24.49]
+DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.78]
+DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 76.37]
+DEAL::, xq=0.25 0.25 0.50, f=[59.35 7.31 5033.90]
+DEAL::, xq=0.50 0.25 0.50, f=[149.04 31.59 76.78]
+DEAL::, xq=0.25 0.50 0.50, f=[248.05 -29.31 198.13]
+DEAL::, xq=0.50 0.50 0.50, f=[248.05 31.59 76.37]
+DEAL::

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.