From 8215f0956f845431c296f7ca7e0806e7bbe298c7 Mon Sep 17 00:00:00 2001 From: Maien Hamed Date: Sun, 6 Sep 2015 14:19:33 +0200 Subject: [PATCH] Implemented computation of shape function 3rd derivatives in real cell, except in FE_PolyTensor --- include/deal.II/fe/fe_poly.h | 62 +- include/deal.II/fe/fe_poly.templates.h | 132 +-- include/deal.II/fe/fe_update_flags.h | 47 +- include/deal.II/fe/fe_values.h | 389 ++++++++- source/fe/fe_poly.cc | 108 ++- source/fe/fe_poly_tensor.cc | 23 + source/fe/fe_system.cc | 8 +- source/fe/fe_values.cc | 174 +++- .../fe_q_3rd_derivative_divergence_theorem.cc | 176 ++++ ...q_3rd_derivative_divergence_theorem.output | 760 ++++++++++++++++++ tests/fe/shapes.h | 15 + 11 files changed, 1746 insertions(+), 148 deletions(-) create mode 100644 tests/fe/fe_q_3rd_derivative_divergence_theorem.cc create mode 100644 tests/fe/fe_q_3rd_derivative_divergence_theorem.output diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h index 0e2916fd84..0b03e783df 100644 --- a/include/deal.II/fe/fe_poly.h +++ b/include/deal.II/fe/fe_poly.h @@ -258,16 +258,18 @@ protected: std::vector > (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 > (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 > (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; idofs_per_cell; ++k) @@ -294,6 +299,10 @@ protected: if (flags & update_hessians) for (unsigned int k=0; kdofs_per_cell; ++k) data->shape_hessians[k][i] = grad_grads[k]; + + if (flags & update_3rd_derivatives) + for (unsigned int k=0; kdofs_per_cell; ++k) + data->shape_3rd_derivatives[k][i] = third_derivatives[k]; } return data; } @@ -420,29 +429,34 @@ protected: std::vector > > 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 > untransformed_shape_hessians; + std::vector > > 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 > > uncorrected_shape_hessians, - const internal::FEValues::MappingRelatedData &mapping_data, - const internal::FEValues::FiniteElementRelatedData &fevalues_data, - const unsigned int n_q_points, - const unsigned int dof) const; + correct_third_derivatives (internal::FEValues::FiniteElementRelatedData &output_data, + const internal::FEValues::MappingRelatedData &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. diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index 6050446263..bdaa97d2e1 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -202,7 +202,12 @@ FE_Poly::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 &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 if (flags & update_hessians) { - // compute the hessians in the unit cell (accounting for the Jacobian gradiant) - for (unsigned int i=0; i > > - ( 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 if (flags & update_hessians) { - // compute the hessians in the unit cell (accounting for the Jacobian gradiant) - for (unsigned int i=0; i > > - (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 -void +inline void FE_Poly:: -correct_untransformed_hessians (VectorSlice< std::vector > > uncorrected_shape_hessians, - const internal::FEValues::MappingRelatedData &mapping_data, - const internal::FEValues::FiniteElementRelatedData &fevalues_data, - const unsigned int n_q_points, - const unsigned int dof) const +correct_third_derivatives (internal::FEValues::FiniteElementRelatedData &output_data, + const internal::FEValues::MappingRelatedData &mapping_data, + const unsigned int n_q_points, + const unsigned int dof) const { for (unsigned int i=0; i 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 Tensor@<3,dim@>. + */ + 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 fe_function at the @@ -339,6 +360,29 @@ namespace FEValuesViews void get_function_laplacians (const InputVector &fe_function, std::vector::type> &laplacians) const; + /** + * Return the third derivatives of the selected scalar component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit 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 + void get_function_third_derivatives (const InputVector &fe_function, + std::vector::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 dim components of the + * finite element, the third derivative is a Tensor@<4,dim@>. + */ + 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 fe_function at the @@ -765,6 +830,29 @@ namespace FEValuesViews void get_function_laplacians (const InputVector &fe_function, std::vector::type> &laplacians) const; + /** + * Return the third derivatives of the selected scalar component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit 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 + void get_function_third_derivatives (const InputVector &fe_function, + std::vector::type> &third_derivatives) const; + private: /** * A reference to the FEValuesBase object we operate on. @@ -1509,14 +1597,14 @@ public: * point_noth 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 function_noth shape function at the + * point_noth 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 third_derivatives[q] will contain the third derivatives + * of the field described by fe_function at the $q$th quadrature point. + * third_derivatives[q][i][j][k] 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<T>, BlockVector<T>, 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 + void + get_function_third_derivatives (const InputVector &fe_function, + std::vector > &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 third_derivatives[q] 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 third_derivatives[q] equals + * the number of components of the finite element, i.e. + * third_derivatives[q][c] returns the third derivative of the + * $c$th vector component at the $q$th quadrature point. + * Consequently, third_derivatives[q][c][i][j][k] 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 + void + get_function_third_derivatives (const InputVector &fe_function, + std::vector > > &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 + void get_function_third_derivatives ( + const InputVector &fe_function, + const VectorSlice > &indices, + std::vector > &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 + void get_function_third_derivatives ( + const InputVector &fe_function, + const VectorSlice > &indices, + VectorSlice > > > 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 - * get_function_value/gradients/hessians/laplacians 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. + * get_function_value/gradients/hessians/laplacians/third_derivatives + * 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::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 - * get_function_value/gradients/hessians 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. + * get_function_value/gradients/hessians/third_derivatives + * 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::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 - * get_function_value/gradients/hessians 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. + * get_function_value/gradients/hessians/third_derivatives + * 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::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 + inline + typename Scalar::third_derivative_type + Scalar::third_derivative (const unsigned int shape_function, + const unsigned int q_point) const + { + typedef FEValuesBase 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 inline typename Vector::value_type @@ -3335,6 +3593,44 @@ namespace FEValuesViews } } + template + inline + typename Vector::third_derivative_type + Vector::third_derivative (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in + // the case above + typedef FEValuesBase 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::shape_hessian_component (const unsigned int i, +template +inline +const Tensor<3,spacedim> & +FEValuesBase::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 (ifinite_element_output.shape_3rd_derivatives.size(), + ExcIndexRange (i, 0, this->finite_element_output.shape_3rd_derivatives.size())); + Assert (jfinite_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 +inline +Tensor<3,spacedim> +FEValuesBase::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 inline const FiniteElement & diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc index 6013677501..7981da4946 100644 --- a/source/fe/fe_poly.cc +++ b/source/fe/fe_poly.cc @@ -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 &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 &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 &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 &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 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 Assert(false, ExcNotImplemented()); } } + + // third derivatives are not implemented + if (flags & update_3rd_derivatives) + { + Assert(false, ExcNotImplemented()) + } } } diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 2be7e80d48..afd837eca0 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1233,6 +1233,11 @@ compute_fill_one_base (const Mapping &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 &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; diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index faa37b8b6e..f5341ebee6 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -1356,6 +1356,30 @@ namespace FEValuesViews + template + template + void + Scalar:: + get_function_third_derivatives (const InputVector &fe_function, + std::vector::type> &third_derivatives) const + { + typedef FEValuesBase 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 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 template void @@ -1527,6 +1551,29 @@ namespace FEValuesViews } + template + template + void + Vector:: + get_function_third_derivatives (const InputVector &fe_function, + std::vector::type> &third_derivatives) const + { + typedef FEValuesBase 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 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 template @@ -2023,7 +2070,7 @@ FEValuesBase::TriaCellIterator::message_string = ("You have previously called the FEValues::reinit function with a\n" "cell iterator of type Triangulation::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::cell_iterator."); @@ -2178,6 +2225,10 @@ namespace internal if (flags & update_hessians) this->shape_hessians.resize (n_nonzero_shape_components, std::vector > (n_quadrature_points)); + + if (flags & update_3rd_derivatives) + this->shape_3rd_derivatives.resize (n_nonzero_shape_components, + std::vector > (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::get_function_laplacians ( +template +template +void +FEValuesBase:: +get_function_third_derivatives (const InputVector &fe_function, + std::vector > &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 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 +template +void FEValuesBase::get_function_third_derivatives ( + const InputVector &fe_function, + const VectorSlice > &indices, + std::vector > &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; ifinite_element_output.shape_3rd_derivatives, + third_derivatives); + } + else + { + Vector dof_values(dofs_per_cell); + for (unsigned int i=0; ifinite_element_output.shape_3rd_derivatives, + third_derivatives); + } +} + + + + +template +template +void +FEValuesBase:: +get_function_third_derivatives (const InputVector &fe_function, + std::vector > > &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 dof_values (dofs_per_cell); + present_cell->get_interpolated_dof_values(fe_function, dof_values); + VectorSlice > > > 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 +template +void FEValuesBase::get_function_third_derivatives ( + const InputVector &fe_function, + const VectorSlice > &indices, + VectorSlice > > > 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; ifinite_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 dof_values(indices.size()); + for (unsigned int i=0; ifinite_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 const typename Triangulation::cell_iterator FEValuesBase::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 index 0000000000..c744720bcb --- /dev/null +++ b/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +template +Tensor<1,dim> ones () +{ + Tensor<1,dim> result; + for (unsigned int i=0; i +void test (const Triangulation &tr, + const FiniteElement &fe, + const double tolerance) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + std::stringstream ss; + + deallog << "FE=" << fe.get_name() << std::endl; + + const QGauss quadrature(6); + FEValues fe_values (fe, quadrature, update_3rd_derivatives + | update_quadrature_points + | update_JxW_values); + + const QGauss face_quadrature(6); + FEFaceValues fe_face_values (fe, face_quadrature, + update_hessians + | update_quadrature_points + | update_normal_vectors + | update_JxW_values); + + for (typename DoFHandler::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::vertices_per_cell; ++i) + { + deallog << i << ": ( "; + for (unsigned int d=0; dvertex(i)[d] << " "; + deallog << ")" << std::endl; + } + + bool cell_ok = true; + + for (unsigned int c=0; c +void test_hyper_ball(const double tolerance) +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + tr.refine_global(1); + + FE_Q 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 index 0000000000..2a7c032e52 --- /dev/null +++ b/tests/fe/fe_q_3rd_derivative_divergence_theorem.output @@ -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... diff --git a/tests/fe/shapes.h b/tests/fe/shapes.h index 7d5cf4774b..4a96a0aa37 100644 --- a/tests/fe/shapes.h +++ b/tests/fe/shapes.h @@ -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()); } }; } -- 2.39.5