From f5af5b1cb4712ecf8665a725fd785bae61c8f74a Mon Sep 17 00:00:00 2001 From: Maien Hamed Date: Mon, 10 Aug 2015 23:16:06 +0200 Subject: [PATCH] Modify hessian computations from numerical to analytical differentiation --- include/deal.II/fe/fe_poly.h | 35 +++++- include/deal.II/fe/fe_poly.templates.h | 141 ++++++++++++++++++------- include/deal.II/fe/fe_system.h | 11 +- source/fe/fe_poly.cc | 141 ++++++++++++++++--------- source/fe/fe_system.cc | 118 +++------------------ source/fe/mapping_cartesian.cc | 77 ++++++++++++++ 6 files changed, 319 insertions(+), 204 deletions(-) diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h index cbbb80b7ae..cd962adbdf 100644 --- a/include/deal.II/fe/fe_poly.h +++ b/include/deal.II/fe/fe_poly.h @@ -248,7 +248,6 @@ protected: */ virtual UpdateFlags update_each (const UpdateFlags flags) const; - /** * Fields of cell-independent data. * @@ -279,8 +278,42 @@ protected: * multiplication) when visiting an actual cell. */ std::vector > > shape_gradients; + + /** + * Array with shape function hessians in quadrature points. There is one + * row for each shape function, containing values for each quadrature + * point. + * + * We store the hessians in the quadrature points on the unit cell. We + * then only have to apply the transformation when visiting an actual cell. + */ + std::vector > > shape_hessians; + + /** + * Scratch array to store temporary values during hessian calculations in + * actual cells. + */ + mutable std::vector > untransformed_shape_hessians; }; + /** + * 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: + * + * \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) + * + * After this correction, the shape hessians are simply a mapping_covariant_gradient + * transformation. + */ + void + correct_untransformed_hessians (VectorSlice< std::vector > > uncorrected_shape_hessians, + const internal::FEValues::MappingRelatedData &mapping_data, + const internal::FEValues::FiniteElementRelatedData &fevalues_data, + const unsigned int n_q_points, + const unsigned int dof) const; + /** * The polynomial space. Its type is given by the template parameter POLY. */ diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index cb751ef28d..d93c690e58 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -150,7 +150,8 @@ FE_Poly::update_each (const UpdateFlags flags) const if (flags & update_gradients) out |= update_gradients | update_covariant_transformation; if (flags & update_hessians) - out |= update_hessians | update_covariant_transformation; + out |= update_hessians | update_covariant_transformation + | update_gradients | update_jacobian_grads; if (flags & update_cell_normal_vectors) out |= update_cell_normal_vectors | update_JxW_values; @@ -166,7 +167,7 @@ FE_Poly::update_each (const UpdateFlags flags) const template typename FiniteElement::InternalDataBase * FE_Poly::get_data (const UpdateFlags update_flags, - const Mapping &mapping, + const Mapping &, const Quadrature &quadrature) const { // generate a new data object and @@ -211,7 +212,12 @@ FE_Poly::get_data (const UpdateFlags update_flags, // then initialize some objects for // that if (flags & update_hessians) - data->initialize_2nd (this, mapping, quadrature); + { + grad_grads.resize (this->dofs_per_cell); + data->shape_hessians.resize (this->dofs_per_cell, + std::vector > (n_q_points)); + data->untransformed_shape_hessians.resize (n_q_points); + } // next already fill those fields // of which we have information by @@ -220,7 +226,7 @@ FE_Poly::get_data (const UpdateFlags update_flags, // unit cell, and need to be // transformed when visiting an // actual cell - if (flags & (update_values | update_gradients)) + if (flags & (update_values | update_gradients | update_hessians)) for (unsigned int i=0; i::get_data (const UpdateFlags update_flags, if (flags & update_gradients) for (unsigned int k=0; kdofs_per_cell; ++k) data->shape_gradients[k][i] = grads[k]; + + if (flags & update_hessians) + for (unsigned int k=0; kdofs_per_cell; ++k) + data->shape_hessians[k][i] = grad_grads[k]; } return data; } @@ -248,14 +258,14 @@ FE_Poly::get_data (const UpdateFlags update_flags, template void FE_Poly:: -fill_fe_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - const typename Mapping::InternalDataBase &mapping_internal, +fill_fe_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &mapping_internal, const typename FiniteElement::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData &, - internal::FEValues::FiniteElementRelatedData &output_data, - const CellSimilarity::Similarity cell_similarity) const + const internal::FEValues::MappingRelatedData &mapping_data, + internal::FEValues::FiniteElementRelatedData &output_data, + const CellSimilarity::Similarity cell_similarity) const { // convert data object to internal // data for this class. fails with @@ -276,12 +286,23 @@ fill_fe_values (const Mapping &mapping, mapping.transform(fe_data.shape_gradients[k], output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, QProjector::DataSetDescriptor::cell(), - mapping_internal, fe_data, - output_data); + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) + { + // compute the hessians in the unit cell (accounting for the Jacobian gradiant) + for (unsigned int i=0; i &mapping, template void FE_Poly:: -fill_fe_face_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int face, - const Quadrature &quadrature, +fill_fe_face_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, const typename Mapping::InternalDataBase &mapping_internal, - const typename FiniteElement::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData &, - internal::FEValues::FiniteElementRelatedData &output_data) const + const typename FiniteElement::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData &mapping_data, + internal::FEValues::FiniteElementRelatedData &output_data) const { // convert data object to internal // data for this class. fails with @@ -328,27 +349,37 @@ fill_fe_face_values (const Mapping &mapping, mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()), output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians) - this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data, - output_data); + if (flags & update_hessians) + { + // compute the hessians in the unit cell (accounting for the Jacobian gradiant) + for (unsigned int i=0; i > > + ( fe_data.untransformed_shape_hessians, offset , quadrature.size()), + mapping_data, output_data, quadrature.size(), k); + + mapping.transform(make_slice(fe_data.untransformed_shape_hessians, offset, quadrature.size()), + output_data.shape_hessians[k], mapping_internal, mapping_covariant_gradient); + } + } } - - template void FE_Poly:: -fill_fe_subface_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int face, - const unsigned int subface, - const Quadrature &quadrature, +fill_fe_subface_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, const typename Mapping::InternalDataBase &mapping_internal, - const typename FiniteElement::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData &, - internal::FEValues::FiniteElementRelatedData &output_data) const + const typename FiniteElement::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData &mapping_data, + internal::FEValues::FiniteElementRelatedData &output_data) const { // convert data object to internal // data for this class. fails with @@ -381,14 +412,42 @@ fill_fe_subface_values (const Mapping &mapping, mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()), output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians) - this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data, - output_data); + if (flags & update_hessians) + { + // compute the hessians in the unit cell (accounting for the Jacobian gradiant) + for (unsigned int i=0; i > > + ( fe_data.untransformed_shape_hessians, offset , quadrature.size()), + mapping_data, output_data, quadrature.size(), k); + + mapping.transform(make_slice(fe_data.untransformed_shape_hessians, offset, quadrature.size()), + output_data.shape_hessians[k], mapping_internal, mapping_covariant_gradient); + } + } } +template +void +FE_Poly:: +correct_untransformed_hessians (VectorSlice< std::vector > > uncorrected_shape_hessians, + const internal::FEValues::MappingRelatedData &mapping_data, + const internal::FEValues::FiniteElementRelatedData &fevalues_data, + const unsigned int n_q_points, + const unsigned int dof) const +{ + for (unsigned int i=0; i #include #include - #include DEAL_II_NAMESPACE_OPEN @@ -30,14 +29,14 @@ DEAL_II_NAMESPACE_OPEN template <> void FE_Poly,1,2>:: -fill_fe_values (const Mapping<1,2> &mapping, - const Triangulation<1,2>::cell_iterator &cell, - const Quadrature<1> &quadrature, - const Mapping<1,2>::InternalDataBase &mapping_internal, - const FiniteElement<1,2>::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData<1,2> &, +fill_fe_values (const Mapping<1,2> &mapping, + const Triangulation<1,2>::cell_iterator &, + const Quadrature<1> &quadrature, + const Mapping<1,2>::InternalDataBase &mapping_internal, + const FiniteElement<1,2>::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData<1,2> &mapping_data, internal::FEValues::FiniteElementRelatedData<1,2> &output_data, - const CellSimilarity::Similarity cell_similarity) const + const CellSimilarity::Similarity cell_similarity) const { // convert data object to internal // data for this class. fails with @@ -58,12 +57,23 @@ fill_fe_values (const Mapping<1,2> &mapping, mapping.transform(fe_data.shape_gradients[k], output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, QProjector<1>::DataSetDescriptor::cell(), - mapping_internal, fe_data, - output_data); + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) + { + // compute the hessians in the unit cell (accounting for the Jacobian gradiant) + for (unsigned int i=0; i &mapping, template <> void FE_Poly,2,3>:: -fill_fe_values (const Mapping<2,3> &mapping, - const Triangulation<2,3>::cell_iterator &cell, - const Quadrature<2> &quadrature, - const Mapping<2,3>::InternalDataBase &mapping_internal, - const FiniteElement<2,3>::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData<2,3> &, +fill_fe_values (const Mapping<2,3> &mapping, + const Triangulation<2,3>::cell_iterator &, + const Quadrature<2> &quadrature, + const Mapping<2,3>::InternalDataBase &mapping_internal, + const FiniteElement<2,3>::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData<2,3> &mapping_data, internal::FEValues::FiniteElementRelatedData<2,3> &output_data, - const CellSimilarity::Similarity cell_similarity) const + const CellSimilarity::Similarity cell_similarity) const { // assert that the following dynamics @@ -98,26 +108,37 @@ fill_fe_values (const Mapping<2,3> &mapping, mapping.transform(fe_data.shape_gradients[k], output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, QProjector<2>::DataSetDescriptor::cell(), - mapping_internal, fe_data, - output_data); + 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 void FE_Poly,1,2>:: -fill_fe_values (const Mapping<1,2> &mapping, - const Triangulation<1,2>::cell_iterator &cell, - const Quadrature<1> &quadrature, - const Mapping<1,2>::InternalDataBase &mapping_internal, - const FiniteElement<1,2>::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData<1,2> &, +fill_fe_values (const Mapping<1,2> &mapping, + const Triangulation<1,2>::cell_iterator &, + const Quadrature<1> &quadrature, + const Mapping<1,2>::InternalDataBase &mapping_internal, + const FiniteElement<1,2>::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData<1,2> &mapping_data, internal::FEValues::FiniteElementRelatedData<1,2> &output_data, - const CellSimilarity::Similarity cell_similarity) const + const CellSimilarity::Similarity cell_similarity) const { // convert data object to internal // data for this class. fails with @@ -139,26 +160,37 @@ fill_fe_values (const Mapping<1,2> &mapping, mapping.transform(fe_data.shape_gradients[k], output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, QProjector<1>::DataSetDescriptor::cell(), - mapping_internal, fe_data, - output_data); + 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 void FE_Poly,2,3>:: -fill_fe_values (const Mapping<2,3> &mapping, - const Triangulation<2,3>::cell_iterator &cell, - const Quadrature<2> &quadrature, - const Mapping<2,3>::InternalDataBase &mapping_internal, - const FiniteElement<2,3>::InternalDataBase &fedata, - const internal::FEValues::MappingRelatedData<2,3> &, +fill_fe_values (const Mapping<2,3> &mapping, + const Triangulation<2,3>::cell_iterator &, + const Quadrature<2> &quadrature, + const Mapping<2,3>::InternalDataBase &mapping_internal, + const FiniteElement<2,3>::InternalDataBase &fedata, + const internal::FEValues::MappingRelatedData<2,3> &mapping_data, internal::FEValues::FiniteElementRelatedData<2,3> &output_data, - const CellSimilarity::Similarity cell_similarity) const + const CellSimilarity::Similarity cell_similarity) const { Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); const InternalData &fe_data = static_cast (fedata); @@ -176,12 +208,23 @@ fill_fe_values (const Mapping<2,3> &mapping, mapping.transform(fe_data.shape_gradients[k], output_data.shape_gradients[k], mapping_internal, mapping_covariant); - } - if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, QProjector<2>::DataSetDescriptor::cell(), - mapping_internal, fe_data, - output_data); + 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 -FESystem::InternalData::InternalData(const unsigned int n_base_elements, - const bool compute_hessians) +FESystem::InternalData::InternalData(const unsigned int n_base_elements) : - compute_hessians (compute_hessians), base_fe_datas(n_base_elements), base_fe_output_objects(n_base_elements) {} @@ -792,18 +790,6 @@ FESystem::update_each (const UpdateFlags flags) const for (unsigned int base_no=0; base_non_base_elements(); ++base_no) out |= base_element(base_no).update_each(flags); - // second derivatives are handled - // by the top-level finite element, - // rather than by the base elements - // since it is generated by finite - // differencing. if second - // derivatives are requested, we - // therefore have to set the - // respective flag since the base - // elements don't have them - if (flags & update_hessians) - out |= update_hessians | update_covariant_transformation; - return out; } @@ -816,27 +802,12 @@ FESystem::get_data (const UpdateFlags flags_, const Quadrature &quadrature) const { UpdateFlags flags = flags_; - InternalData *data = new InternalData(this->n_base_elements(), - flags & update_hessians); + InternalData *data = new InternalData(this->n_base_elements()); data->update_once = update_once (flags); data->update_each = update_each (flags); flags = data->update_once | data->update_each; - UpdateFlags sub_flags = flags; - // if second derivatives through - // finite differencing are required, - // then initialize some objects for - // that - if (data->compute_hessians) - { - // delete - // update_hessians - // from flags list - sub_flags = UpdateFlags (sub_flags ^ update_hessians); - data->initialize_2nd (this, mapping, quadrature); - } - // get data objects from each of // the base elements and store them for (unsigned int base_no=0; base_non_base_elements(); ++base_no) @@ -844,7 +815,7 @@ FESystem::get_data (const UpdateFlags flags_, // create internal objects for base elements and initialize their // respective output objects typename FiniteElement::InternalDataBase *base_fe_data = - base_element(base_no).get_data(sub_flags, mapping, quadrature); + base_element(base_no).get_data(flags, mapping, quadrature); internal::FEValues::FiniteElementRelatedData &base_fe_output_object = data->get_fe_output_object(base_no); @@ -853,14 +824,6 @@ FESystem::get_data (const UpdateFlags flags_, // then store the pointer to the base internal object data->set_fe_data(base_no, base_fe_data); - - // make sure that *we* compute - // second derivatives, base - // elements should not do it - Assert (!(base_fe_data->update_each & update_hessians), - ExcInternalError()); - Assert (!(base_fe_data->update_once & update_hessians), - ExcInternalError()); } data->update_flags = data->update_once | data->update_each; @@ -878,26 +841,18 @@ FESystem::get_face_data ( const Quadrature &quadrature) const { UpdateFlags flags = flags_; - InternalData *data = new InternalData(this->n_base_elements(), - flags & update_hessians); + InternalData *data = new InternalData(this->n_base_elements()); data->update_once = update_once (flags); data->update_each = update_each (flags); flags = data->update_once | data->update_each; - UpdateFlags sub_flags = flags; - if (data->compute_hessians) - { - sub_flags = UpdateFlags (sub_flags ^ update_hessians); - data->initialize_2nd (this, mapping, QProjector::project_to_all_faces(quadrature)); - } - for (unsigned int base_no=0; base_non_base_elements(); ++base_no) { // create internal objects for base elements and initialize their // respective output objects typename FiniteElement::InternalDataBase *base_fe_data = - base_element(base_no).get_face_data(sub_flags, mapping, quadrature); + base_element(base_no).get_face_data(flags, mapping, quadrature); internal::FEValues::FiniteElementRelatedData &base_fe_output_object = data->get_fe_output_object(base_no); @@ -906,11 +861,6 @@ FESystem::get_face_data ( // then store the pointer to the base internal object data->set_fe_data(base_no, base_fe_data); - - Assert (!(base_fe_data->update_each & update_hessians), - ExcInternalError()); - Assert (!(base_fe_data->update_once & update_hessians), - ExcInternalError()); } data->update_flags = data->update_once | data->update_each; @@ -930,26 +880,18 @@ FESystem::get_subface_data ( const Quadrature &quadrature) const { UpdateFlags flags = flags_; - InternalData *data = new InternalData(this->n_base_elements(), - flags & update_hessians); + InternalData *data = new InternalData(this->n_base_elements()); data->update_once = update_once (flags); data->update_each = update_each (flags); flags = data->update_once | data->update_each; - UpdateFlags sub_flags = flags; - if (data->compute_hessians) - { - sub_flags = UpdateFlags (sub_flags ^ update_hessians); - data->initialize_2nd (this, mapping, QProjector::project_to_all_subfaces(quadrature)); - } - for (unsigned int base_no=0; base_non_base_elements(); ++base_no) { // create internal objects for base elements and initialize their // respective output objects typename FiniteElement::InternalDataBase *base_fe_data = - base_element(base_no).get_subface_data(sub_flags, mapping, quadrature); + base_element(base_no).get_subface_data(flags, mapping, quadrature); internal::FEValues::FiniteElementRelatedData &base_fe_output_object = data->get_fe_output_object(base_no); @@ -958,11 +900,6 @@ FESystem::get_subface_data ( // then store the pointer to the base internal object data->set_fe_data(base_no, base_fe_data); - - Assert (!(base_fe_data->update_each & update_hessians), - ExcInternalError()); - Assert (!(base_fe_data->update_once & update_hessians), - ExcInternalError()); } data->update_flags = data->update_once | data->update_each; @@ -1161,10 +1098,11 @@ compute_fill_one_base (const Mapping &mapping output_data.shape_gradients[out_index+s][q] = base_data.shape_gradients[in_index+s][q]; - // _we_ handle computation of second derivatives, so the - // base elements should not have computed them! - Assert (!(base_flags & update_hessians), - ExcInternalError()); + if (base_flags & update_hessians) + for (unsigned int q=0; q &mapping, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, - const CellSimilarity::Similarity cell_similarity, + const CellSimilarity::Similarity , const typename Mapping::InternalDataBase &mapping_internal, const typename FiniteElement::InternalDataBase &fedata, const internal::FEValues::MappingRelatedData &mapping_data, internal::FEValues::FiniteElementRelatedData &output_data) const { - const unsigned int n_q_points = quadrature.size(); - // convert data object to internal // data for this class. fails with // an exception if that is not @@ -1204,7 +1140,7 @@ compute_fill (const Mapping &mapping, fe_data.update_flags); - if (flags & (update_values | update_gradients)) + if (flags & (update_values | update_gradients | update_hessians)) { // let base elements update the necessary data Threads::TaskGroup<> task_group; @@ -1223,32 +1159,6 @@ compute_fill (const Mapping &mapping, std_cxx11::ref(output_data)))); task_group.join_all(); } - - if (fe_data.compute_hessians && cell_similarity != CellSimilarity::translation) - { - unsigned int offset = 0; - if (face_no != invalid_face_number) - { - if (sub_no == invalid_face_number) - offset=QProjector::DataSetDescriptor - ::face(face_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - n_q_points); - else - offset=QProjector::DataSetDescriptor - ::subface(face_no, sub_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - n_q_points, - cell->subface_case(face_no)); - } - - this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data, - output_data); - } } diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 1df48ce7a4..89df6a5fb0 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -657,6 +657,21 @@ MappingCartesian::transform ( } case mapping_piola: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); @@ -695,6 +710,68 @@ MappingCartesian::transform ( switch (mapping_type) { + case mapping_covariant: + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + + for (unsigned int i=0; i