]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implemented computation of shape function 3rd derivatives in real cell, except in...
authorMaien Hamed <tomaien@hotmail.com>
Sun, 6 Sep 2015 12:19:33 +0000 (14:19 +0200)
committerMaien Hamed <tomaien@hotmail.com>
Tue, 8 Sep 2015 10:09:36 +0000 (12:09 +0200)
include/deal.II/fe/fe_poly.h
include/deal.II/fe/fe_poly.templates.h
include/deal.II/fe/fe_update_flags.h
include/deal.II/fe/fe_values.h
source/fe/fe_poly.cc
source/fe/fe_poly_tensor.cc
source/fe/fe_system.cc
source/fe/fe_values.cc
tests/fe/fe_q_3rd_derivative_divergence_theorem.cc [new file with mode: 0644]
tests/fe/fe_q_3rd_derivative_divergence_theorem.output [new file with mode: 0644]
tests/fe/shapes.h

index 0e2916fd845620c634679ee0eb865c316623109d..0b03e783df6d4d5ca22e2b8ed056460808eebbf6 100644 (file)
@@ -258,16 +258,18 @@ protected:
                                       std::vector<Tensor<1,dim> > (n_q_points));
       }
 
-    // if second derivatives through
-    // finite differencing is required,
-    // then initialize some objects for
-    // that
     if (flags & update_hessians)
       {
         grad_grads.resize (this->dofs_per_cell);
         data->shape_hessians.resize (this->dofs_per_cell,
                                      std::vector<Tensor<2,dim> > (n_q_points));
-        data->untransformed_shape_hessians.resize (n_q_points);
+      }
+
+    if (flags & update_3rd_derivatives)
+      {
+        third_derivatives.resize (this->dofs_per_cell);
+        data->shape_3rd_derivatives.resize (this->dofs_per_cell,
+                                            std::vector<Tensor<3,dim> > (n_q_points));
       }
 
     // next already fill those fields
@@ -277,11 +279,14 @@ protected:
     // unit cell, and need to be
     // transformed when visiting an
     // actual cell
-    if (flags & (update_values | update_gradients | update_hessians))
+    if (flags & (update_values | update_gradients
+                 | update_hessians | update_3rd_derivatives) )
       for (unsigned int i=0; i<n_q_points; ++i)
         {
           poly_space.compute(quadrature.point(i),
-                             values, grads, grad_grads, third_derivatives, fourth_derivatives);
+                             values, grads, grad_grads,
+                             third_derivatives,
+                             fourth_derivatives);
 
           if (flags & update_values)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
@@ -294,6 +299,10 @@ protected:
           if (flags & update_hessians)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_hessians[k][i] = grad_grads[k];
+
+          if (flags & update_3rd_derivatives)
+            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+              data->shape_3rd_derivatives[k][i] = third_derivatives[k];
         }
     return data;
   }
@@ -420,29 +429,34 @@ protected:
     std::vector<std::vector<Tensor<2,dim> > > shape_hessians;
 
     /**
-     * Scratch array to store temporary values during hessian calculations in
-     * actual cells.
+     * Array with shape function third derivatives in quadrature points. There
+     * is one row for each shape function, containing values for each
+     * quadrature point.
+     *
+     * We store the third derivatives in the quadrature points on the unit
+     * cell. We then only have to apply the transformation when visiting an
+     * actual cell.
      */
-    mutable std::vector<Tensor<2,dim> > untransformed_shape_hessians;
+    std::vector<std::vector<Tensor<3,dim> > > shape_3rd_derivatives;
   };
 
   /**
-   * Correct the hessian in the reference cell by subtracting the term corresponding
-   * to the Jacobian gradient for one degree of freedom. The result being given by:
-   * @f[
-   * \frac{\partial^2 \phi_i}{\partial\hat{x}_J\partial\hat{x}_K}
-   * - \frac{\partial \phi_i}{\partial {x}_l}
-   * \left( \frac{\partial^2{x}_l}{\partial\hat{x}_J\partial\hat{x}_K} \right).
-   * @f]
-   * After this correction, the shape hessians are simply a mapping_covariant_gradient
-   * transformation.
+   * Correct the shape third derivatives by subtracting the terms corresponding
+   * to the Jacobian pushed forward gradient and second derivative.
+   *
+   * Before the correction, the third derivatives would be given by
+   * D_{ijkl} = \frac{d^3\phi_i}{d \hat x_J d \hat x_K d \hat x_L} (J_{jJ})^{-1} (J_{kK})^{-1} (J_{lL})^{-1},
+   * where J_{iI}=\frac{d x_i}{d \hat x_I}. After the correction, the correct
+   * third derivative would be given by
+   * \frac{d^3\phi_i}{d x_j d x_k d x_l} = D_{ijkl} - H_{mjl} \frac{d^2 \phi_i}{d x_k d x_m}  - H_{mkl} \frac{d^2 \phi_i}{d x_j d x_m} - H_{mjk} \frac{d^2 \phi_i}{d x_l d x_m} - K_{mjkl} \frac{d \phi_i}{d x_m},
+   * where H_{ijk} = \frac{d^2 x_i}{d \hat x_J d \hat x_K} (J_{jJ})^{-1} (J_{kK})^{-1},
+   * and K_{ijkl} = \frac{d^3 x_i}{d \hat x_J d \hat x_K d \hat x_L} (J_{jJ})^{-1} (J_{kK})^{-1} (J_{lL})^{-1}
    */
   void
-  correct_untransformed_hessians (VectorSlice< std::vector<Tensor<2, dim> > >                       uncorrected_shape_hessians,
-                                  const internal::FEValues::MappingRelatedData<dim,spacedim>       &mapping_data,
-                                  const internal::FEValues::FiniteElementRelatedData<dim,spacedim> &fevalues_data,
-                                  const unsigned int                                                n_q_points,
-                                  const unsigned int                                                dof) const;
+  correct_third_derivatives (internal::FEValues::FiniteElementRelatedData<dim,spacedim>       &output_data,
+                             const internal::FEValues::MappingRelatedData<dim,spacedim>       &mapping_data,
+                             const unsigned int                                                n_q_points,
+                             const unsigned int                                                dof) const;
 
   /**
    * The polynomial space. Its type is given by the template parameter POLY.
index 6050446263676b93b051866e2025b00b89d27bac..bdaa97d2e19109b4f14c5e4df154fc038796d272 100644 (file)
@@ -202,7 +202,12 @@ FE_Poly<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const
     out |= update_gradients | update_covariant_transformation;
   if (flags & update_hessians)
     out |= update_hessians | update_covariant_transformation
-           | update_gradients | update_jacobian_grads;
+           | update_gradients | update_jacobian_pushed_forward_grads;
+  if (flags & update_3rd_derivatives)
+    out |= update_3rd_derivatives | update_covariant_transformation
+           | update_hessians | update_gradients
+           | update_jacobian_pushed_forward_grads
+           | update_jacobian_pushed_forward_2nd_derivatives;
   if (flags & update_cell_normal_vectors)
     out |= update_cell_normal_vectors | update_JxW_values;
 
@@ -251,19 +256,26 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
 
       if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
-            }
-
-          correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
-                                          mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (fe_data.untransformed_shape_hessians,
+          mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<spacedim; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+        {
+          mapping.transform (fe_data.shape_3rd_derivatives[k],
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
@@ -316,22 +328,30 @@ fill_fe_face_values (const Mapping<dim,spacedim>
 
       if (flags & update_hessians)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i+offset] = fe_data.shape_hessians[k][i+offset];
-            }
-
-          correct_untransformed_hessians(VectorSlice< std::vector<Tensor<2,dim> > >
-                                         ( fe_data.untransformed_shape_hessians, offset , quadrature.size()),
-                                         mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (make_slice(fe_data.untransformed_shape_hessians,
+          mapping.transform (make_slice(fe_data.shape_hessians[k],
                                         offset,
                                         quadrature.size()),
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<spacedim; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives)
+        {
+          mapping.transform (make_slice(fe_data.shape_3rd_derivatives[k],
+                                        offset,
+                                        quadrature.size()),
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
@@ -384,47 +404,59 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
 
       if (flags & update_hessians)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i+offset] = fe_data.shape_hessians[k][i+offset];
-            }
-
-          correct_untransformed_hessians(VectorSlice< std::vector<Tensor<2,dim> > >
-                                         (fe_data.untransformed_shape_hessians,
-                                          offset,
-                                          quadrature.size()),
-                                         mapping_data,
-                                         output_data,
-                                         quadrature.size(),
-                                         k);
-
-          mapping.transform (make_slice(fe_data.untransformed_shape_hessians,
+          mapping.transform (make_slice(fe_data.shape_hessians[k],
                                         offset,
                                         quadrature.size()),
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<spacedim; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives)
+        {
+          mapping.transform (make_slice(fe_data.shape_3rd_derivatives[k],
+                                        offset,
+                                        quadrature.size()),
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
 
-
 template <class POLY, int dim, int spacedim>
-void
+inline void
 FE_Poly<POLY,dim,spacedim>::
-correct_untransformed_hessians (VectorSlice< std::vector<Tensor<2, dim> > >                       uncorrected_shape_hessians,
-                                const internal::FEValues::MappingRelatedData<dim,spacedim>       &mapping_data,
-                                const internal::FEValues::FiniteElementRelatedData<dim,spacedim> &fevalues_data,
-                                const unsigned int                                                n_q_points,
-                                const unsigned int                                                dof) const
+correct_third_derivatives (internal::FEValues::FiniteElementRelatedData<dim,spacedim>       &output_data,
+                           const internal::FEValues::MappingRelatedData<dim,spacedim>       &mapping_data,
+                           const unsigned int                                                n_q_points,
+                           const unsigned int                                                dof) const
 {
   for (unsigned int i=0; i<n_q_points; ++i)
-    for (unsigned int j=0; j<dim; ++j)
-      for (unsigned int l=0; l<dim; ++l)
-        for (unsigned int n=0; n<spacedim; ++n)
-          uncorrected_shape_hessians[i][j][l] -= fevalues_data.shape_gradients[dof][i][n]
-                                                 * mapping_data.jacobian_grads[i][n][l][j];
+    for (unsigned int j=0; j<spacedim; ++j)
+      for (unsigned int k=0; k<spacedim; ++k)
+        for (unsigned int l=0; l<spacedim; ++l)
+          for (unsigned int m=0; m<spacedim; ++m)
+            {
+              output_data.shape_3rd_derivatives[dof][i][j][k][l] -=
+                (mapping_data.jacobian_pushed_forward_grads[i][m][j][l] *
+                 output_data.shape_hessians[dof][i][k][m])
+                + (mapping_data.jacobian_pushed_forward_grads[i][m][k][l] *
+                   output_data.shape_hessians[dof][i][j][m])
+                + (mapping_data.jacobian_pushed_forward_grads[i][m][j][k] *
+                   output_data.shape_hessians[dof][i][l][m])
+                + (mapping_data.jacobian_pushed_forward_2nd_derivatives[i][m][j][k][l] *
+                   output_data.shape_gradients[dof][i][m]);
+            }
+
 }
 
 namespace internal
index 67b43b00dc6fae94eddab132848692f2b4d5fba6..305ef0ae15ac0b75280a64778ea883eb68aa95ec 100644 (file)
@@ -111,32 +111,38 @@ enum UpdateFlags
    * the real cell.
    */
   update_hessians = 0x0004,
+  //! Third derivatives of shape functions
+  /**
+   * Compute the third derivatives of the shape functions in coordinates of
+   * the real cell
+   */
+  update_3rd_derivatives = 0x0008,
   //! Outer normal vector, not normalized
   /**
    * Vector product of tangential vectors, yielding a normal vector with a
    * length corresponding to the surface element; may be more efficient than
    * computing both.
    */
-  update_boundary_forms = 0x0008,
+  update_boundary_forms = 0x0010,
   //! Transformed quadrature points
   /**
    * Compute the quadrature points transformed into real cell coordinates.
    */
-  update_quadrature_points = 0x0010,
+  update_quadrature_points = 0x0020,
   //! Transformed quadrature weights
   /**
    * Compute the quadrature weights on the real cell, i.e. the weights of the
    * quadrature rule multiplied with the determinant of the Jacobian of the
    * transformation from reference to real cell.
    */
-  update_JxW_values = 0x0020,
+  update_JxW_values = 0x0040,
   //! Normal vectors
   /**
    * Compute the normal vectors, either for a face or for a cell of
    * codimension one. Setting this flag for any other object will raise an
    * error.
    */
-  update_normal_vectors = 0x0040,
+  update_normal_vectors = 0x0080,
   /**
    * @deprecated Use #update_normal_vectors instead.
    */
@@ -150,86 +156,86 @@ enum UpdateFlags
    * Compute the Jacobian of the transformation from the reference cell to the
    * real cell.
    */
-  update_jacobians = 0x0080,
+  update_jacobians = 0x0100,
   //! Gradient of volume element
   /**
    * Compute the derivatives of the Jacobian of the transformation.
    */
-  update_jacobian_grads = 0x0100,
+  update_jacobian_grads = 0x0200,
   //! Volume element
   /**
    * Compute the inverse Jacobian of the transformation from the reference
    * cell to the real cell.
    */
-  update_inverse_jacobians = 0x0200,
+  update_inverse_jacobians = 0x0400,
   //! Covariant transformation
   /**
    * Compute all values the Mapping needs to perform a contravariant
    * transformation of vectors. For special mappings like MappingCartesian
    * this may be simpler than #update_inverse_jacobians.
    */
-  update_covariant_transformation = 0x0400,
+  update_covariant_transformation = 0x0800,
   //! Contravariant transformation
   /**
    * Compute all values the Mapping needs to perform a contravariant
    * transformation of vectors. For special mappings like MappingCartesian
    * this may be simpler than #update_jacobians.
    */
-  update_contravariant_transformation = 0x0800,
+  update_contravariant_transformation = 0x1000,
   //! Shape function values of transformation
   /**
    * Compute the shape function values of the transformation defined by the
    * Mapping.
    */
-  update_transformation_values = 0x1000,
+  update_transformation_values = 0x2000,
   //! Shape function gradients of transformation
   /**
    * Compute the shape function gradients of the transformation defined by the
    * Mapping.
    */
-  update_transformation_gradients = 0x2000,
+  update_transformation_gradients = 0x4000,
   //! Determinant of the Jacobian
   /**
    * Compute the volume element in each quadrature point.
    */
-  update_volume_elements = 0x4000,
+  update_volume_elements = 0x10000,
   /**
    * @deprecated This flag has no effect.
    */
-  update_support_points = 0x10000,
+  update_support_points = 0x20000,
   //! Jacobian at generalized support points
   /**
    * Update the Jacobian of the mapping in generalized support points.
    */
-  update_support_jacobians = 0x20000,
+  update_support_jacobians = 0x40000,
   //! inverse Jacobian at generalized support points
   /**
    * Update the inverse Jacobian of the mapping in generalized support points.
    */
-  update_support_inverse_jacobians = 0x40000,
+  update_support_inverse_jacobians = 0x80000,
   /**
    * Compute the derivatives of the Jacobian of the transformation pushed
    * forward to the real cell coordinates.
    */
-  update_jacobian_pushed_forward_grads = 0x80000,
+  update_jacobian_pushed_forward_grads = 0x100000,
   /**
    * Compute the second derivatives of the Jacobian of the transformation.
    */
-  update_jacobian_2nd_derivatives = 0x100000,
+  update_jacobian_2nd_derivatives = 0x200000,
   /**
    * Compute the second derivatives of the Jacobian of the transformation
    * pushed forward to the real cell coordinates.
    */
-  update_jacobian_pushed_forward_2nd_derivatives = 0x200000,
+  update_jacobian_pushed_forward_2nd_derivatives = 0x400000,
   /**
    * Compute the third derivatives of the Jacobian of the transformation.
    */
-  update_jacobian_3rd_derivatives = 0x400000,
+  update_jacobian_3rd_derivatives = 0x800000,
   /**
    * Compute the third derivatives of the Jacobian of the transformation
    * pushed forward to the real cell coordinates.
    */
-  update_jacobian_pushed_forward_3rd_derivatives = 0x800000,
+  update_jacobian_pushed_forward_3rd_derivatives = 0x1000000,
   /**
    * @deprecated Update quadrature points
    */
@@ -259,6 +265,7 @@ STREAM &operator << (STREAM &s, UpdateFlags u)
   if (u & update_values)                                  s << "values|";
   if (u & update_gradients)                               s << "gradients|";
   if (u & update_hessians)                                s << "hessians|";
+  if (u & update_3rd_derivatives)                         s << "3rd_derivatives|";
   if (u & update_quadrature_points)                       s << "quadrature_points|";
   if (u & update_JxW_values)                              s << "JxW_values|";
   if (u & update_normal_vectors)                          s << "normal_vectors|";
index 123984957af7aafb1653efb05ae66cef5fc60cdb..62a7fa76a820eea2ba3304bfa3a0534d44085a38 100644 (file)
@@ -163,6 +163,13 @@ namespace FEValuesViews
      */
     typedef dealii::Tensor<2,spacedim> hessian_type;
 
+    /**
+     * A typedef for the type of third derivatives of the view this class
+     * represents. Here, for a scalar component of the finite element, the
+     * Third derivative is a <code>Tensor@<3,dim@></code>.
+     */
+    typedef dealii::Tensor<3,spacedim> third_derivative_type;
+
     /**
      * A structure where for each shape function we pre-compute a bunch of
      * data that will make later accesses much cheaper.
@@ -254,6 +261,20 @@ namespace FEValuesViews
     hessian (const unsigned int shape_function,
              const unsigned int q_point) const;
 
+    /**
+     * Return the tensor of rank 3 of all third derivatives of the vector
+     * component selected by this view, for the shape function and quadrature
+     * point selected by the arguments.
+     *
+     * @note The meaning of the arguments is as documented for the value()
+     * function.
+     *
+     * @dealiiRequiresUpdateFlags{update_third_derivatives}
+     */
+    third_derivative_type
+    third_derivative (const unsigned int shape_function,
+                      const unsigned int q_point) const;
+
     /**
      * Return the values of the selected scalar component of the finite
      * element function characterized by <tt>fe_function</tt> at the
@@ -339,6 +360,29 @@ namespace FEValuesViews
     void get_function_laplacians (const InputVector &fe_function,
                                   std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const;
 
+    /**
+     * Return the third derivatives of the selected scalar component of the
+     * finite element function characterized by <tt>fe_function</tt> at the
+     * quadrature points of the cell, face or subface selected the last time
+     * the <tt>reinit</tt> function of the FEValues object was called.
+     *
+     * This function is the equivalent of the
+     * FEValuesBase::get_function_third_derivatives function but it only works
+     * on the selected scalar component.
+     *
+     * The data type stored by the output vector must be what you get when you
+     * multiply the third derivatives of shape functions
+     * (i.e., @p third_derivative_type) times the type used to store the values
+     * of the unknowns $U_j$ of your finite element vector $U$ (represented by
+     * the @p fe_function argument).
+     *
+     * @dealiiRequiresUpdateFlags{update_third_derivatives}
+     */
+    template <class InputVector>
+    void get_function_third_derivatives (const InputVector &fe_function,
+                                         std::vector<typename ProductType<third_derivative_type,
+                                         typename InputVector::value_type>::type> &third_derivatives) const;
+
   private:
     /**
      * A reference to the FEValuesBase object we operate on.
@@ -443,6 +487,13 @@ namespace FEValuesViews
      */
     typedef dealii::Tensor<3,spacedim>          hessian_type;
 
+    /**
+     * A typedef for the type of third derivatives of the view this class
+     * represents. Here, for a set of <code>dim</code> components of the
+     * finite element, the third derivative is a <code>Tensor@<4,dim@></code>.
+     */
+    typedef dealii::Tensor<4,spacedim>          third_derivative_type;
+
     /**
      * A structure where for each shape function we pre-compute a bunch of
      * data that will make later accesses much cheaper.
@@ -608,6 +659,20 @@ namespace FEValuesViews
     hessian (const unsigned int shape_function,
              const unsigned int q_point) const;
 
+    /**
+     * Return the tensor of rank 3 of all third derivatives of
+     * the vector components selected by this view, for the shape function and
+     * quadrature point selected by the arguments.
+     *
+     * @note The meaning of the arguments is as documented for the value()
+     * function.
+     *
+     * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+     */
+    third_derivative_type
+    third_derivative (const unsigned int shape_function,
+                      const unsigned int q_point) const;
+
     /**
      * Return the values of the selected vector components of the finite
      * element function characterized by <tt>fe_function</tt> at the
@@ -765,6 +830,29 @@ namespace FEValuesViews
     void get_function_laplacians (const InputVector &fe_function,
                                   std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const;
 
+    /**
+     * Return the third derivatives of the selected scalar component of the
+     * finite element function characterized by <tt>fe_function</tt> at the
+     * quadrature points of the cell, face or subface selected the last time
+     * the <tt>reinit</tt> function of the FEValues object was called.
+     *
+     * This function is the equivalent of the
+     * FEValuesBase::get_function_third_derivatives function but it only works
+     * on the selected scalar component.
+     *
+     * The data type stored by the output vector must be what you get when you
+     * multiply the third derivatives of shape functions
+     * (i.e., @p third_derivative_type) times the type used to store the values
+     * of the unknowns $U_j$ of your finite element vector $U$ (represented by
+     * the @p fe_function argument).
+     *
+     * @dealiiRequiresUpdateFlags{update_third_derivatives}
+     */
+    template <class InputVector>
+    void get_function_third_derivatives (const InputVector &fe_function,
+                                         std::vector<typename ProductType<third_derivative_type,
+                                         typename InputVector::value_type>::type> &third_derivatives) const;
+
   private:
     /**
      * A reference to the FEValuesBase object we operate on.
@@ -1509,14 +1597,14 @@ public:
    * <tt>point_no</tt>th quadrature point with respect to real cell
    * coordinates. If you want to get the derivatives in one of the coordinate
    * directions, use the appropriate function of the Tensor class to extract
-   * one component. Since only a reference to the derivative values is
+   * one component. Since only a reference to the hessian values is
    * returned, there should be no major performance drawback.
    *
    * If the shape function is vector-valued, then this returns the only non-
    * zero component. If the shape function has more than one non-zero
    * component (i.e. it is not primitive), then throw an exception of type
    * ExcShapeFunctionNotPrimitive. In that case, use the
-   * shape_grad_grad_component() function.
+   * shape_hessian_component() function.
    *
    * The same holds for the arguments of this function as for the
    * shape_value() function.
@@ -1548,6 +1636,50 @@ public:
                            const unsigned int point_no,
                            const unsigned int component) const;
 
+  /**
+   * Third derivatives of the <tt>function_no</tt>th shape function at the
+   * <tt>point_no</tt>th quadrature point with respect to real cell
+   * coordinates. If you want to get the 3rd derivatives in one of the coordinate
+   * directions, use the appropriate function of the Tensor class to extract
+   * one component. Since only a reference to the 3rd derivative values is
+   * returned, there should be no major performance drawback.
+   *
+   * If the shape function is vector-valued, then this returns the only non-
+   * zero component. If the shape function has more than one non-zero
+   * component (i.e. it is not primitive), then throw an exception of type
+   * ExcShapeFunctionNotPrimitive. In that case, use the
+   * shape_3rdderivative_component() function.
+   *
+   * The same holds for the arguments of this function as for the
+   * shape_value() function.
+   *
+   * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+   */
+  const Tensor<3,spacedim> &
+  shape_3rd_derivative (const unsigned int function_no,
+                        const unsigned int point_no) const;
+
+  /**
+   * Return one vector component of the third derivative of a shape function at a
+   * quadrature point. If the finite element is scalar, then only component
+   * zero is allowed and the return value equals that of the shape_3rdderivative()
+   * function. If the finite element is vector valued but all shape functions
+   * are primitive (i.e. they are non-zero in only one component), then the
+   * value returned by shape_3rdderivative() equals that of this function for
+   * exactly one component. This function is therefore only of greater
+   * interest if the shape function is not primitive, but then it is necessary
+   * since the other function cannot be used.
+   *
+   * The same holds for the arguments of this function as for the
+   * shape_value_component() function.
+   *
+   * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+   */
+  Tensor<3,spacedim>
+  shape_3rd_derivative_component (const unsigned int function_no,
+                                  const unsigned int point_no,
+                                  const unsigned int component) const;
+
   //@}
   /// @name Access to values of global finite element fields
   //@{
@@ -1991,6 +2123,103 @@ public:
     bool quadrature_points_fastest = false) const;
   //@}
 
+  //@}
+  /// @name Access to third derivatives of global finite element fields
+  //@{
+
+  /**
+   * Compute the tensor of third derivatives of a finite element at the
+   * quadrature points of a cell. This function is the equivalent of the
+   * corresponding get_function_values() function (see there for more
+   * information) but evaluates the finite element field's third derivatives
+   * instead of its value.
+   *
+   * This function may only be used if the finite element in use is a scalar
+   * one, i.e. has only one vector component. There is a corresponding
+   * function of the same name for vector-valued finite elements.
+   *
+   * @param[in] fe_function A vector of values that describes (globally) the
+   * finite element function that this function should evaluate at the
+   * quadrature points of the current cell.
+   *
+   * @param[out] third_derivatives The third derivatives of the function
+   * specified by fe_function at the quadrature points of the current cell.
+   * The third derivatives are computed in real space (as opposed to on the
+   * unit cell).  The object is assumed to already have the correct size. The
+   * data type stored by this output vector must be what you get when you
+   * multiply the third derivatives of shape function times the type used to
+   * store the values of the unknowns $U_j$ of your finite element vector $U$
+   * (represented by the @p fe_function argument).
+   *
+   * @post <code>third_derivatives[q]</code> will contain the third derivatives
+   * of the field described by fe_function at the $q$th quadrature point.
+   * <code>third_derivatives[q][i][j][k]</code> represents the $(i,j,k)$th
+   * component of the 3rd order tensor of third derivatives at quadrature
+   * point $q$.
+   *
+   * @note The actual data type of the input vector may be either a
+   * Vector&lt;T&gt;, BlockVector&lt;T&gt;, or one of the sequential PETSc or
+   * Trilinos vector wrapper classes. It represents a global vector of DoF
+   * values associated with the DofHandler object with which this FEValues
+   * object was last initialized. Alternatively, if the vector argument is of
+   * type IndexSet, then the function is represented as one that is either
+   * zero or one, depending on whether a DoF index is in the set or not.
+   *
+   * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+   */
+  template <class InputVector>
+  void
+  get_function_third_derivatives (const InputVector &fe_function,
+                                  std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const;
+
+  /**
+   * This function does the same as the other get_function_third_derivatives(),
+   * but applied to multi-component (vector-valued) elements. The meaning of
+   * the arguments is as explained there.
+   *
+   * @post <code>third_derivatives[q]</code> is a vector of third derivatives
+   * of the field described by fe_function at the $q$th quadrature point. The
+   * size of the vector accessed by <code>third_derivatives[q]</code> equals
+   * the number of components of the finite element, i.e.
+   * <code>third_derivatives[q][c]</code> returns the third derivative of the
+   * $c$th vector component at the $q$th quadrature point.
+   * Consequently, <code>third_derivatives[q][c][i][j][k]</code> is
+   * the $(i,j,k)$th component of the tensor of third derivatives of the $c$th
+   * vector component of the vector field at quadrature point $q$ of the
+   * current cell.
+   *
+   * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+   */
+  template <class InputVector>
+  void
+  get_function_third_derivatives (const InputVector      &fe_function,
+                                  std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > &third_derivatives,
+                                  bool quadrature_points_fastest = false) const;
+
+  /**
+   * Access to the third derivatives of a function with more flexibility. See
+   * get_function_values() with corresponding arguments.
+   */
+  template <class InputVector>
+  void get_function_third_derivatives (
+    const InputVector &fe_function,
+    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
+    std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const;
+
+  /**
+   * Access to the third derivatives of a function with more flexibility. See
+   * get_function_values() with corresponding arguments.
+   *
+   * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
+   */
+  template <class InputVector>
+  void get_function_third_derivatives (
+    const InputVector &fe_function,
+    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
+    VectorSlice<std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > > third_derivatives,
+    bool quadrature_points_fastest = false) const;
+  //@}
+
   /// @name Geometry of the cell
   //@{
 
@@ -2590,9 +2819,10 @@ public:
    * associated with this cell, you will not be able to call some functions of
    * this class if they need information about degrees of freedom. These
    * functions are, above all, the
-   * <tt>get_function_value/gradients/hessians/laplacians</tt> functions. If
-   * you want to call these functions, you have to call the @p reinit variants
-   * that take iterators into DoFHandler or other DoF handler type objects.
+   * <tt>get_function_value/gradients/hessians/laplacians/third_derivatives</tt>
+   * functions. If you want to call these functions, you have to call the
+   * @p reinit variants that take iterators into DoFHandler or other DoF handler
+   * type objects.
    */
   void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell);
 
@@ -2797,9 +3027,10 @@ public:
    * freedom possibly associated with this cell, you will not be able to call
    * some functions of this class if they need information about degrees of
    * freedom. These functions are, above all, the
-   * <tt>get_function_value/gradients/hessians</tt> functions. If you want to
-   * call these functions, you have to call the @p reinit variants that take
-   * iterators into DoFHandler or other DoF handler type objects.
+   * <tt>get_function_value/gradients/hessians/third_derivatives</tt>
+   * functions. If you want to call these functions, you have to call the
+   * @p reinit variants that take iterators into DoFHandler or other
+   * DoF handler type objects.
    */
   void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                const unsigned int                                         face_no);
@@ -2907,9 +3138,10 @@ public:
    * freedom possibly associated with this cell, you will not be able to call
    * some functions of this class if they need information about degrees of
    * freedom. These functions are, above all, the
-   * <tt>get_function_value/gradients/hessians</tt> functions. If you want to
-   * call these functions, you have to call the @p reinit variants that take
-   * iterators into DoFHandler or other DoF handler type objects.
+   * <tt>get_function_value/gradients/hessians/third_derivatives</tt>
+   * functions. If you want to call these functions, you have to call the
+   * @p reinit variants that take iterators into DoFHandler or other
+   * DoF handler type objects.
    */
   void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                const unsigned int                    face_no,
@@ -3035,7 +3267,7 @@ namespace FEValuesViews
             typename FVB::ExcAccessToUninitializedField("update_hessians"));
 
     // an adaptation of the
-    // FEValuesBase::shape_grad_component
+    // FEValuesBase::shape_hessian_component
     // function except that here we know the
     // component as fixed and we have
     // pre-computed and cached a bunch of
@@ -3048,6 +3280,32 @@ namespace FEValuesViews
 
 
 
+  template <int dim, int spacedim>
+  inline
+  typename Scalar<dim,spacedim>::third_derivative_type
+  Scalar<dim,spacedim>::third_derivative (const unsigned int shape_function,
+                                          const unsigned int q_point) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (shape_function < fe_values.fe->dofs_per_cell,
+            ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+    Assert (fe_values.update_flags & update_3rd_derivatives,
+            typename FVB::ExcAccessToUninitializedField("update_3rd_derivatives"));
+
+    // an adaptation of the
+    // FEValuesBase::shape_3rdderivative_component
+    // function except that here we know the
+    // component as fixed and we have
+    // pre-computed and cached a bunch of
+    // information. See the comments there.
+    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
+      return fe_values.finite_element_output.shape_3rd_derivatives[shape_function_data[shape_function].row_index][q_point];
+    else
+      return third_derivative_type();
+  }
+
+
+
   template <int dim, int spacedim>
   inline
   typename Vector<dim,spacedim>::value_type
@@ -3335,6 +3593,44 @@ namespace FEValuesViews
       }
   }
 
+  template <int dim, int spacedim>
+  inline
+  typename Vector<dim,spacedim>::third_derivative_type
+  Vector<dim,spacedim>::third_derivative (const unsigned int shape_function,
+                                          const unsigned int q_point) const
+  {
+    // this function works like in
+    // the case above
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (shape_function < fe_values.fe->dofs_per_cell,
+            ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+    Assert (fe_values.update_flags & update_3rd_derivatives,
+            typename FVB::ExcAccessToUninitializedField("update_3rd_derivatives"));
+
+    // same as for the scalar case except
+    // that we have one more index
+    const int snc = shape_function_data[shape_function].single_nonzero_component;
+    if (snc == -2)
+      return third_derivative_type();
+    else if (snc != -1)
+      {
+        third_derivative_type return_value;
+        return_value[shape_function_data[shape_function].single_nonzero_component_index]
+          = fe_values.finite_element_output.shape_3rd_derivatives[snc][q_point];
+        return return_value;
+      }
+    else
+      {
+        third_derivative_type return_value;
+        for (unsigned int d=0; d<dim; ++d)
+          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+            return_value[d]
+              = fe_values.finite_element_output.shape_3rd_derivatives[shape_function_data[shape_function].row_index[d]][q_point];
+
+        return return_value;
+      }
+  }
+
 
   namespace
   {
@@ -3964,6 +4260,75 @@ FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
 
 
 
+template <int dim, int spacedim>
+inline
+const Tensor<3,spacedim> &
+FEValuesBase<dim,spacedim>::shape_3rd_derivative (const unsigned int i,
+                                                  const unsigned int j) const
+{
+  Assert (i < fe->dofs_per_cell,
+          ExcIndexRange (i, 0, fe->dofs_per_cell));
+  Assert (this->update_flags & update_hessians,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  Assert (fe->is_primitive (i),
+          ExcShapeFunctionNotPrimitive(i));
+  Assert (i<this->finite_element_output.shape_3rd_derivatives.size(),
+          ExcIndexRange (i, 0, this->finite_element_output.shape_3rd_derivatives.size()));
+  Assert (j<this->finite_element_output.shape_3rd_derivatives[0].size(),
+          ExcIndexRange (j, 0, this->finite_element_output.shape_3rd_derivatives[0].size()));
+
+  // if the entire FE is primitive,
+  // then we can take a short-cut:
+  if (fe->is_primitive())
+    return this->finite_element_output.shape_3rd_derivatives[i][j];
+  else
+    {
+      // otherwise, use the mapping
+      // between shape function
+      // numbers and rows. note that
+      // by the assertions above, we
+      // know that this particular
+      // shape function is primitive,
+      // so we can call
+      // system_to_component_index
+      const unsigned int
+      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+      return this->finite_element_output.shape_3rd_derivatives[row][j];
+    }
+}
+
+
+
+template <int dim, int spacedim>
+inline
+Tensor<3,spacedim>
+FEValuesBase<dim,spacedim>::shape_3rd_derivative_component (const unsigned int i,
+                                                            const unsigned int j,
+                                                            const unsigned int component) const
+{
+  Assert (i < fe->dofs_per_cell,
+          ExcIndexRange (i, 0, fe->dofs_per_cell));
+  Assert (this->update_flags & update_hessians,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  Assert (component < fe->n_components(),
+          ExcIndexRange(component, 0, fe->n_components()));
+
+  // check whether the shape function
+  // is non-zero at all within
+  // this component:
+  if (fe->get_nonzero_components(i)[component] == false)
+    return Tensor<3,spacedim>();
+
+  // look up the right row in the
+  // table and take the data from
+  // there
+  const unsigned int
+  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
+  return this->finite_element_output.shape_3rd_derivatives[row][j];
+}
+
+
+
 template <int dim, int spacedim>
 inline
 const FiniteElement<dim,spacedim> &
index 6013677501b0adb54e3dde1d0439851b9d94724a..7981da494635b73e89d092b52c0b3b0fde572a1d 100644 (file)
@@ -63,19 +63,26 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
 
       if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
-            }
-
-          correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
-                                          mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (fe_data.untransformed_shape_hessians,
+          mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<2; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+        {
+          mapping.transform (fe_data.shape_3rd_derivatives[k],
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
@@ -116,19 +123,26 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
 
       if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
-            }
-
-          correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
-                                          mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (fe_data.untransformed_shape_hessians,
+          mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<3; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+        {
+          mapping.transform (fe_data.shape_3rd_derivatives[k],
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
@@ -170,19 +184,26 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
 
       if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
-            }
-
-          correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
-                                          mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (fe_data.untransformed_shape_hessians,
+          mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<2; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+        {
+          mapping.transform (fe_data.shape_3rd_derivatives[k],
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
@@ -220,19 +241,26 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
 
       if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
         {
-          // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
-          for (unsigned int i=0; i<quadrature.size(); ++i)
-            {
-              fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
-            }
-
-          correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
-                                          mapping_data, output_data, quadrature.size(), k);
-
-          mapping.transform (fe_data.untransformed_shape_hessians,
+          mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
                              mapping_internal,
                              output_data.shape_hessians[k]);
+
+          for (unsigned int i=0; i<quadrature.size(); ++i)
+            for (unsigned int j=0; j<3; ++j)
+              output_data.shape_hessians[k][i] -=
+                mapping_data.jacobian_pushed_forward_grads[i][j]
+                * output_data.shape_gradients[k][i][j];
+        }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+        {
+          mapping.transform (fe_data.shape_3rd_derivatives[k],
+                             mapping_covariant_hessian,
+                             mapping_internal,
+                             output_data.shape_3rd_derivatives[k]);
+
+          correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
         }
     }
 }
index 09ab842d3553d654494553261f58ab786e46151e..8939e9d05d0329d2433fbcd99b18a86667e147f1 100644 (file)
@@ -740,6 +740,17 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
               Assert(false, ExcNotImplemented());
             }
         }
+
+      // third derivatives are not implemented
+      if (flags & update_3rd_derivatives
+          &&
+          ((cell_similarity != CellSimilarity::translation)
+           ||
+           ((mapping_type == mapping_piola) || (mapping_type == mapping_raviart_thomas)
+            || (mapping_type == mapping_nedelec))))
+        {
+          Assert(false, ExcNotImplemented())
+        }
     }
 }
 
@@ -1195,6 +1206,12 @@ fill_fe_face_values (const Mapping<dim,spacedim>
               Assert(false, ExcNotImplemented());
             }
         }
+
+      // third derivatives are not implemented
+      if (flags & update_3rd_derivatives)
+        {
+          Assert(false, ExcNotImplemented())
+        }
     }
 }
 
@@ -1653,6 +1670,12 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
               Assert(false, ExcNotImplemented());
             }
         }
+
+      // third derivatives are not implemented
+      if (flags & update_3rd_derivatives)
+        {
+          Assert(false, ExcNotImplemented())
+        }
     }
 }
 
index 2be7e80d4813d18900ade3c437167c35f31ab72b..afd837eca0bb05816e481768a64c7a4796f55bac 100644 (file)
@@ -1233,6 +1233,11 @@ compute_fill_one_base (const Mapping<dim,spacedim>                      &mapping
                   output_data.shape_hessians[out_index+s][q] =
                     base_data.shape_hessians[in_index+s][q];
 
+              if (base_flags & update_3rd_derivatives)
+                for (unsigned int q=0; q<n_q_points; ++q)
+                  output_data.shape_3rd_derivatives[out_index+s][q] =
+                    base_data.shape_3rd_derivatives[in_index+s][q];
+
             }
         }
 }
@@ -1270,7 +1275,8 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
                           fe_data.update_flags);
 
 
-  if (flags & (update_values | update_gradients | update_hessians))
+  if (flags & (update_values | update_gradients
+               | update_hessians | update_3rd_derivatives ))
     {
       // let base elements update the necessary data
       Threads::TaskGroup<> task_group;
index faa37b8b6e2c4c60e7d0fe78b60baf3b69d04ec6..f5341ebee6c6b58b1b11c817666f2f30bdf75cc7 100644 (file)
@@ -1356,6 +1356,30 @@ namespace FEValuesViews
 
 
 
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Scalar<dim,spacedim>::
+  get_function_third_derivatives (const InputVector &fe_function,
+                                  std::vector<typename ProductType<third_derivative_type,typename InputVector::value_type>::type> &third_derivatives) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (fe_values.update_flags & update_3rd_derivatives,
+            typename FVB::ExcAccessToUninitializedField("update_3rd_derivatives"));
+    Assert (fe_values.present_cell.get() != 0,
+            ExcMessage ("FEValues object is not reinit'ed to any cell"));
+    AssertDimension (fe_function.size(),
+                     fe_values.present_cell->n_dofs_for_dof_handler());
+
+    // get function values of dofs on this cell
+    dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+    fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+    internal::do_function_derivatives<3,dim,spacedim>
+    (dof_values, fe_values.finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
+  }
+
+
+
   template <int dim, int spacedim>
   template <class InputVector>
   void
@@ -1527,6 +1551,29 @@ namespace FEValuesViews
   }
 
 
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::
+  get_function_third_derivatives (const InputVector &fe_function,
+                                  std::vector<typename ProductType<third_derivative_type,typename InputVector::value_type>::type> &third_derivatives) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (fe_values.update_flags & update_3rd_derivatives,
+            typename FVB::ExcAccessToUninitializedField("update_3rd_derivatives"));
+    Assert (fe_values.present_cell.get() != 0,
+            ExcMessage ("FEValues object is not reinit'ed to any cell"));
+    AssertDimension (fe_function.size(),
+                     fe_values.present_cell->n_dofs_for_dof_handler());
+
+    // get function values of dofs on this cell
+    dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+    fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+    internal::do_function_derivatives<3,dim,spacedim>
+    (dof_values, fe_values.finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
+  }
+
+
 
   template <int dim, int spacedim>
   template <class InputVector>
@@ -2023,7 +2070,7 @@ FEValuesBase<dim,spacedim>::TriaCellIterator::message_string
   = ("You have previously called the FEValues::reinit function with a\n"
      "cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,\n"
      "when you do this, you cannot call some functions in the FEValues\n"
-     "class, such as the get_function_values/gradients/hessians\n"
+     "class, such as the get_function_values/gradients/hessians/third_derivatives\n"
      "functions. If you need these functions, then you need to call\n"
      "FEValues::reinit with an iterator type that allows to extract\n"
      "degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.");
@@ -2178,6 +2225,10 @@ namespace internal
       if (flags & update_hessians)
         this->shape_hessians.resize (n_nonzero_shape_components,
                                      std::vector<Tensor<2,spacedim> > (n_quadrature_points));
+
+      if (flags & update_3rd_derivatives)
+        this->shape_3rd_derivatives.resize (n_nonzero_shape_components,
+                                            std::vector<Tensor<3,spacedim> > (n_quadrature_points));
     }
 
 
@@ -2190,6 +2241,7 @@ namespace internal
       return (MemoryConsumption::memory_consumption (shape_values) +
               MemoryConsumption::memory_consumption (shape_gradients) +
               MemoryConsumption::memory_consumption (shape_hessians) +
+              MemoryConsumption::memory_consumption (shape_3rd_derivatives) +
               MemoryConsumption::memory_consumption (shape_function_to_row_table));
     }
   }
@@ -3188,6 +3240,126 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
 
 
 
+template <int dim, int spacedim>
+template <class InputVector>
+void
+FEValuesBase<dim,spacedim>::
+get_function_third_derivatives (const InputVector                &fe_function,
+                                std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const
+{
+  typedef typename InputVector::value_type Number;
+  AssertDimension (fe->n_components(), 1);
+  Assert (this->update_flags & update_3rd_derivatives,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  Assert (present_cell.get() != 0,
+          ExcMessage ("FEValues object is not reinit'ed to any cell"));
+  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
+
+  // get function values of dofs on this cell
+  Vector<Number> dof_values (dofs_per_cell);
+  present_cell->get_interpolated_dof_values(fe_function, dof_values);
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_3rd_derivatives,
+                                    third_derivatives);
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector>
+void FEValuesBase<dim,spacedim>::get_function_third_derivatives (
+  const InputVector &fe_function,
+  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
+  std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const
+{
+  typedef typename InputVector::value_type Number;
+  Assert (this->update_flags & update_3rd_derivatives,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
+  AssertDimension (indices.size(), dofs_per_cell);
+  if (dofs_per_cell <= 100)
+    {
+      Number dof_values[100];
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        dof_values[i] = get_vector_element (fe_function, indices[i]);
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_3rd_derivatives,
+                                        third_derivatives);
+    }
+  else
+    {
+      Vector<Number> dof_values(dofs_per_cell);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        dof_values[i] = get_vector_element (fe_function, indices[i]);
+      internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_3rd_derivatives,
+                                        third_derivatives);
+    }
+}
+
+
+
+
+template <int dim, int spacedim>
+template <class InputVector>
+void
+FEValuesBase<dim,spacedim>::
+get_function_third_derivatives (const InputVector                         &fe_function,
+                                std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > &third_derivatives,
+                                bool quadrature_points_fastest) const
+{
+  typedef typename InputVector::value_type Number;
+  Assert (this->update_flags & update_3rd_derivatives,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  Assert (present_cell.get() != 0,
+          ExcMessage ("FEValues object is not reinit'ed to any cell"));
+  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
+
+  // get function values of dofs on this cell
+  Vector<Number> dof_values (dofs_per_cell);
+  present_cell->get_interpolated_dof_values(fe_function, dof_values);
+  VectorSlice<std::vector<std::vector<Tensor<3,spacedim,Number> > > > third(third_derivatives);
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_3rd_derivatives,
+                                    *fe, this->finite_element_output.shape_function_to_row_table,
+                                    third, quadrature_points_fastest);
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector>
+void FEValuesBase<dim, spacedim>::get_function_third_derivatives (
+  const InputVector &fe_function,
+  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
+  VectorSlice<std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > > third_derivatives,
+  bool quadrature_points_fastest) const
+{
+  typedef typename InputVector::value_type Number;
+  Assert (this->update_flags & update_3rd_derivatives,
+          ExcAccessToUninitializedField("update_3rd_derivatives"));
+  Assert (indices.size() % dofs_per_cell == 0,
+          ExcNotMultiple(indices.size(), dofs_per_cell));
+  if (indices.size() <= 100)
+    {
+      Number dof_values[100];
+      for (unsigned int i=0; i<indices.size(); ++i)
+        dof_values[i] = get_vector_element (fe_function, indices[i]);
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_3rd_derivatives,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
+                                        third_derivatives, quadrature_points_fastest,
+                                        indices.size()/dofs_per_cell);
+    }
+  else
+    {
+      Vector<Number> dof_values(indices.size());
+      for (unsigned int i=0; i<indices.size(); ++i)
+        dof_values[i] = get_vector_element (fe_function, indices[i]);
+      internal::do_function_derivatives(dof_values.begin(),this->finite_element_output.shape_3rd_derivatives,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
+                                        third_derivatives, quadrature_points_fastest,
+                                        indices.size()/dofs_per_cell);
+    }
+}
+
+
+
 template <int dim, int spacedim>
 const typename Triangulation<dim,spacedim>::cell_iterator
 FEValuesBase<dim,spacedim>::get_cell () const
diff --git a/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc
new file mode 100644 (file)
index 0000000..c744720
--- /dev/null
@@ -0,0 +1,176 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check the correctness of fe_values.shape_3rd_derivative for FE_Q by comparing
+// the integral of all shape third derivative components with the flux of the
+// hessian over the boundary by the divergence theorem
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <sstream>
+#include <fstream>
+
+template<int dim>
+Tensor<1,dim> ones ()
+{
+  Tensor<1,dim> result;
+  for (unsigned int i=0; i<dim; ++i)
+    result[i] = 1.0;
+  return result;
+}
+
+template<int dim>
+void test (const Triangulation<dim> &tr,
+           const FiniteElement<dim> &fe,
+           const double tolerance)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  std::stringstream ss;
+
+  deallog << "FE=" << fe.get_name() << std::endl;
+
+  const QGauss<dim> quadrature(6);
+  FEValues<dim> fe_values (fe, quadrature, update_3rd_derivatives
+                           | update_quadrature_points
+                           | update_JxW_values);
+
+  const QGauss<dim-1> face_quadrature(6);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature,
+                                    update_hessians
+                                    | update_quadrature_points
+                                    | update_normal_vectors
+                                    | update_JxW_values);
+
+  for (typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+       cell != dof.end();
+       ++cell)
+    {
+      fe_values.reinit (cell);
+
+      deallog << "Cell nodes:" << std::endl;
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+        {
+          deallog << i << ": ( ";
+          for (unsigned int d=0; d<dim; ++d)
+            deallog << cell->vertex(i)[d] << " ";
+          deallog << ")" << std::endl;
+        }
+
+      bool cell_ok = true;
+
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+        {
+          FEValuesExtractors::Scalar single_component (c);
+
+          for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+            {
+              ss << "component=" << c
+                 << ", dof=" << i
+                 << std::endl;
+
+              Tensor<3,dim> bulk_integral;
+              for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+                {
+                  bulk_integral += fe_values[single_component].third_derivative(i,q) * fe_values.JxW(q);
+
+                  Tensor<3,dim> third_derivative = fe_values[single_component].third_derivative(i,q);
+                }
+
+              Tensor<3,dim> boundary_integral;
+              for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+                {
+                  fe_face_values.reinit(cell,face);
+                  for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
+                    {
+                      Tensor<2,dim> hessian = fe_face_values[single_component].hessian (i,q);
+                      Tensor<3,dim> hessian_normal_outer_prod;
+
+                      outer_product(hessian_normal_outer_prod, hessian, fe_face_values.normal_vector(q));
+                      boundary_integral += hessian_normal_outer_prod * fe_face_values.JxW(q);
+                    }
+                }
+
+              if ((bulk_integral-boundary_integral).norm_square() > tolerance * (bulk_integral.norm() + boundary_integral.norm()))
+                {
+                  deallog << "Failed:" << std::endl;
+                  deallog << ss.str() << std::endl;
+                  deallog << "    bulk integral=" << bulk_integral << std::endl;
+                  deallog << "boundary integral=" << boundary_integral << std::endl;
+                  deallog << "Error! difference between bulk and surface integrals is "
+                          << (bulk_integral-boundary_integral).norm_square()
+                          << " and greater than "
+                          << tolerance * (bulk_integral.norm() + boundary_integral.norm())
+                          << "!\n\n" << std::endl;
+                  cell_ok = false;
+                }
+
+              ss.str("");
+            }
+        }
+
+      deallog << (cell_ok? "OK: cell bulk and boundary integrals match...\n" : "Failed divergence test...\n") << std::endl;
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_ball(const double tolerance)
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const HyperBallBoundary<dim> boundary;
+  tr.set_boundary (0, boundary);
+
+  tr.refine_global(1);
+
+  FE_Q<dim>  fe(3);
+  test(tr, fe, tolerance);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (3);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-7);
+
+  test_hyper_ball<2>(1e-6);
+  test_hyper_ball<3>(1e-6);
+
+  deallog << "done..." << std::endl;
+}
+
diff --git a/tests/fe/fe_q_3rd_derivative_divergence_theorem.output b/tests/fe/fe_q_3rd_derivative_divergence_theorem.output
new file mode 100644 (file)
index 0000000..2a7c032
--- /dev/null
@@ -0,0 +1,760 @@
+
+DEAL::FE=FE_Q<2>(3)
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 -0.707 )
+DEAL::1: ( 0 -1.00 )
+DEAL::2: ( -0.500 -0.500 )
+DEAL::3: ( 0 -0.646 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -1.00 )
+DEAL::1: ( 0.707 -0.707 )
+DEAL::2: ( 0 -0.646 )
+DEAL::3: ( 0.500 -0.500 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.500 -0.500 )
+DEAL::1: ( 0 -0.646 )
+DEAL::2: ( -0.293 -0.293 )
+DEAL::3: ( 0 -0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.646 )
+DEAL::1: ( 0.500 -0.500 )
+DEAL::2: ( 0 -0.293 )
+DEAL::3: ( 0.293 -0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 -0.707 )
+DEAL::1: ( -0.500 -0.500 )
+DEAL::2: ( -1.00 0 )
+DEAL::3: ( -0.646 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.500 -0.500 )
+DEAL::1: ( -0.293 -0.293 )
+DEAL::2: ( -0.646 0 )
+DEAL::3: ( -0.293 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -1.00 0 )
+DEAL::1: ( -0.646 0 )
+DEAL::2: ( -0.707 0.707 )
+DEAL::3: ( -0.500 0.500 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.646 0 )
+DEAL::1: ( -0.293 0 )
+DEAL::2: ( -0.500 0.500 )
+DEAL::3: ( -0.293 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.293 -0.293 )
+DEAL::1: ( 0 -0.293 )
+DEAL::2: ( -0.293 0 )
+DEAL::3: ( 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.293 )
+DEAL::1: ( 0.293 -0.293 )
+DEAL::2: ( 0 0 )
+DEAL::3: ( 0.293 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.293 0 )
+DEAL::1: ( 0 0 )
+DEAL::2: ( -0.293 0.293 )
+DEAL::3: ( 0 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 )
+DEAL::1: ( 0.293 0 )
+DEAL::2: ( 0 0.293 )
+DEAL::3: ( 0.293 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.707 -0.707 )
+DEAL::1: ( 1.00 0 )
+DEAL::2: ( 0.500 -0.500 )
+DEAL::3: ( 0.646 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 1.00 0 )
+DEAL::1: ( 0.707 0.707 )
+DEAL::2: ( 0.646 0 )
+DEAL::3: ( 0.500 0.500 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.500 -0.500 )
+DEAL::1: ( 0.646 0 )
+DEAL::2: ( 0.293 -0.293 )
+DEAL::3: ( 0.293 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.646 0 )
+DEAL::1: ( 0.500 0.500 )
+DEAL::2: ( 0.293 0 )
+DEAL::3: ( 0.293 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 0.707 )
+DEAL::1: ( -0.500 0.500 )
+DEAL::2: ( 0 1.00 )
+DEAL::3: ( 0 0.646 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.500 0.500 )
+DEAL::1: ( -0.293 0.293 )
+DEAL::2: ( 0 0.646 )
+DEAL::3: ( 0 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 1.00 )
+DEAL::1: ( 0 0.646 )
+DEAL::2: ( 0.707 0.707 )
+DEAL::3: ( 0.500 0.500 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0.646 )
+DEAL::1: ( 0 0.293 )
+DEAL::2: ( 0.500 0.500 )
+DEAL::3: ( 0.293 0.293 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::FE=FE_Q<3>(3)
+DEAL::Cell nodes:
+DEAL::0: ( -0.211 -0.211 -0.211 )
+DEAL::1: ( 0 -0.211 -0.211 )
+DEAL::2: ( -0.211 0 -0.211 )
+DEAL::3: ( 0 0 -0.211 )
+DEAL::4: ( -0.211 -0.211 0 )
+DEAL::5: ( 0 -0.211 0 )
+DEAL::6: ( -0.211 0 0 )
+DEAL::7: ( 0 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.211 -0.211 )
+DEAL::1: ( 0.211 -0.211 -0.211 )
+DEAL::2: ( 0 0 -0.211 )
+DEAL::3: ( 0.211 0 -0.211 )
+DEAL::4: ( 0 -0.211 0 )
+DEAL::5: ( 0.211 -0.211 0 )
+DEAL::6: ( 0 0 0 )
+DEAL::7: ( 0.211 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.211 0 -0.211 )
+DEAL::1: ( 0 0 -0.211 )
+DEAL::2: ( -0.211 0.211 -0.211 )
+DEAL::3: ( 0 0.211 -0.211 )
+DEAL::4: ( -0.211 0 0 )
+DEAL::5: ( 0 0 0 )
+DEAL::6: ( -0.211 0.211 0 )
+DEAL::7: ( 0 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 -0.211 )
+DEAL::1: ( 0.211 0 -0.211 )
+DEAL::2: ( 0 0.211 -0.211 )
+DEAL::3: ( 0.211 0.211 -0.211 )
+DEAL::4: ( 0 0 0 )
+DEAL::5: ( 0.211 0 0 )
+DEAL::6: ( 0 0.211 0 )
+DEAL::7: ( 0.211 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.211 -0.211 0 )
+DEAL::1: ( 0 -0.211 0 )
+DEAL::2: ( -0.211 0 0 )
+DEAL::3: ( 0 0 0 )
+DEAL::4: ( -0.211 -0.211 0.211 )
+DEAL::5: ( 0 -0.211 0.211 )
+DEAL::6: ( -0.211 0 0.211 )
+DEAL::7: ( 0 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.211 0 )
+DEAL::1: ( 0.211 -0.211 0 )
+DEAL::2: ( 0 0 0 )
+DEAL::3: ( 0.211 0 0 )
+DEAL::4: ( 0 -0.211 0.211 )
+DEAL::5: ( 0.211 -0.211 0.211 )
+DEAL::6: ( 0 0 0.211 )
+DEAL::7: ( 0.211 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.211 0 0 )
+DEAL::1: ( 0 0 0 )
+DEAL::2: ( -0.211 0.211 0 )
+DEAL::3: ( 0 0.211 0 )
+DEAL::4: ( -0.211 0 0.211 )
+DEAL::5: ( 0 0 0.211 )
+DEAL::6: ( -0.211 0.211 0.211 )
+DEAL::7: ( 0 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 0 )
+DEAL::1: ( 0.211 0 0 )
+DEAL::2: ( 0 0.211 0 )
+DEAL::3: ( 0.211 0.211 0 )
+DEAL::4: ( 0 0 0.211 )
+DEAL::5: ( 0.211 0 0.211 )
+DEAL::6: ( 0 0.211 0.211 )
+DEAL::7: ( 0.211 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.577 -0.577 -0.577 )
+DEAL::1: ( 0 -0.707 -0.707 )
+DEAL::2: ( -0.707 0 -0.707 )
+DEAL::3: ( 0 0 -1.00 )
+DEAL::4: ( -0.394 -0.394 -0.394 )
+DEAL::5: ( 0 -0.419 -0.419 )
+DEAL::6: ( -0.419 0 -0.419 )
+DEAL::7: ( 0 0 -0.457 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.707 -0.707 )
+DEAL::1: ( 0.577 -0.577 -0.577 )
+DEAL::2: ( 0 0 -1.00 )
+DEAL::3: ( 0.707 0 -0.707 )
+DEAL::4: ( 0 -0.419 -0.419 )
+DEAL::5: ( 0.394 -0.394 -0.394 )
+DEAL::6: ( 0 0 -0.457 )
+DEAL::7: ( 0.419 0 -0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 0 -0.707 )
+DEAL::1: ( 0 0 -1.00 )
+DEAL::2: ( -0.577 0.577 -0.577 )
+DEAL::3: ( 0 0.707 -0.707 )
+DEAL::4: ( -0.419 0 -0.419 )
+DEAL::5: ( 0 0 -0.457 )
+DEAL::6: ( -0.394 0.394 -0.394 )
+DEAL::7: ( 0 0.419 -0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 -1.00 )
+DEAL::1: ( 0.707 0 -0.707 )
+DEAL::2: ( 0 0.707 -0.707 )
+DEAL::3: ( 0.577 0.577 -0.577 )
+DEAL::4: ( 0 0 -0.457 )
+DEAL::5: ( 0.419 0 -0.419 )
+DEAL::6: ( 0 0.419 -0.419 )
+DEAL::7: ( 0.394 0.394 -0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.394 -0.394 -0.394 )
+DEAL::1: ( 0 -0.419 -0.419 )
+DEAL::2: ( -0.419 0 -0.419 )
+DEAL::3: ( 0 0 -0.457 )
+DEAL::4: ( -0.211 -0.211 -0.211 )
+DEAL::5: ( 0 -0.211 -0.211 )
+DEAL::6: ( -0.211 0 -0.211 )
+DEAL::7: ( 0 0 -0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.419 -0.419 )
+DEAL::1: ( 0.394 -0.394 -0.394 )
+DEAL::2: ( 0 0 -0.457 )
+DEAL::3: ( 0.419 0 -0.419 )
+DEAL::4: ( 0 -0.211 -0.211 )
+DEAL::5: ( 0.211 -0.211 -0.211 )
+DEAL::6: ( 0 0 -0.211 )
+DEAL::7: ( 0.211 0 -0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 0 -0.419 )
+DEAL::1: ( 0 0 -0.457 )
+DEAL::2: ( -0.394 0.394 -0.394 )
+DEAL::3: ( 0 0.419 -0.419 )
+DEAL::4: ( -0.211 0 -0.211 )
+DEAL::5: ( 0 0 -0.211 )
+DEAL::6: ( -0.211 0.211 -0.211 )
+DEAL::7: ( 0 0.211 -0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 -0.457 )
+DEAL::1: ( 0.419 0 -0.419 )
+DEAL::2: ( 0 0.419 -0.419 )
+DEAL::3: ( 0.394 0.394 -0.394 )
+DEAL::4: ( 0 0 -0.211 )
+DEAL::5: ( 0.211 0 -0.211 )
+DEAL::6: ( 0 0.211 -0.211 )
+DEAL::7: ( 0.211 0.211 -0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.577 -0.577 -0.577 )
+DEAL::1: ( 0.707 0 -0.707 )
+DEAL::2: ( 0.394 -0.394 -0.394 )
+DEAL::3: ( 0.419 0 -0.419 )
+DEAL::4: ( 0.707 -0.707 0 )
+DEAL::5: ( 1.00 0 0 )
+DEAL::6: ( 0.419 -0.419 0 )
+DEAL::7: ( 0.457 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.707 0 -0.707 )
+DEAL::1: ( 0.577 0.577 -0.577 )
+DEAL::2: ( 0.419 0 -0.419 )
+DEAL::3: ( 0.394 0.394 -0.394 )
+DEAL::4: ( 1.00 0 0 )
+DEAL::5: ( 0.707 0.707 0 )
+DEAL::6: ( 0.457 0 0 )
+DEAL::7: ( 0.419 0.419 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.394 -0.394 -0.394 )
+DEAL::1: ( 0.419 0 -0.419 )
+DEAL::2: ( 0.211 -0.211 -0.211 )
+DEAL::3: ( 0.211 0 -0.211 )
+DEAL::4: ( 0.419 -0.419 0 )
+DEAL::5: ( 0.457 0 0 )
+DEAL::6: ( 0.211 -0.211 0 )
+DEAL::7: ( 0.211 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.419 0 -0.419 )
+DEAL::1: ( 0.394 0.394 -0.394 )
+DEAL::2: ( 0.211 0 -0.211 )
+DEAL::3: ( 0.211 0.211 -0.211 )
+DEAL::4: ( 0.457 0 0 )
+DEAL::5: ( 0.419 0.419 0 )
+DEAL::6: ( 0.211 0 0 )
+DEAL::7: ( 0.211 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.707 -0.707 0 )
+DEAL::1: ( 1.00 0 0 )
+DEAL::2: ( 0.419 -0.419 0 )
+DEAL::3: ( 0.457 0 0 )
+DEAL::4: ( 0.577 -0.577 0.577 )
+DEAL::5: ( 0.707 0 0.707 )
+DEAL::6: ( 0.394 -0.394 0.394 )
+DEAL::7: ( 0.419 0 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 1.00 0 0 )
+DEAL::1: ( 0.707 0.707 0 )
+DEAL::2: ( 0.457 0 0 )
+DEAL::3: ( 0.419 0.419 0 )
+DEAL::4: ( 0.707 0 0.707 )
+DEAL::5: ( 0.577 0.577 0.577 )
+DEAL::6: ( 0.419 0 0.419 )
+DEAL::7: ( 0.394 0.394 0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.419 -0.419 0 )
+DEAL::1: ( 0.457 0 0 )
+DEAL::2: ( 0.211 -0.211 0 )
+DEAL::3: ( 0.211 0 0 )
+DEAL::4: ( 0.394 -0.394 0.394 )
+DEAL::5: ( 0.419 0 0.419 )
+DEAL::6: ( 0.211 -0.211 0.211 )
+DEAL::7: ( 0.211 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.457 0 0 )
+DEAL::1: ( 0.419 0.419 0 )
+DEAL::2: ( 0.211 0 0 )
+DEAL::3: ( 0.211 0.211 0 )
+DEAL::4: ( 0.419 0 0.419 )
+DEAL::5: ( 0.394 0.394 0.394 )
+DEAL::6: ( 0.211 0 0.211 )
+DEAL::7: ( 0.211 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.577 -0.577 0.577 )
+DEAL::1: ( 0 -0.707 0.707 )
+DEAL::2: ( -0.394 -0.394 0.394 )
+DEAL::3: ( 0 -0.419 0.419 )
+DEAL::4: ( -0.707 0 0.707 )
+DEAL::5: ( 0 0 1.00 )
+DEAL::6: ( -0.419 0 0.419 )
+DEAL::7: ( 0 0 0.457 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.707 0.707 )
+DEAL::1: ( 0.577 -0.577 0.577 )
+DEAL::2: ( 0 -0.419 0.419 )
+DEAL::3: ( 0.394 -0.394 0.394 )
+DEAL::4: ( 0 0 1.00 )
+DEAL::5: ( 0.707 0 0.707 )
+DEAL::6: ( 0 0 0.457 )
+DEAL::7: ( 0.419 0 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.394 -0.394 0.394 )
+DEAL::1: ( 0 -0.419 0.419 )
+DEAL::2: ( -0.211 -0.211 0.211 )
+DEAL::3: ( 0 -0.211 0.211 )
+DEAL::4: ( -0.419 0 0.419 )
+DEAL::5: ( 0 0 0.457 )
+DEAL::6: ( -0.211 0 0.211 )
+DEAL::7: ( 0 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.419 0.419 )
+DEAL::1: ( 0.394 -0.394 0.394 )
+DEAL::2: ( 0 -0.211 0.211 )
+DEAL::3: ( 0.211 -0.211 0.211 )
+DEAL::4: ( 0 0 0.457 )
+DEAL::5: ( 0.419 0 0.419 )
+DEAL::6: ( 0 0 0.211 )
+DEAL::7: ( 0.211 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 0 0.707 )
+DEAL::1: ( 0 0 1.00 )
+DEAL::2: ( -0.419 0 0.419 )
+DEAL::3: ( 0 0 0.457 )
+DEAL::4: ( -0.577 0.577 0.577 )
+DEAL::5: ( 0 0.707 0.707 )
+DEAL::6: ( -0.394 0.394 0.394 )
+DEAL::7: ( 0 0.419 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 1.00 )
+DEAL::1: ( 0.707 0 0.707 )
+DEAL::2: ( 0 0 0.457 )
+DEAL::3: ( 0.419 0 0.419 )
+DEAL::4: ( 0 0.707 0.707 )
+DEAL::5: ( 0.577 0.577 0.577 )
+DEAL::6: ( 0 0.419 0.419 )
+DEAL::7: ( 0.394 0.394 0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 0 0.419 )
+DEAL::1: ( 0 0 0.457 )
+DEAL::2: ( -0.211 0 0.211 )
+DEAL::3: ( 0 0 0.211 )
+DEAL::4: ( -0.394 0.394 0.394 )
+DEAL::5: ( 0 0.419 0.419 )
+DEAL::6: ( -0.211 0.211 0.211 )
+DEAL::7: ( 0 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0 0.457 )
+DEAL::1: ( 0.419 0 0.419 )
+DEAL::2: ( 0 0 0.211 )
+DEAL::3: ( 0.211 0 0.211 )
+DEAL::4: ( 0 0.419 0.419 )
+DEAL::5: ( 0.394 0.394 0.394 )
+DEAL::6: ( 0 0.211 0.211 )
+DEAL::7: ( 0.211 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.577 -0.577 -0.577 )
+DEAL::1: ( -0.394 -0.394 -0.394 )
+DEAL::2: ( -0.707 0 -0.707 )
+DEAL::3: ( -0.419 0 -0.419 )
+DEAL::4: ( -0.707 -0.707 0 )
+DEAL::5: ( -0.419 -0.419 0 )
+DEAL::6: ( -1.00 0 0 )
+DEAL::7: ( -0.457 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.394 -0.394 -0.394 )
+DEAL::1: ( -0.211 -0.211 -0.211 )
+DEAL::2: ( -0.419 0 -0.419 )
+DEAL::3: ( -0.211 0 -0.211 )
+DEAL::4: ( -0.419 -0.419 0 )
+DEAL::5: ( -0.211 -0.211 0 )
+DEAL::6: ( -0.457 0 0 )
+DEAL::7: ( -0.211 0 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 0 -0.707 )
+DEAL::1: ( -0.419 0 -0.419 )
+DEAL::2: ( -0.577 0.577 -0.577 )
+DEAL::3: ( -0.394 0.394 -0.394 )
+DEAL::4: ( -1.00 0 0 )
+DEAL::5: ( -0.457 0 0 )
+DEAL::6: ( -0.707 0.707 0 )
+DEAL::7: ( -0.419 0.419 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 0 -0.419 )
+DEAL::1: ( -0.211 0 -0.211 )
+DEAL::2: ( -0.394 0.394 -0.394 )
+DEAL::3: ( -0.211 0.211 -0.211 )
+DEAL::4: ( -0.457 0 0 )
+DEAL::5: ( -0.211 0 0 )
+DEAL::6: ( -0.419 0.419 0 )
+DEAL::7: ( -0.211 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 -0.707 0 )
+DEAL::1: ( -0.419 -0.419 0 )
+DEAL::2: ( -1.00 0 0 )
+DEAL::3: ( -0.457 0 0 )
+DEAL::4: ( -0.577 -0.577 0.577 )
+DEAL::5: ( -0.394 -0.394 0.394 )
+DEAL::6: ( -0.707 0 0.707 )
+DEAL::7: ( -0.419 0 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 -0.419 0 )
+DEAL::1: ( -0.211 -0.211 0 )
+DEAL::2: ( -0.457 0 0 )
+DEAL::3: ( -0.211 0 0 )
+DEAL::4: ( -0.394 -0.394 0.394 )
+DEAL::5: ( -0.211 -0.211 0.211 )
+DEAL::6: ( -0.419 0 0.419 )
+DEAL::7: ( -0.211 0 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -1.00 0 0 )
+DEAL::1: ( -0.457 0 0 )
+DEAL::2: ( -0.707 0.707 0 )
+DEAL::3: ( -0.419 0.419 0 )
+DEAL::4: ( -0.707 0 0.707 )
+DEAL::5: ( -0.419 0 0.419 )
+DEAL::6: ( -0.577 0.577 0.577 )
+DEAL::7: ( -0.394 0.394 0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.457 0 0 )
+DEAL::1: ( -0.211 0 0 )
+DEAL::2: ( -0.419 0.419 0 )
+DEAL::3: ( -0.211 0.211 0 )
+DEAL::4: ( -0.419 0 0.419 )
+DEAL::5: ( -0.211 0 0.211 )
+DEAL::6: ( -0.394 0.394 0.394 )
+DEAL::7: ( -0.211 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.577 -0.577 -0.577 )
+DEAL::1: ( 0 -0.707 -0.707 )
+DEAL::2: ( -0.394 -0.394 -0.394 )
+DEAL::3: ( 0 -0.419 -0.419 )
+DEAL::4: ( -0.707 -0.707 0 )
+DEAL::5: ( 0 -1.00 0 )
+DEAL::6: ( -0.419 -0.419 0 )
+DEAL::7: ( 0 -0.457 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.707 -0.707 )
+DEAL::1: ( 0.577 -0.577 -0.577 )
+DEAL::2: ( 0 -0.419 -0.419 )
+DEAL::3: ( 0.394 -0.394 -0.394 )
+DEAL::4: ( 0 -1.00 0 )
+DEAL::5: ( 0.707 -0.707 0 )
+DEAL::6: ( 0 -0.457 0 )
+DEAL::7: ( 0.419 -0.419 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.394 -0.394 -0.394 )
+DEAL::1: ( 0 -0.419 -0.419 )
+DEAL::2: ( -0.211 -0.211 -0.211 )
+DEAL::3: ( 0 -0.211 -0.211 )
+DEAL::4: ( -0.419 -0.419 0 )
+DEAL::5: ( 0 -0.457 0 )
+DEAL::6: ( -0.211 -0.211 0 )
+DEAL::7: ( 0 -0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.419 -0.419 )
+DEAL::1: ( 0.394 -0.394 -0.394 )
+DEAL::2: ( 0 -0.211 -0.211 )
+DEAL::3: ( 0.211 -0.211 -0.211 )
+DEAL::4: ( 0 -0.457 0 )
+DEAL::5: ( 0.419 -0.419 0 )
+DEAL::6: ( 0 -0.211 0 )
+DEAL::7: ( 0.211 -0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 -0.707 0 )
+DEAL::1: ( 0 -1.00 0 )
+DEAL::2: ( -0.419 -0.419 0 )
+DEAL::3: ( 0 -0.457 0 )
+DEAL::4: ( -0.577 -0.577 0.577 )
+DEAL::5: ( 0 -0.707 0.707 )
+DEAL::6: ( -0.394 -0.394 0.394 )
+DEAL::7: ( 0 -0.419 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -1.00 0 )
+DEAL::1: ( 0.707 -0.707 0 )
+DEAL::2: ( 0 -0.457 0 )
+DEAL::3: ( 0.419 -0.419 0 )
+DEAL::4: ( 0 -0.707 0.707 )
+DEAL::5: ( 0.577 -0.577 0.577 )
+DEAL::6: ( 0 -0.419 0.419 )
+DEAL::7: ( 0.394 -0.394 0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 -0.419 0 )
+DEAL::1: ( 0 -0.457 0 )
+DEAL::2: ( -0.211 -0.211 0 )
+DEAL::3: ( 0 -0.211 0 )
+DEAL::4: ( -0.394 -0.394 0.394 )
+DEAL::5: ( 0 -0.419 0.419 )
+DEAL::6: ( -0.211 -0.211 0.211 )
+DEAL::7: ( 0 -0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 -0.457 0 )
+DEAL::1: ( 0.419 -0.419 0 )
+DEAL::2: ( 0 -0.211 0 )
+DEAL::3: ( 0.211 -0.211 0 )
+DEAL::4: ( 0 -0.419 0.419 )
+DEAL::5: ( 0.394 -0.394 0.394 )
+DEAL::6: ( 0 -0.211 0.211 )
+DEAL::7: ( 0.211 -0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.577 0.577 -0.577 )
+DEAL::1: ( -0.394 0.394 -0.394 )
+DEAL::2: ( 0 0.707 -0.707 )
+DEAL::3: ( 0 0.419 -0.419 )
+DEAL::4: ( -0.707 0.707 0 )
+DEAL::5: ( -0.419 0.419 0 )
+DEAL::6: ( 0 1.00 0 )
+DEAL::7: ( 0 0.457 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.394 0.394 -0.394 )
+DEAL::1: ( -0.211 0.211 -0.211 )
+DEAL::2: ( 0 0.419 -0.419 )
+DEAL::3: ( 0 0.211 -0.211 )
+DEAL::4: ( -0.419 0.419 0 )
+DEAL::5: ( -0.211 0.211 0 )
+DEAL::6: ( 0 0.457 0 )
+DEAL::7: ( 0 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0.707 -0.707 )
+DEAL::1: ( 0 0.419 -0.419 )
+DEAL::2: ( 0.577 0.577 -0.577 )
+DEAL::3: ( 0.394 0.394 -0.394 )
+DEAL::4: ( 0 1.00 0 )
+DEAL::5: ( 0 0.457 0 )
+DEAL::6: ( 0.707 0.707 0 )
+DEAL::7: ( 0.419 0.419 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0.419 -0.419 )
+DEAL::1: ( 0 0.211 -0.211 )
+DEAL::2: ( 0.394 0.394 -0.394 )
+DEAL::3: ( 0.211 0.211 -0.211 )
+DEAL::4: ( 0 0.457 0 )
+DEAL::5: ( 0 0.211 0 )
+DEAL::6: ( 0.419 0.419 0 )
+DEAL::7: ( 0.211 0.211 0 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.707 0.707 0 )
+DEAL::1: ( -0.419 0.419 0 )
+DEAL::2: ( 0 1.00 0 )
+DEAL::3: ( 0 0.457 0 )
+DEAL::4: ( -0.577 0.577 0.577 )
+DEAL::5: ( -0.394 0.394 0.394 )
+DEAL::6: ( 0 0.707 0.707 )
+DEAL::7: ( 0 0.419 0.419 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.419 0.419 0 )
+DEAL::1: ( -0.211 0.211 0 )
+DEAL::2: ( 0 0.457 0 )
+DEAL::3: ( 0 0.211 0 )
+DEAL::4: ( -0.394 0.394 0.394 )
+DEAL::5: ( -0.211 0.211 0.211 )
+DEAL::6: ( 0 0.419 0.419 )
+DEAL::7: ( 0 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 1.00 0 )
+DEAL::1: ( 0 0.457 0 )
+DEAL::2: ( 0.707 0.707 0 )
+DEAL::3: ( 0.419 0.419 0 )
+DEAL::4: ( 0 0.707 0.707 )
+DEAL::5: ( 0 0.419 0.419 )
+DEAL::6: ( 0.577 0.577 0.577 )
+DEAL::7: ( 0.394 0.394 0.394 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0 0.457 0 )
+DEAL::1: ( 0 0.211 0 )
+DEAL::2: ( 0.419 0.419 0 )
+DEAL::3: ( 0.211 0.211 0 )
+DEAL::4: ( 0 0.419 0.419 )
+DEAL::5: ( 0 0.211 0.211 )
+DEAL::6: ( 0.394 0.394 0.394 )
+DEAL::7: ( 0.211 0.211 0.211 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::done...
index 7d5cf4774b073997bb1126e4bfae461853f73f05..4a96a0aa3756df7b9d669ec7fde7db768e134554 100644 (file)
@@ -163,6 +163,9 @@ plot_face_shape_functions(
                               if (uflags & update_hessians)
                                 AssertThrow((fe.shape_hessian(i,k) == fe.shape_hessian_component(i,k,c)),
                                             ExcInternalError());
+                              if (uflags & update_3rd_derivatives)
+                                AssertThrow((fe.shape_3rd_derivative(i,k) == fe.shape_3rd_derivative_component(i,k,c)),
+                                            ExcInternalError());
                             }
                           else
                             {
@@ -175,6 +178,9 @@ plot_face_shape_functions(
                               if (uflags & update_hessians)
                                 AssertThrow ((fe.shape_hessian_component(i,k,c) == Tensor<2,dim>()),
                                              ExcInternalError());
+                              if (uflags & update_3rd_derivatives)
+                                AssertThrow ((fe.shape_3rd_derivative_component(i,k,c) == Tensor<3,dim>()),
+                                             ExcInternalError());
                             }
                         }
                     }
@@ -227,6 +233,12 @@ plot_face_shape_functions(
                                                           s2 = sub.shape_hessian_component(i,k,c);
                                       Assert (s1 == s2, ExcInternalError());
                                     }
+                                  if (uflags & update_3rd_derivatives)
+                                    {
+                                      const Tensor<3,dim> t1 = sub.shape_3rd_derivative(i,k),
+                                                          t2 = sub.shape_3rd_derivative_component(i,k,c);
+                                      Assert (t1 == t2, ExcInternalError());
+                                    }
                                 }
                               else
                                 {
@@ -239,6 +251,9 @@ plot_face_shape_functions(
                                   if (uflags & update_hessians)
                                     Assert ((sub.shape_hessian_component(i,k,c) == Tensor<2,dim>()),
                                             ExcInternalError());
+                                  if (uflags & update_3rd_derivatives)
+                                    Assert ((sub.shape_3rd_derivative_component(i,k,c) == Tensor<3,dim>()),
+                                            ExcInternalError());
                                 }
                             };
                         }

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.