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
// 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)
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;
}
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.
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;
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);
}
}
}
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);
}
}
}
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
* 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.
*/
* 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
*/
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|";
*/
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.
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
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.
*/
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.
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
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.
* <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.
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
//@{
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<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 <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
//@{
* 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);
* 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);
* 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,
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
+ 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
}
}
+ 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
{
+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> &
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);
}
}
}
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);
}
}
}
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);
}
}
}
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);
}
}
}
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())
+ }
}
}
Assert(false, ExcNotImplemented());
}
}
+
+ // third derivatives are not implemented
+ if (flags & update_3rd_derivatives)
+ {
+ Assert(false, ExcNotImplemented())
+ }
}
}
Assert(false, ExcNotImplemented());
}
}
+
+ // third derivatives are not implemented
+ if (flags & update_3rd_derivatives)
+ {
+ Assert(false, ExcNotImplemented())
+ }
}
}
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];
+
}
}
}
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;
+ 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
}
+ 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>
= ("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.");
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));
}
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));
}
}
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
+
--- /dev/null
+
+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...
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
{
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());
}
}
}
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
{
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());
}
};
}
projects / dealii.git / commitdiff