From af4207dfbd92f31179be1c93c498e37e6fd896dd Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 13 Aug 2008 22:05:24 +0000 Subject: [PATCH] Add the functionality inverse_jacobian to the FEValues class, which gives the inverse Jacobian of the transformation from the real to the unit cell. git-svn-id: https://svn.dealii.org/trunk@16539 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_update_flags.h | 27 ++++++---- deal.II/deal.II/include/fe/fe_values.h | 51 +++++++++++++++++-- deal.II/deal.II/include/fe/mapping.h | 11 ++-- .../deal.II/include/fe/mapping_cartesian.h | 3 +- deal.II/deal.II/include/fe/mapping_q.h | 3 +- deal.II/deal.II/include/fe/mapping_q1.h | 3 +- deal.II/deal.II/source/fe/fe_values.cc | 6 ++- .../deal.II/source/fe/mapping_cartesian.cc | 13 ++++- deal.II/deal.II/source/fe/mapping_q.cc | 5 +- deal.II/deal.II/source/fe/mapping_q1.cc | 30 ++++++++--- 10 files changed, 122 insertions(+), 30 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h index 5d6711e033..027ee86bd1 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -165,6 +165,15 @@ enum UpdateFlags * transformation. */ update_jacobian_grads = 0x0100, + //! Volume element + /** + * Compute the inverse + * Jacobian of the + * transformation from the + * reference cell to the real + * cell. + */ + update_inverse_jacobians = 0x0200, //! Covariant transformation /** * Compute all values the @@ -174,9 +183,9 @@ enum UpdateFlags * mappings like * MappingCartesian this may be * simpler than - * #update_jacobians. + * #update_inverse_jacobians. */ - update_covariant_transformation = 0x0200, + update_covariant_transformation = 0x0400, //! Contravariant transformation /** * Compute all values the @@ -188,14 +197,14 @@ enum UpdateFlags * simpler than * #update_jacobians. */ - update_contravariant_transformation = 0x0400, + update_contravariant_transformation = 0x0800, //! Shape function values of transformation /** * Compute the shape function * values of the transformation * defined by the Mapping. */ - update_transformation_values = 0x0800, + update_transformation_values = 0x1000, //! Shape function gradients of transformation /** * Compute the shape function @@ -203,7 +212,7 @@ enum UpdateFlags * transformation defined by * the Mapping. */ - update_transformation_gradients = 0x1000, + update_transformation_gradients = 0x2000, /** * Compute the JxW values @@ -214,25 +223,25 @@ enum UpdateFlags * is used in conjunction with * H_div subspaces like RT and ABF. */ - update_cell_JxW_values = 0x2000, + update_cell_JxW_values = 0x4000, /** * Update the location of the * mapped generalized support * points of the element. */ - update_support_points = 0x4000, + update_support_points = 0x8000, /** * Update the Jacobian of the * mapping in generalized * support points. */ - update_support_jacobians = 0x8000, + update_support_jacobians = 0x10000, /** * Update the inverse Jacobian * of the mapping in * generalized support points. */ - update_support_inverse_jacobians = 0x10000, + update_support_inverse_jacobians = 0x20000, /** * @deprecated Update * quadrature points diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 54edbc94e9..cc843b40ca 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -562,13 +562,19 @@ class FEValuesData * quadrature points. */ std::vector > jacobians; - + /** * Array of the derivatives of the Jacobian * matrices at the quadrature points. */ std::vector > jacobian_grads; - + + /** + * Array of the inverse Jacobian matrices + * at the quadrature points. + */ + std::vector > inverse_jacobians; + /** * Store an array of weights * times the Jacobi determinant @@ -1633,7 +1639,22 @@ class FEValuesBase : protected FEValuesData, * jacobian_grads(). */ const std::vector > & get_jacobian_grads () const; - + + /** + * Return the inverse Jacobian of the + * transformation at the specified + * quadrature point, i.e. + * $J_{ij}=d\hat x_i/dx_j$ + */ + const Tensor<2,dim> & inverse_jacobian (const unsigned int quadrature_point) const; + + /** + * Pointer to the array holding + * the values returned by + * inverse_jacobian(). + */ + const std::vector > & get_inverse_jacobians () const; + /** * Constant reference to the * selected mapping object. @@ -4038,6 +4059,17 @@ FEValuesBase::get_jacobian_grads () const +template +inline +const std::vector >& +FEValuesBase::get_inverse_jacobians () const +{ + Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField()); + return this->inverse_jacobians; +} + + + template inline const Point & @@ -4091,6 +4123,19 @@ FEValuesBase::jacobian_grad (const unsigned int i) const +template +inline +const Tensor<2,dim> & +FEValuesBase::inverse_jacobian (const unsigned int i) const +{ + Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField()); + Assert (iinverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size())); + + return this->inverse_jacobians[i]; +} + + + template template inline diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 88c3992a00..efc90904bf 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -484,9 +484,11 @@ class Mapping : public Subscriptor * matrices needed to transform * vector-valued functions, * namely - * @p jacobians - * and the derivatives - * @p jacobian_grads. + * @p jacobians, + * the derivatives + * @p jacobian_grads, + * and the inverse operations in + * @p inverse_jacobians. */ virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, @@ -495,7 +497,8 @@ class Mapping : public Subscriptor std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const = 0; + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const = 0; /** * Performs the same as @p fill_fe_values diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index d9bf66f544..2ef9a7d79f 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -61,7 +61,8 @@ class MappingCartesian : public Mapping std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const ; + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const ; virtual void fill_fe_face_values (const typename Triangulation::cell_iterator &cell, diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 2f903a3e4f..53b3290ba6 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -233,7 +233,8 @@ class MappingQ : public MappingQ1 typename std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const ; + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const ; /** * Implementation of the interface in diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index c600375b5e..293609b83c 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -309,7 +309,8 @@ class MappingQ1 : public Mapping typename std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const ; + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const ; /** * Implementation of the interface in diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 991cd04fd8..8c619ec4c8 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -314,6 +314,9 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, if (flags & update_jacobian_grads) this->jacobian_grads.resize(n_quadrature_points); + if (flags & update_inverse_jacobians) + this->inverse_jacobians.resize(n_quadrature_points); + if (flags & update_boundary_forms) this->boundary_forms.resize(n_quadrature_points); @@ -1386,7 +1389,8 @@ void FEValues::do_reinit () this->quadrature_points, this->JxW_values, this->jacobians, - this->jacobian_grads); + this->jacobian_grads, + this->inverse_jacobians); this->get_fe().fill_fe_values(this->get_mapping(), *this->present_cell, diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 323b5afc24..b36b23d004 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -303,7 +303,8 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, std::vector >& quadrature_points, std::vector& JxW_values, std::vector >& jacobians, - std::vector >& jacobian_grads) const + std::vector >& jacobian_grads, + std::vector >& inverse_jacobians) const { // convert data object to internal // data for this class. fails with @@ -346,6 +347,16 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, if (data.current_update_flags() & update_jacobian_grads) for (unsigned int i=0; i(); + // "compute" inverse Jacobian at the + // quadrature points, which are all + // the same + if (data.current_update_flags() & update_inverse_jacobians) + for (unsigned int i=0; i(); + for (unsigned int j=0; j::fill_fe_values (const typename Triangulation::cell_iterator std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const { // convert data object to internal // data for this class. fails with @@ -334,7 +335,7 @@ MappingQ::fill_fe_values (const typename Triangulation::cell_iterator MappingQ1::fill_fe_values(cell, q, *p_data, quadrature_points, JxW_values, - jacobians, jacobian_grads); + jacobians, jacobian_grads, inverse_jacobians); } diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 877748f947..fec6d69eb4 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -343,7 +343,8 @@ MappingQ1::update_once (const UpdateFlags in) const | update_boundary_forms | update_normal_vectors | update_jacobians - | update_jacobian_grads)) + | update_jacobian_grads + | update_inverse_jacobians)) out |= update_transformation_gradients; return out; @@ -365,7 +366,8 @@ MappingQ1::update_each (const UpdateFlags in) const | update_boundary_forms | update_normal_vectors | update_jacobians - | update_jacobian_grads)); + | update_jacobian_grads + | update_inverse_jacobians)); // add a few flags. note that some // flags appear in both conditions @@ -373,10 +375,10 @@ MappingQ1::update_each (const UpdateFlags in) const // operations. this leads to some // circular logic. the only way to // treat this is to iterate. since - // there are 3 if-clauses in the + // there are 4 if-clauses in the // loop, it will take at most 3 // iterations to converge. do them: - for (unsigned int i=0; i<3; ++i) + for (unsigned int i=0; i<4; ++i) { // The following is a little incorrect: // If not applied on a face, @@ -398,6 +400,9 @@ MappingQ1::update_each (const UpdateFlags in) const | update_boundary_forms | update_normal_vectors)) out |= update_contravariant_transformation; + + if (out & (update_inverse_jacobians)) + out |= update_covariant_transformation; // The contravariant transformation // is a Piola transformation, which @@ -672,12 +677,13 @@ MappingQ1::compute_mapping_support_points( template void MappingQ1::fill_fe_values (const typename Triangulation::cell_iterator &cell, - const Quadrature &q, - typename Mapping::InternalDataBase &mapping_data, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &jacobians, - std::vector > &jacobian_grads) const + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians) const { InternalData *data_ptr = dynamic_cast (&mapping_data); Assert(data_ptr!=0, ExcInternalError()); @@ -733,6 +739,16 @@ MappingQ1::fill_fe_values (const typename Triangulation::cell_iterator += (data.second_derivative(point+DataSetDescriptor::cell (), k)[j][l] * data.mapping_support_points[k][i]); + } + // copy values from InternalData to vector + // given by reference + if (update_flags & update_inverse_jacobians) + { + Assert (inverse_jacobians.size() == n_q_points, + ExcDimensionMismatch(inverse_jacobians.size(), n_q_points)); + for (unsigned int point=0; point