From: kronbichler Date: Mon, 15 Dec 2008 23:16:45 +0000 (+0000) Subject: Add get_function_laplacians() methods in FEValues class, improve performance of get_f... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=71c4b07466fa66b8693f4f7b14ed57873727fac7;p=dealii-svn.git Add get_function_laplacians() methods in FEValues class, improve performance of get_function_xxx in views, use new functions in step-31 and step-32. git-svn-id: https://svn.dealii.org/trunk@17954 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index a3b8331285..61b36bce72 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -281,7 +281,7 @@ namespace FEValuesViews */ template void get_function_gradients (const InputVector& fe_function, - std::vector& values) const; + std::vector& gradients) const; /** * Return the Hessians of the selected @@ -302,7 +302,33 @@ namespace FEValuesViews */ template void get_function_hessians (const InputVector& fe_function, - std::vector& values) const; + std::vector& hessians) const; + + /** + * Return the Laplacians 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, at the quadrature + * points. The Laplacians are + * the trace of the Hessians. + * + * This function is the + * equivalent of the + * FEValuesBase::get_function_laplacians + * function but it only works + * on the selected scalar + * component. + */ + template + void get_function_laplacians (const InputVector& fe_function, + std::vector& laplacians) const; private: @@ -431,7 +457,7 @@ namespace FEValuesViews * finite element, the Hessian is a * Tensor@<3,dim@>. */ - typedef Tensor<3,dim> hessian_type; + typedef Tensor<3,spacedim> hessian_type; /** * Default constructor. Creates an @@ -520,7 +546,7 @@ namespace FEValuesViews divergence_type divergence (const unsigned int shape_function, const unsigned int q_point) const; - + /** * Return the Hessian (the tensor of * rank 2 of all second derivatives) of @@ -573,7 +599,7 @@ namespace FEValuesViews */ template void get_function_gradients (const InputVector& fe_function, - std::vector& values) const; + std::vector& gradients) const; /** * Return the symmetrized gradients of @@ -597,7 +623,7 @@ namespace FEValuesViews */ template void get_function_symmetric_gradients (const InputVector& fe_function, - std::vector& values) const; + std::vector& symmetric_gradients) const; /** * Return the divergence of the selected @@ -620,7 +646,7 @@ namespace FEValuesViews */ template void get_function_divergences (const InputVector& fe_function, - std::vector& values) const; + std::vector& divergences) const; /** * Return the Hessians of the selected @@ -641,7 +667,33 @@ namespace FEValuesViews */ template void get_function_hessians (const InputVector& fe_function, - std::vector& values) const; + std::vector& hessians) const; + + /** + * Return the Laplacians of the + * selected vector components + * 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, at the quadrature + * points. The Laplacians are + * the trace of the Hessians. + * + * This function is the + * equivalent of the + * FEValuesBase::get_function_laplacians + * function but it only works + * on the selected vector + * components. + */ + template + void get_function_laplacians (const InputVector& fe_function, + std::vector& laplacians) const; private: /** @@ -972,7 +1024,7 @@ class FEValuesData */ std::vector > normal_vectors; - /** + /** * List of outward vectors normal to the cell * surface (line) at the quadrature points * for the codimension 1 case, @@ -1848,18 +1900,19 @@ class FEValuesBase : protected FEValuesData, * BlockVector<T>, or one * of the sequential PETSc or * Trilinos vector wrapper - * classes..It represents a + * classes. It represents a * global vector of DoF values * associated with the DofHandler * object with which this * FEValues object was last * initialized. * - * The output are the second derivatives - * of the function represented by - * these DoF values, as computed - * in real space (as opposed to - * on the unit cell). + * The output are the second + * derivatives of the function + * represented by these DoF + * values, as computed in real + * space (as opposed to on the + * unit cell). */ template void @@ -1919,7 +1972,7 @@ class FEValuesBase : protected FEValuesData, void get_function_hessians ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector >& result) const; + std::vector >& hessians) const; /** * Access to the second @@ -1932,7 +1985,7 @@ class FEValuesBase : protected FEValuesData, void get_function_hessians ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector > >& result, + std::vector > >& hessians, bool quadrature_points_fastest = false) const; /** @@ -1953,6 +2006,151 @@ class FEValuesBase : protected FEValuesData, bool = false) const; + /** + * Compute the (scalar) Laplacian + * of the finite element function + * characterized by @p + * fe_function restricted to @p + * cell at the quadrature + * points. The Laplacian output + * vector is equivalent to + * getting + * trace(hessians), + * where hessian would + * be the output of the + * get_function_hessians() + * function. + * + * The function assumes that the + * @p laplacians object + * already has the correct size. + * + * This function may only be used if the + * finite element in use is a scalar one, + * i.e. has only one vector component. If + * it is a vector-valued one, then use + * the other + * get_function_laplacians() + * function. + * + * 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. + * + * The output are the traces of + * the second derivatives + * (i.e. Laplacians) of the + * function represented by these + * DoF values, as computed in + * real space (as opposed to on + * the unit cell). + */ + template + void + get_function_laplacians (const InputVector& fe_function, + std::vector& laplacians) const; + + + /** + * Compute the (scalar) Laplacian + * of the finite element function + * characterized by @p + * fe_function restricted to @p + * cell at the quadrature + * points. The Laplacian output + * vector is equivalent to + * getting + * trace(hessians), with + * hessian corresponding + * to the output of the + * get_function_hessians() + * function. + * + * The function assumes that the + * @p laplacians object + * already has the correct size. + * + * This function does the same as + * the other one with the same + * name, but applies to + * vector-valued finite elements. + * + * 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. + * + * The output are the traces of + * the second derivatives (i.e. + * Laplacians) of the function + * represented by these DoF + * values, as computed in real + * space (as opposed to on the + * unit cell). + */ + template + void + get_function_laplacians (const InputVector &fe_function, + std::vector > &laplacians) const; + + /** + * Access to the second + * derivatives of a function with + * more flexibility. see + * get_function_values() with + * corresponding arguments. + */ + template + void get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector& laplacians) const; + + /** + * Access to the second + * derivatives of a function with + * more flexibility. see + * get_function_values() with + * corresponding arguments. + */ + template + void get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector >& laplacians) const; + + /** + * Access to the second + * derivatives of a function with + * more flexibility. see + * get_function_values() with + * corresponding arguments. + */ + template + void get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector >& laplacians, + bool quadrature_points_fastest = false) const; + + //@} /** @@ -2573,7 +2771,7 @@ class FEValues : public FEValuesBase * about degrees of * freedom. These functions are, * above all, the - * get_function_value/gradients/hessians + * get_function_value/gradients/hessians/laplacians * functions. If you want to call * these functions, you have to * call the @p reinit variants @@ -2771,7 +2969,7 @@ class FEFaceValuesBase : public FEValuesBase * function was called. */ unsigned int present_face_index; - + /** * Store a copy of the quadrature * formula here. @@ -3269,122 +3467,6 @@ namespace FEValuesViews - template - template - inline - void - Scalar:: - get_function_values (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_values, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), value_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - if (shape_function_data[shape_function].is_nonzero_shape_function_component) - for (unsigned int q_point=0; q_point - template - inline - void - Scalar:: - get_function_gradients (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_gradients, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), gradient_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - if (shape_function_data[shape_function].is_nonzero_shape_function_component) - for (unsigned int q_point=0; q_point - template - inline - void - Scalar:: - get_function_hessians (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_hessians, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), hessian_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - for (unsigned int q_point=0; q_point inline typename Vector::value_type @@ -3646,270 +3728,6 @@ namespace FEValuesViews return symmetrize(return_value); } } - - - - template - template - inline - void - Vector:: - get_function_values (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_values, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), value_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - { - const int snc = shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the - // selected components - continue; - else if (snc != -1) - for (unsigned int q_point=0; q_point - template - inline - void - Vector:: - get_function_gradients (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_gradients, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), gradient_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - { - const int snc = shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the - // selected components - continue; - else if (snc != -1) - for (unsigned int q_point=0; q_point - template - inline - void - Vector:: - get_function_symmetric_gradients (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_gradients, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), symmetric_gradient_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - { - const int snc = shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the - // selected components - continue; - else if (snc != -1) - for (unsigned int q_point=0; q_point - template - inline - void - Vector:: - get_function_divergences (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_gradients, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), divergence_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - { - const int snc = shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the - // selected components - continue; - else if (snc != -1) - for (unsigned int q_point=0; q_point - template - inline - void - Vector:: - get_function_hessians (const InputVector &fe_function, - std::vector &values) const - { - typedef FEValuesBase FVB; - Assert (fe_values.update_flags & update_hessians, - typename FVB::ExcAccessToUninitializedField()); - Assert (values.size() == fe_values.n_quadrature_points, - ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); - Assert (fe_values.present_cell.get() != 0, - ExcMessage ("FEValues object is not reinit'ed to any cell")); - Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), - ExcDimensionMismatch(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); - - std::fill (values.begin(), values.end(), hessian_type()); - - for (unsigned int shape_function=0; - shape_functiondofs_per_cell; ++shape_function) - { - const int snc = shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the - // selected components - continue; - else if (snc != -1) - for (unsigned int q_point=0; q_point::shape_grad (const unsigned int i, -template +template inline Tensor<1,spacedim> FEValuesBase::shape_grad_component (const unsigned int i, diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 235a58f586..fad84c5413 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -57,7 +57,7 @@ namespace FEValuesViews { template Scalar::Scalar (const FEValuesBase &fe_values, - const unsigned int component) + const unsigned int component) : fe_values (fe_values), component (component), @@ -222,6 +222,525 @@ namespace FEValuesViews Assert (false, ExcInternalError()); return *this; } + + + + template + template + void + Scalar:: + get_function_values (const InputVector &fe_function, + std::vector &values) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_values, + typename FVB::ExcAccessToUninitializedField()); + Assert (values.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (values.begin(), values.end(), value_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + if (shape_function_data[shape_function].is_nonzero_shape_function_component) + { + const double value = dof_values(shape_function); + const double * shape_value_ptr = + &fe_values.shape_values(shape_function_data[shape_function].row_index, 0); + for (unsigned int q_point=0; q_point + template + void + Scalar:: + get_function_gradients (const InputVector &fe_function, + std::vector &gradients) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + Assert (gradients.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (gradients.begin(), gradients.end(), gradient_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + if (shape_function_data[shape_function].is_nonzero_shape_function_component) + { + const double value = dof_values(shape_function); + const Tensor<1,spacedim> * shape_gradient_ptr = + &fe_values.shape_gradients[shape_function_data[shape_function]. + row_index][0]; + for (unsigned int q_point=0; q_point + template + void + Scalar:: + get_function_hessians (const InputVector &fe_function, + std::vector &hessians) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_hessians, + typename FVB::ExcAccessToUninitializedField()); + Assert (hessians.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(hessians.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (hessians.begin(), hessians.end(), hessian_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + if (shape_function_data[shape_function].is_nonzero_shape_function_component) + { + const double value = dof_values(shape_function); + const Tensor<2,spacedim> * shape_hessian_ptr = + &fe_values.shape_hessians[shape_function_data[shape_function]. + row_index][0]; + for (unsigned int q_point=0; q_point + template + void + Scalar:: + get_function_laplacians (const InputVector &fe_function, + std::vector &laplacians) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_hessians, + typename FVB::ExcAccessToUninitializedField()); + Assert (laplacians.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (laplacians.begin(), laplacians.end(), value_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + if (shape_function_data[shape_function].is_nonzero_shape_function_component) + { + const double value = dof_values(shape_function); + const unsigned int row_index = shape_function_data[shape_function].row_index; + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_values (const InputVector &fe_function, + std::vector &values) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_values, + typename FVB::ExcAccessToUninitializedField()); + Assert (values.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (values.begin(), values.end(), value_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const double * shape_value_ptr = &fe_values.shape_values(snc,0); + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_gradients (const InputVector &fe_function, + std::vector &gradients) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + Assert (gradients.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (gradients.begin(), gradients.end(), gradient_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<1,spacedim> * shape_gradient_ptr = + &fe_values.shape_gradients[snc][0]; + for (unsigned int q_point=0; q_point * shape_gradient_ptr = + &fe_values.shape_gradients[shape_function_data[shape_function]. + row_index[d]][0]; + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_symmetric_gradients (const InputVector &fe_function, + std::vector &symmetric_gradients) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + Assert (symmetric_gradients.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(symmetric_gradients.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (symmetric_gradients.begin(), symmetric_gradients.end(), + symmetric_gradient_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<1,spacedim> * shape_gradient_ptr = + &fe_values.shape_gradients[snc][0]; + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_divergences (const InputVector &fe_function, + std::vector &divergences) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + Assert (divergences.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (divergences.begin(), divergences.end(), divergence_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<1,spacedim> * shape_gradient_ptr = + &fe_values.shape_gradients[snc][0]; + for (unsigned int q_point=0; q_point * shape_gradient_ptr = + &fe_values.shape_gradients[shape_function_data[shape_function]. + row_index[d]][0]; + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_hessians (const InputVector &fe_function, + std::vector &hessians) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_hessians, + typename FVB::ExcAccessToUninitializedField()); + Assert (hessians.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(hessians.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (hessians.begin(), hessians.end(), hessian_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<2,spacedim> * shape_hessian_ptr = + &fe_values.shape_hessians[snc][0]; + for (unsigned int q_point=0; q_point * shape_hessian_ptr = + &fe_values.shape_hessians[shape_function_data[shape_function]. + row_index[d]][0]; + for (unsigned int q_point=0; q_point + template + void + Vector:: + get_function_laplacians (const InputVector &fe_function, + std::vector &laplacians) const + { + typedef FEValuesBase FVB; + Assert (fe_values.update_flags & update_hessians, + typename FVB::ExcAccessToUninitializedField()); + Assert (laplacians.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + std::fill (laplacians.begin(), laplacians.end(), value_type()); + + for (unsigned int shape_function=0; + shape_functiondofs_per_cell; ++shape_function) + { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the + // selected components + continue; + + const double value = dof_values(shape_function); + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<2,spacedim> * shape_hessian_ptr = + &fe_values.shape_hessians[snc][0]; + for (unsigned int q_point=0; q_point * shape_hessian_ptr = + &fe_values.shape_hessians[shape_function_data[shape_function]. + row_index[d]][0]; + for (unsigned int q_point=0; q_point::get_function_values ( // deal with scalar finite // elements, so no need to check // for non-primitivity of shape - // functions + // functions. in order to increase + // the speed of this function, we + // directly access the data in the + // shape_values array, and + // increment pointers for accessing + // the data. this saves some lookup + // time and indexing. moreover, the + // order of the loops is such that + // we can access the shape_values + // data stored contiguously (which + // is also advantageous because + // access to dof_values is + // generally more expensive than + // access to the std::vector values + // - so we do the cheaper operation + // in the innermost loop) for (unsigned int shape_func=0; shape_func::get_function_values ( // deal with scalar finite // elements, so no need to check // for non-primitivity of shape - // functions + // functions. in order to increase + // the speed of this function, we + // directly access the data in the + // shape_values array, and + // increment pointers for accessing + // the data. this saves some lookup + // time and indexing. moreover, the + // order of the loops is such that + // we can access the shape_values + // data stored contiguously (which + // is also advantageous because + // access to the global vector + // fe_function is more expensive + // than access to the small + // std::vector values - so we do + // the cheaper operation in the + // innermost loop) for (unsigned int shape_func=0; shape_func::get_function_values ( // shape function is primitive or // not. if it is, then set its only // non-zero component, otherwise - // loop over components + // loop over components. in order + // to increase the speed of this + // function, we directly access the + // data in the shape_values array, + // and increment pointers for + // accessing the data. this saves + // some lookup time and + // indexing. moreover, in order of + // the loops is such that we can + // access the shape_values data + // stored contiguously (which is + // also advantageous because access + // to the global vector fe_function + // is more expensive than access to + // the small std::vector values - + // so we do the cheaper operation + // in the innermost loop) for (unsigned int shape_func=0; shape_funcis_primitive(shape_func)) { + // compared to the scalar + // functions, finding the correct + // index in the shape_value table + // is more involved, since we have + // to find the row in shape_values + // that corresponds to the present + // shape_func. this is done + // manually in the same way as in + // shape_value_component() (that + // function can't be used because + // it doesn't return us a pointer + // to the data). const unsigned int row = fe->is_primitive() ? shape_func : this->shape_function_to_row_table[shape_func]; @@ -1334,6 +1912,8 @@ void FEValuesBase::get_function_values ( for (unsigned int point=0; point::get_function_values ( const double *shape_value_ptr = &this->shape_values(row, 0); const unsigned int comp = fe->system_to_component_index(shape_func).first + mc * n_components; + if (quadrature_points_fastest) for (unsigned int point=0; point::get_function_gradients ( // deal with scalar finite // elements, so no need to check // for non-primitivity of shape - // functions + // functions. in order to increase + // the speed of this function, we + // directly access the data in the + // shape_gradients array, and + // increment pointers for accessing + // the data. this saves some lookup + // time and indexing. moreover, the + // order of the loops is such that + // we can access the + // shape_gradients data stored + // contiguously (which is also + // advantageous because access to + // the vector dof_values is + // gerenally more expensive than + // access to the std::vector + // gradients - so we do the cheaper + // operation in the innermost loop) for (unsigned int shape_func=0; shape_func *shape_gradient_ptr = &this->shape_gradients[shape_func][0]; + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[shape_func][0]; + for (unsigned int point=0; point::get_function_gradients ( // deal with scalar finite // elements, so no need to check // for non-primitivity of shape - // functions + // functions. in order to increase + // the speed of this function, we + // directly access the data in the + // shape_gradients array, and + // increment pointers for accessing + // the data. this saves some lookup + // time and indexing. moreover, the + // order of the loops is such that + // we can access the + // shape_gradients data stored + // contiguously (which is also + // advantageous because access to + // the global vector fe_function is + // more expensive than access to + // the small std::vector gradients + // - so we do the cheaper operation + // in the innermost loop) for (unsigned int shape_func=0; shape_func *shape_gradient_ptr = &this->shape_gradients[shape_func][0]; + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[shape_func][0]; + for (unsigned int point=0; point *shape_hessians_ptr = &this->shape_hessians[shape_func][0]; + const Tensor<2,spacedim> *shape_hessians_ptr + = &this->shape_hessians[shape_func][0]; + for (unsigned int point=0; point::get_function_hessians ( for (unsigned int shape_func=0; shape_func *shape_hessians_ptr = &this->shape_hessians[shape_func][0]; + const Tensor<2,spacedim> *shape_hessians_ptr + = &this->shape_hessians[shape_func][0]; + for (unsigned int point=0; point::get_function_hessians ( +template +template +void FEValuesBase::get_function_laplacians ( + const InputVector &fe_function, + std::vector &laplacians) const +{ + Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField()); + Assert (fe->n_components() == 1, + ExcDimensionMismatch(fe->n_components(), 1)); + Assert (laplacians.size() == n_quadrature_points, + ExcDimensionMismatch(laplacians.size(), n_quadrature_points)); + Assert (present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + // initialize with zero + std::fill_n (laplacians.begin(), n_quadrature_points, 0); + + // add up contributions of trial + // functions. note that here we + // deal with scalar finite + // elements, so no need to check + // for non-primitivity of shape + // functions. note that the + // laplacian is the trace of the + // hessian, so we use a pointer to + // the hessians and get their trace + for (unsigned int shape_func=0; shape_func *shape_hessian_ptr + = &this->shape_hessians[shape_func][0]; + + for (unsigned int point=0; point +template +void FEValuesBase::get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector &laplacians) const +{ + Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField()); + // This function fills a single + // component only + Assert (fe->n_components() == 1, + ExcDimensionMismatch(fe->n_components(), 1)); + // One index for each dof + Assert (indices.size() == dofs_per_cell, + ExcDimensionMismatch(indices.size(), dofs_per_cell)); + // This vector has one entry for + // each quadrature point + Assert (laplacians.size() == n_quadrature_points, + ExcDimensionMismatch(laplacians.size(), n_quadrature_points)); + + // initialize with zero + std::fill_n (laplacians.begin(), n_quadrature_points, 0); + + // add up contributions of trial + // functions. note that here we + // deal with scalar finite + // elements, so no need to check + // for non-primitivity of shape + // functions. note that the + // laplacian is the trace of the + // hessian, so we use a pointer to + // the hessians and get their trace + for (unsigned int shape_func=0; shape_func *shape_hessian_ptr + = &this->shape_hessians[shape_func][0]; + + for (unsigned int point=0; point +template +void FEValuesBase::get_function_laplacians ( + const InputVector& fe_function, + std::vector >& laplacians) const +{ +//TODO: Find out how to do this assertion. + // This vector must correspond to a + // complete discretization +// Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), +// ExcDimensionMismatch(fe_function.size(), +// present_cell->get_dof_handler().n_dofs())); + // One entry per quadrature point + Assert (present_cell.get() != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (laplacians.size() == n_quadrature_points, + ExcDimensionMismatch(laplacians.size(), n_quadrature_points)); + + const unsigned int n_components = fe->n_components(); + // Assert that we can write all + // components into the result + // vectors + for (unsigned i=0;iupdate_flags & update_hessians, ExcAccessToUninitializedField()); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcDimensionMismatch(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); + + // initialize with zero + for (unsigned i=0;iis_primitive(shape_func)) + { + const unsigned int + row = fe->is_primitive() ? + shape_func : this->shape_function_to_row_table[shape_func]; + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first; + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + + for (unsigned int point=0; point +template +void FEValuesBase::get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector >& laplacians) const +{ + // One value per quadrature point + Assert (n_quadrature_points == laplacians.size(), + ExcDimensionMismatch(laplacians.size(), n_quadrature_points)); + + const unsigned int n_components = fe->n_components(); + + // Size of indices must be a + // multiple of dofs_per_cell such + // that an integer number of + // function values is generated in + // each point. + Assert (indices.size() % dofs_per_cell == 0, + ExcNotMultiple(indices.size(), dofs_per_cell)); + + // The number of components of the + // result may be a multiple of the + // number of components of the + // finite element + const unsigned int result_components = indices.size() * n_components / dofs_per_cell; + + for (unsigned i=0;iupdate_flags & update_hessians, ExcAccessToUninitializedField()); + + // initialize with zero + for (unsigned i=0;iis_primitive(shape_func)) + { + const unsigned int + row = fe->is_primitive() ? + shape_func : this->shape_function_to_row_table[shape_func]; + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = c + mc * n_components; + + for (unsigned int point=0; point +template +void FEValuesBase::get_function_laplacians ( + const InputVector& fe_function, + const VectorSlice >& indices, + std::vector >& laplacians, + bool quadrature_points_fastest) const +{ + const unsigned int n_components = fe->n_components(); + + // Size of indices must be a + // multiple of dofs_per_cell such + // that an integer number of + // function values is generated in + // each point. + Assert (indices.size() % dofs_per_cell == 0, + ExcNotMultiple(indices.size(), dofs_per_cell)); + + // The number of components of the + // result may be a multiple of the + // number of components of the + // finite element + const unsigned int result_components = indices.size() * n_components / dofs_per_cell; + + // Check if the value argument is + // initialized to the correct sizes + if (quadrature_points_fastest) + { + Assert (laplacians.size() == result_components, + ExcDimensionMismatch(laplacians.size(), result_components)); + for (unsigned i=0;iupdate_flags & update_hessians, ExcAccessToUninitializedField()); + + // initialize with zero + for (unsigned i=0;iis_primitive(shape_func)) + { + const unsigned int + row = fe->is_primitive() ? + shape_func : this->shape_function_to_row_table[shape_func]; + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + if (quadrature_points_fastest) + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = c + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; point const std::vector > & FEValuesBase::get_cell_normal_vectors () const diff --git a/deal.II/deal.II/source/fe/fe_values.inst.in b/deal.II/deal.II/source/fe/fe_values.inst.in index a85fcd5ae1..5009baf9a5 100644 --- a/deal.II/deal.II/source/fe/fe_values.inst.in +++ b/deal.II/deal.II/source/fe/fe_values.inst.in @@ -11,10 +11,97 @@ // //----------------------------------------------------------------------- +// Instantiations of functions in FEValuesViews (explicitly write +// dealii:: namespace in order to not confuse the compiler with +// FEValuesViews::Vector). +for (VEC : SERIAL_VECTORS) + { + template + void FEValuesViews::Scalar::get_function_values + (const dealii::VEC&, std::vector&) const; + template + void FEValuesViews::Scalar::get_function_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Scalar::get_function_hessians + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Scalar::get_function_laplacians + (const dealii::VEC&, std::vector&) const; + + template + void FEValuesViews::Vector::get_function_values + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector::get_function_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector::get_function_symmetric_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector::get_function_divergences + (const dealii::VEC&, std::vector&) const; + template + void FEValuesViews::Vector::get_function_hessians + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector::get_function_laplacians + (const dealii::VEC&, std::vector >&) const; + + +#if deal_II_dimension != 3 + template + void FEValuesViews::Scalar + ::get_function_values + (const dealii::VEC&, std::vector&) const; + template + void FEValuesViews::Scalar + ::get_function_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Scalar + ::get_function_hessians + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Scalar + ::get_function_laplacians + (const dealii::VEC&, std::vector&) const; + + template + void FEValuesViews::Vector + ::get_function_values + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector + ::get_function_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector + ::get_function_symmetric_gradients + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector + ::get_function_divergences + (const dealii::VEC&, std::vector&) const; + template + void FEValuesViews::Vector + ::get_function_hessians + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Vector + ::get_function_laplacians + (const dealii::VEC&, std::vector >&) const; + +#endif + + } + + + // Instantiations of functions in FEValuesBase for (VEC : SERIAL_VECTORS) - { + { template void FEValuesBase::get_function_values (const VEC&, std::vector&) const; @@ -85,6 +172,44 @@ for (VEC : SERIAL_VECTORS) (const VEC&, const VectorSlice >&, std::vector > > &, bool) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, std::vector&) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector > &) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector > &) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &, bool) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &, bool) const; + #if deal_II_dimension != 3 @@ -158,6 +283,45 @@ for (VEC : SERIAL_VECTORS) (const VEC&, const VectorSlice >&, std::vector > > &, bool) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, std::vector&) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, std::vector&) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector > &) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, std::vector > &) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &) const; + + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &, bool) const; + template + void FEValuesBase::get_function_laplacians + (const VEC&, const VectorSlice >&, + std::vector > &, bool) const; + + #endif } diff --git a/deal.II/doc/doxygen/headers/vector_valued.h b/deal.II/doc/doxygen/headers/vector_valued.h index 08a9800db6..8c7f8b7045 100644 --- a/deal.II/doc/doxygen/headers/vector_valued.h +++ b/deal.II/doc/doxygen/headers/vector_valued.h @@ -868,6 +868,11 @@ scalar_product (const Tensor<2,dim> &u, * FEValuesViews::Vector::get_function_symmetric_gradients and * FEValuesViews::Vector::get_function_divergences. * + * Moreover, there is a shortcut available in case when only the + * Laplacians of the solution (which is the trace of the hessians) is + * needed, usable for both scalar and vector-valued problems as + * FEValuesViews::Scalar::get_function_laplacians and + * FEValuesViews::Vector::get_function_laplacians. * * @ingroup feall feaccess */ diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 70a46561be..4020dd77eb 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -598,6 +598,16 @@ inconvenience this causes.

deal.II

    +
  1. +

    + New: FEValues objects can now directly calculate the Laplacians of + solution functions at quadrature points. This is equivalent to the + trace of the Hessians, but faster and more convenient in user code. + @ref step_31 "step-31" program has also been updated to use this. +
    + (Martin Kronbichler 2008/12/15) +

    +
  2. Improved: The FEValuesViews objects that one gets when writing things diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 2d1b94c356..7ef6eb8f98 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -597,8 +597,8 @@ class BoussinesqFlowProblem const std::vector &old_old_temperature, const std::vector > &old_temperature_grads, const std::vector > &old_old_temperature_grads, - const std::vector > &old_temperature_hessians, - const std::vector > &old_old_temperature_hessians, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, const std::vector > &old_velocity_values, const std::vector > &old_old_velocity_values, const std::vector &gamma_values, @@ -963,8 +963,8 @@ compute_viscosity (const std::vector &old_temperature, const std::vector &old_old_temperature, const std::vector > &old_temperature_grads, const std::vector > &old_old_temperature_grads, - const std::vector > &old_temperature_hessians, - const std::vector > &old_old_temperature_hessians, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, const std::vector > &old_velocity_values, const std::vector > &old_old_velocity_values, const std::vector &gamma_values, @@ -994,8 +994,8 @@ compute_viscosity (const std::vector &old_temperature, old_old_temperature_grads[q]) / 2; const double kappa_Delta_T = EquationData::kappa - * (trace(old_temperature_hessians[q]) + - trace(old_old_temperature_hessians[q])) / 2; + * (old_temperature_laplacians[q] + + old_old_temperature_laplacians[q]) / 2; const double residual = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) * @@ -1987,8 +1987,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector old_old_temperature_values(n_q_points); std::vector > old_temperature_grads(n_q_points); std::vector > old_old_temperature_grads(n_q_points); - std::vector > old_temperature_hessians(n_q_points); - std::vector > old_old_temperature_hessians(n_q_points); + std::vector old_temperature_laplacians(n_q_points); + std::vector old_old_temperature_laplacians(n_q_points); EquationData::TemperatureRightHandSide temperature_right_hand_side; std::vector gamma_values (n_q_points); @@ -2039,10 +2039,10 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads); - temperature_fe_values.get_function_hessians (old_temperature_solution, - old_temperature_hessians); - temperature_fe_values.get_function_hessians (old_old_temperature_solution, - old_old_temperature_hessians); + temperature_fe_values.get_function_laplacians (old_temperature_solution, + old_temperature_laplacians); + temperature_fe_values.get_function_laplacians (old_old_temperature_solution, + old_old_temperature_laplacians); temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); @@ -2079,8 +2079,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_values, old_temperature_grads, old_old_temperature_grads, - old_temperature_hessians, - old_old_temperature_hessians, + old_temperature_laplacians, + old_old_temperature_laplacians, old_velocity_values, old_old_velocity_values, gamma_values, diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index affab19e45..73ee39434f 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -284,8 +284,8 @@ class BoussinesqFlowProblem const std::vector &old_old_temperature, const std::vector > &old_temperature_grads, const std::vector > &old_old_temperature_grads, - const std::vector > &old_temperature_hessians, - const std::vector > &old_old_temperature_hessians, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, const std::vector > &present_velocity_values, const std::vector &gamma_values, const double global_u_infty, @@ -465,7 +465,7 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const double min_temperature, max_temperature; trilinos_communicator.MaxAll(&max_local_temperature, &max_temperature, 1); - trilinos_communicator.MaxAll(&min_local_temperature, &min_temperature, 1); + trilinos_communicator.MinAll(&min_local_temperature, &min_temperature, 1); return std::make_pair(min_temperature, max_temperature); } @@ -497,7 +497,7 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const double min_temperature, max_temperature; trilinos_communicator.MaxAll(&max_local_temperature, &max_temperature, 1); - trilinos_communicator.MaxAll(&min_local_temperature, &min_temperature, 1); + trilinos_communicator.MinAll(&min_local_temperature, &min_temperature, 1); return std::make_pair(min_temperature, max_temperature); } @@ -512,8 +512,8 @@ compute_viscosity(const std::vector &old_temperature, const std::vector &old_old_temperature, const std::vector > &old_temperature_grads, const std::vector > &old_old_temperature_grads, - const std::vector > &old_temperature_hessians, - const std::vector > &old_old_temperature_hessians, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplaceians, const std::vector > &present_velocity_values, const std::vector &gamma_values, const double global_u_infty, @@ -544,8 +544,8 @@ compute_viscosity(const std::vector &old_temperature, old_old_temperature_grads[q]) / 2; const double kappa_Delta_T = EquationData::kappa - * (trace(old_temperature_hessians[q]) + - trace(old_old_temperature_hessians[q])) / 2; + * (old_temperature_laplacians[q] + + old_old_temperature_laplacians[q]) / 2; const double residual = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) * @@ -1069,8 +1069,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector old_old_temperature_values(n_q_points); std::vector > old_temperature_grads(n_q_points); std::vector > old_old_temperature_grads(n_q_points); - std::vector > old_temperature_hessians(n_q_points); - std::vector > old_old_temperature_hessians(n_q_points); + std::vector old_temperature_laplacians(n_q_points); + std::vector old_old_temperature_laplacians(n_q_points); EquationData::TemperatureRightHandSide temperature_right_hand_side; @@ -1113,9 +1113,9 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_grads); temperature_fe_values.get_function_hessians (old_temperature_solution, - old_temperature_hessians); + old_temperature_laplacians); temperature_fe_values.get_function_hessians (old_old_temperature_solution, - old_old_temperature_hessians); + old_old_temperature_laplacians); temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); @@ -1128,8 +1128,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_values, old_temperature_grads, old_old_temperature_grads, - old_temperature_hessians, - old_old_temperature_hessians, + old_temperature_laplacians, + old_old_temperature_laplacians, present_velocity_values, gamma_values, global_u_infty,