From 10cdecff505246b6d9dc87a534356b2c93d93e71 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 7 Mar 2009 15:16:42 +0000 Subject: [PATCH] Make the FEValues::reinit functions (much) more efficient in case we repeatedly work on the same cell - in such a case, the derivatives will not change and need to be regenerated. Introduce a new enum CellSimilarity for doing that and pass it down to the respective functions. git-svn-id: https://svn.dealii.org/trunk@18464 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 11 +- .../deal.II/include/fe/fe_dgp_nonparametric.h | 3 +- deal.II/deal.II/include/fe/fe_nedelec.h | 9 +- deal.II/deal.II/include/fe/fe_poly.h | 7 +- .../deal.II/include/fe/fe_poly.templates.h | 138 +++++++------ deal.II/deal.II/include/fe/fe_poly_tensor.h | 5 +- deal.II/deal.II/include/fe/fe_system.h | 14 +- deal.II/deal.II/include/fe/fe_values.h | 37 ++++ deal.II/deal.II/include/fe/mapping.h | 35 +++- .../deal.II/include/fe/mapping_cartesian.h | 16 +- deal.II/deal.II/include/fe/mapping_q.h | 1 + deal.II/deal.II/include/fe/mapping_q1.h | 4 +- .../deal.II/include/fe/mapping_q1_eulerian.h | 19 +- .../deal.II/include/fe/mapping_q_eulerian.h | 18 ++ .../deal.II/source/fe/fe_dgp_nonparametric.cc | 19 +- deal.II/deal.II/source/fe/fe_nedelec.cc | 18 +- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 19 +- deal.II/deal.II/source/fe/fe_system.cc | 160 ++++++++------- deal.II/deal.II/source/fe/fe_values.cc | 92 +++++++-- .../deal.II/source/fe/mapping_cartesian.cc | 87 ++++---- deal.II/deal.II/source/fe/mapping_q.cc | 49 +++-- deal.II/deal.II/source/fe/mapping_q1.cc | 188 ++++++++++-------- .../deal.II/source/fe/mapping_q1_eulerian.cc | 33 ++- .../deal.II/source/fe/mapping_q_eulerian.cc | 26 +++ 24 files changed, 655 insertions(+), 353 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 8be9203217..7fef35499f 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -2588,12 +2588,13 @@ class FiniteElement : public Subscriptor, * called for the same cell first! */ virtual void - fill_fe_values (const Mapping &mapping, + fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_internal, - typename Mapping::InternalDataBase &fe_internal, - FEValuesData &data) const = 0; + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData &data) const = 0; /** * Fill the fields of diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h index 41e92926eb..98b0c4f9e0 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -480,7 +480,8 @@ class FE_DGPNonparametric : public FiniteElement virtual void fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index 62cc2404f5..5270ded777 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -412,11 +412,12 @@ class FE_Nedelec : public FiniteElement * FiniteElement. */ virtual void - fill_fe_values (const Mapping &mapping, + fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_internal, - typename Mapping::InternalDataBase &fe_internal, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; /** diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h index ed1b4652ff..c014e6819a 100644 --- a/deal.II/deal.II/include/fe/fe_poly.h +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -218,12 +218,13 @@ class FE_Poly : public FiniteElement const Quadrature& quadrature) const ; virtual void - fill_fe_values (const Mapping &mapping, + fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, - FEValuesData& data) const; + FEValuesData& data) const; virtual void fill_fe_face_values (const Mapping &mapping, diff --git a/deal.II/deal.II/include/fe/fe_poly.templates.h b/deal.II/deal.II/include/fe/fe_poly.templates.h index 9c54d47526..9f979ba926 100644 --- a/deal.II/deal.II/include/fe/fe_poly.templates.h +++ b/deal.II/deal.II/include/fe/fe_poly.templates.h @@ -156,8 +156,8 @@ FE_Poly::update_each (const UpdateFlags flags) const template typename Mapping::InternalDataBase * FE_Poly::get_data (const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature &quadrature) const + const Mapping &mapping, + const Quadrature &quadrature) const { // generate a new data object and // initialize some fields @@ -238,12 +238,14 @@ FE_Poly::get_data (const UpdateFlags update_flags, template <> inline void -FE_Poly,1,2>::fill_fe_values (const Mapping<1,2> &mapping, - const Triangulation<1,2>::cell_iterator &/*cell*/, - const Quadrature<1> &quadrature, - Mapping<1,2>::InternalDataBase &mapping_data, - Mapping<1,2>::InternalDataBase &fedata, - FEValuesData<1,2> &data) const +FE_Poly,1,2>::fill_fe_values + (const Mapping<1,2> &mapping, + const Triangulation<1,2>::cell_iterator &/*cell*/, + const Quadrature<1> &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + Mapping<1,2>::InternalDataBase &mapping_data, + Mapping<1,2>::InternalDataBase &fedata, + FEValuesData<1,2> &data) const { // convert data object to internal // data for this class. fails with @@ -260,13 +262,13 @@ FE_Poly,1,2>::fill_fe_values (const Mapping<1,2> &ma for (unsigned int i=0; i::DataSetDescriptor dsd; - if (flags & update_hessians) + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) { AssertThrow(false, ExcNotImplemented()); /* this->compute_2nd (mapping, cell, dsd.cell(), */ @@ -279,12 +281,14 @@ FE_Poly,1,2>::fill_fe_values (const Mapping<1,2> &ma template <> inline void -FE_Poly,2,3>::fill_fe_values (const Mapping<2,3> &mapping, - const Triangulation<2,3>::cell_iterator &/* cell */, - const Quadrature<2> &quadrature, - Mapping<2,3>::InternalDataBase & mapping_data, - Mapping<2,3>::InternalDataBase &fedata, - FEValuesData<2,3> &data) const +FE_Poly,2,3>::fill_fe_values + (const Mapping<2,3> &mapping, + const Triangulation<2,3>::cell_iterator &/* cell */, + const Quadrature<2> &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + Mapping<2,3>::InternalDataBase & mapping_data, + Mapping<2,3>::InternalDataBase &fedata, + FEValuesData<2,3> &data) const { // assert that the following dynamics @@ -300,13 +304,13 @@ FE_Poly,2,3>::fill_fe_values (const Mapping<2,3> &ma for (unsigned int i=0; i::DataSetDescriptor dsd; */ - if (flags & update_hessians) + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) { AssertThrow(false, ExcNotImplemented()); /* this->compute_2nd (mapping, cell, dsd.cell(), */ @@ -323,6 +327,7 @@ FE_Poly,1,2>::fill_fe_values ( const Mapping<1,2>& mapping, const Triangulation<1,2>::cell_iterator& /*cell*/, const Quadrature<1>& quadrature, + const enum CellSimilarity::Similarity cell_similarity, Mapping<1,2>::InternalDataBase& mapping_data, Mapping<1,2>::InternalDataBase& fedata, FEValuesData<1,2>& data) const @@ -345,13 +350,13 @@ FE_Poly,1,2>::fill_fe_values ( // TODO: I would think this should work. Guido - if (flags & update_gradients) + if (flags & update_gradients && cell_similarity != CellSimilarity::translation) mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], mapping_data, mapping_covariant); } // const typename QProjector<1>::DataSetDescriptor dsd; - if (flags & update_hessians) + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) { AssertThrow(false, ExcNotImplemented()); /* this->compute_2nd (mapping, cell, dsd.cell(), */ @@ -363,12 +368,14 @@ FE_Poly,1,2>::fill_fe_values ( template <> inline void -FE_Poly,2,3>::fill_fe_values (const Mapping<2,3>& mapping, - const Triangulation<2,3>::cell_iterator &/* cell */, - const Quadrature<2>& quadrature, - Mapping<2,3>::InternalDataBase& mapping_data, - Mapping<2,3>::InternalDataBase &fedata, - FEValuesData<2,3> &data) const +FE_Poly,2,3>::fill_fe_values + (const Mapping<2,3>& mapping, + const Triangulation<2,3>::cell_iterator &/* cell */, + const Quadrature<2>& quadrature, + const enum CellSimilarity::Similarity cell_similarity, + Mapping<2,3>::InternalDataBase& mapping_data, + Mapping<2,3>::InternalDataBase &fedata, + FEValuesData<2,3> &data) const { Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); InternalData &fe_data = static_cast (fedata); @@ -382,13 +389,13 @@ FE_Poly,2,3>::fill_fe_values (const Mapping<2,3>& mapping, data.shape_values(k,i) = fe_data.shape_values[k][i]; - if (flags & update_gradients) + if (flags & update_gradients && cell_similarity != CellSimilarity::translation) mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], mapping_data, mapping_covariant); } /* const QProjector<2>::DataSetDescriptor dsd; */ - if (flags & update_hessians) + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) { AssertThrow(false, ExcNotImplemented()); /* this->compute_2nd (mapping, cell, dsd.cell(), */ @@ -400,18 +407,19 @@ FE_Poly,2,3>::fill_fe_values (const Mapping<2,3>& mapping, template void -FE_Poly::fill_fe_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const +FE_Poly::fill_fe_values + (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with // an exception if that is not - // possible - + // possible Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); InternalData &fe_data = static_cast (fedata); @@ -423,13 +431,13 @@ FE_Poly::fill_fe_values (const Mapping for (unsigned int i=0; i::DataSetDescriptor dsd; - if (flags & update_hessians) + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) this->compute_2nd (mapping, cell, dsd.cell(), mapping_data, fe_data, data); } @@ -439,13 +447,14 @@ FE_Poly::fill_fe_values (const Mapping template <> inline void -FE_Poly,1,2>::fill_fe_face_values (const Mapping<1,2> &, - const Triangulation<1,2>::cell_iterator &, - const unsigned int, - const Quadrature<0> &, - Mapping<1,2>::InternalDataBase &, - Mapping<1,2>::InternalDataBase &, - FEValuesData<1,2> &) const +FE_Poly,1,2>::fill_fe_face_values + (const Mapping<1,2> &, + const Triangulation<1,2>::cell_iterator &, + const unsigned int, + const Quadrature<0> &, + Mapping<1,2>::InternalDataBase &, + Mapping<1,2>::InternalDataBase &, + FEValuesData<1,2> &) const { AssertThrow(false, ExcNotImplemented()); } @@ -472,12 +481,12 @@ template <> inline void FE_Poly,1,2>::fill_fe_face_values (const Mapping<1,2> &, - const Triangulation<1,2>::cell_iterator &, - const unsigned int, - const Quadrature<0> &, - Mapping<1,2>::InternalDataBase &, - Mapping<1,2>::InternalDataBase &, - FEValuesData<1,2> &) const + const Triangulation<1,2>::cell_iterator &, + const unsigned int, + const Quadrature<0> &, + Mapping<1,2>::InternalDataBase &, + Mapping<1,2>::InternalDataBase &, + FEValuesData<1,2> &) const { AssertThrow(false, ExcNotImplemented()); } @@ -488,12 +497,12 @@ template <> inline void FE_Poly,2,3>::fill_fe_face_values (const Mapping<2,3> &, - const Triangulation<2,3>::cell_iterator &, - const unsigned int, - const Quadrature<1> &, - Mapping<2,3>::InternalDataBase &, - Mapping<2,3>::InternalDataBase &, - FEValuesData<2,3> &) const + const Triangulation<2,3>::cell_iterator &, + const unsigned int, + const Quadrature<1> &, + Mapping<2,3>::InternalDataBase &, + Mapping<2,3>::InternalDataBase &, + FEValuesData<2,3> &) const { AssertThrow(false, ExcNotImplemented()); } @@ -501,13 +510,14 @@ FE_Poly,2,3>::fill_fe_face_values (const Mapping<2,3> &, template void -FE_Poly::fill_fe_face_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int face, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const +FE_Poly:: +fill_fe_face_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h index e4dfabb33a..7d641c4450 100644 --- a/deal.II/deal.II/include/fe/fe_poly_tensor.h +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -234,9 +234,10 @@ class FE_PolyTensor : public FiniteElement const Quadrature& quadrature) const ; virtual void - fill_fe_values (const Mapping &mapping, + fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 28a1f90548..336f74e7cc 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -611,7 +611,8 @@ class FESystem : public FiniteElement virtual void fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, FEValuesData &data) const; @@ -679,14 +680,15 @@ class FESystem : public FiniteElement * sub_no!=invalid_face_no. */ template - void compute_fill (const Mapping &mapping, + void compute_fill (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const ; + FEValuesData &data) const ; private: diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 3e8b5e70c1..4d77e95177 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -2224,6 +2224,12 @@ class FEValuesBase : protected FEValuesData, */ const FiniteElement & get_fe () const; + /** + * Return the update flags set + * for this object. + */ + UpdateFlags get_update_flags () const; + /** * Return a triangulation * iterator to the current cell. @@ -2809,6 +2815,28 @@ class FEValues : public FEValuesBase */ void initialize (const UpdateFlags update_flags); + /** + * An enum variable that can store + * different states of the current cell + * in comparison to the previously + * visited cell. If wanted, additional + * states can be checked here and used + * in one of the methods used during + * reinit. + */ + enum CellSimilarity::Similarity cell_similarity; + + /** + * A function that checks whether the + * new cell is similar to the one + * previously used. Then, a significant + * amount of the data can be reused, + * e.g. the derivatives of the basis + * functions in real space, shape_grad. + */ + void + check_cell_similarity (const typename Triangulation::cell_iterator &cell); + /** * The reinit() functions do * only that part of the work @@ -4083,6 +4111,15 @@ FEValuesBase::get_mapping () const template inline +UpdateFlags +FEValuesBase::get_update_flags () const +{ + return this->update_flags; +} + + + +template const typename Triangulation::cell_iterator FEValuesBase::get_cell () const { diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 442bd5441b..cfbaebbc29 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -32,6 +32,24 @@ template class FEValues; template class FEFaceValues; template class FESubfaceValues; + /** + * An enum variable that can store + * different states of the current cell + * in comparison to the previously + * visited cell. If wanted, additional + * states can be checked here and used + * in one of the methods used during + * reinit. + */ +namespace CellSimilarity +{ + enum Similarity + { + no_similarity, + translation + }; +} + /** * The transformation type used * for the Mapping::transform() functions. @@ -545,14 +563,15 @@ class Mapping : public Subscriptor */ virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - InternalDataBase &internal, - std::vector > &quadrature_points, - std::vector &JxW_values, - std::vector > &jacobians, - std::vector > &jacobian_grads, - std::vector > &inverse_jacobians, - std::vector > &cell_normal_vectors + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + InternalDataBase &internal, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors ) const=0; diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index dff1212763..9b094b6ef8 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -60,13 +60,14 @@ class MappingCartesian : public Mapping virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, - std::vector &JxW_values, - std::vector > &jacobians, - std::vector > &jacobian_grads, - std::vector > &inverse_jacobians, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, std::vector > &) const ; @@ -180,6 +181,7 @@ class MappingCartesian : public Mapping void compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, + const enum CellSimilarity::Similarity cell_similarity, InternalData& data, std::vector > &quadrature_points, std::vector >& normal_vectors) const; diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 24c8d11e60..038144041b 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -205,6 +205,7 @@ class MappingQ : public MappingQ1 virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, + enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename std::vector > &quadrature_points, std::vector &JxW_values, diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index cfc67ab395..f21b8eee9a 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -294,7 +294,8 @@ class MappingQ1 : public Mapping */ virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, - const Quadrature & quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename std::vector > &quadrature_points, std::vector &JxW_values, @@ -384,6 +385,7 @@ class MappingQ1 : public Mapping void compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int npts, const DataSetDescriptor data_set, + const enum CellSimilarity::Similarity cell_similarity, InternalData &data, std::vector > &quadrature_points) const; diff --git a/deal.II/deal.II/include/fe/mapping_q1_eulerian.h b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h index 949f37c599..540c377bfc 100644 --- a/deal.II/deal.II/include/fe/mapping_q1_eulerian.h +++ b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -143,6 +143,23 @@ class MappingQ1Eulerian : public MappingQ1 protected: + /** + * Implementation of the interface in + * MappingQ1. Overrides the function in + * the base class, since we cannot use + * any cell similarity for this class. + */ + virtual void + fill_fe_values (const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors) const; /** * Reference to the vector of diff --git a/deal.II/deal.II/include/fe/mapping_q_eulerian.h b/deal.II/deal.II/include/fe/mapping_q_eulerian.h index c5517c31f6..809bab67f8 100644 --- a/deal.II/deal.II/include/fe/mapping_q_eulerian.h +++ b/deal.II/deal.II/include/fe/mapping_q_eulerian.h @@ -146,6 +146,24 @@ class MappingQEulerian : public MappingQ << "-- expected size " << arg2); protected: + /** + * Implementation of the interface in + * MappingQ. Overrides the function in + * the base class, since we cannot use + * any cell similarity for this class. + */ + virtual void + fill_fe_values (const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors) const; + /** * Reference to the vector of * shifts. diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc index 5cf951458a..e80074c73a 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -289,6 +289,7 @@ FE_DGPNonparametric::fill_fe_values ( const Mapping&, const typename Triangulation::cell_iterator&, const Quadrature&, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase&, typename Mapping::InternalDataBase& fedata, FEValuesData& data) const @@ -297,13 +298,15 @@ FE_DGPNonparametric::fill_fe_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); - + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); + const UpdateFlags flags(fe_data.current_update_flags()); Assert (flags & update_quadrature_points, ExcInternalError()); - + const unsigned int n_q_points = data.quadrature_points.size(); - + if (flags & (update_values | update_gradients)) for (unsigned int i=0; i::fill_fe_face_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); Assert (flags & update_quadrature_points, ExcInternalError()); @@ -380,7 +385,9 @@ FE_DGPNonparametric::fill_fe_subface_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); Assert (flags & update_quadrature_points, ExcInternalError()); diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 71901c7b08..75ad0b8d5a 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -990,18 +990,22 @@ FE_Nedelec::get_data (const UpdateFlags update_flags, template void -FE_Nedelec::fill_fe_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const +FE_Nedelec::fill_fe_values + (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); // get the flags indicating the // fields that have to be filled diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index 04aa19b501..ebc257d682 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -359,18 +359,21 @@ FE_PolyTensor::get_data ( template void FE_PolyTensor::fill_fe_values ( - const Mapping &mapping, + const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + FEValuesData &data) const { // convert data object to internal // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); const unsigned int n_q_points = quadrature.size(); const UpdateFlags flags(fe_data.current_update_flags()); @@ -513,7 +516,9 @@ FE_PolyTensor::fill_fe_face_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); const unsigned int n_q_points = quadrature.size(); // offset determines which data set @@ -666,7 +671,9 @@ FE_PolyTensor::fill_fe_subface_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, + ExcInternalError()); + InternalData &fe_data = static_cast (fedata); const unsigned int n_q_points = quadrature.size(); // offset determines which data set diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 0290b28e69..9989fb6ea8 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -146,7 +146,7 @@ const unsigned int FESystem::invalid_face_number; template FESystem::FESystem (const FiniteElement &fe, - const unsigned int n_elements) : + const unsigned int n_elements) : FiniteElement (multiply_dof_numbers(fe, n_elements), compute_restriction_is_additive_flags (fe, n_elements), compute_nonzero_components(fe, n_elements)), @@ -161,9 +161,9 @@ FESystem::FESystem (const FiniteElement &fe, template FESystem::FESystem (const FiniteElement &fe1, - const unsigned int n1, - const FiniteElement &fe2, - const unsigned int n2) : + const unsigned int n1, + const FiniteElement &fe2, + const unsigned int n2) : FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2), compute_restriction_is_additive_flags (fe1, n1, fe2, n2), @@ -183,12 +183,12 @@ FESystem::FESystem (const FiniteElement &fe1, template FESystem::FESystem (const FiniteElement &fe1, - const unsigned int n1, - const FiniteElement &fe2, - const unsigned int n2, - const FiniteElement &fe3, - const unsigned int n3) : - FiniteElement (multiply_dof_numbers(fe1, n1, + const unsigned int n1, + const FiniteElement &fe2, + const unsigned int n2, + const FiniteElement &fe3, + const unsigned int n3) : + FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2, fe3, n3), compute_restriction_is_additive_flags (fe1, n1, @@ -288,7 +288,7 @@ FESystem::clone() const template double FESystem::shape_value (const unsigned int i, - const Point &p) const + const Point &p) const { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), @@ -303,8 +303,8 @@ FESystem::shape_value (const unsigned int i, template double FESystem::shape_value_component (const unsigned int i, - const Point &p, - const unsigned int component) const + const Point &p, + const unsigned int component) const { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component < this->n_components(), @@ -340,7 +340,7 @@ FESystem::shape_value_component (const unsigned int i, template Tensor<1,dim> FESystem::shape_grad (const unsigned int i, - const Point &p) const + const Point &p) const { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), @@ -392,7 +392,7 @@ FESystem::shape_grad_component (const unsigned int i, template Tensor<2,dim> FESystem::shape_grad_grad (const unsigned int i, - const Point &p) const + const Point &p) const { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), @@ -577,8 +577,8 @@ FESystem::update_each (const UpdateFlags flags) const template typename Mapping::InternalDataBase * FESystem::get_data (const UpdateFlags flags_, - const Mapping &mapping, - const Quadrature &quadrature) const + const Mapping &mapping, + const Quadrature &quadrature) const { UpdateFlags flags = flags_; @@ -664,15 +664,16 @@ FESystem::get_data (const UpdateFlags flags_, template void FESystem:: -fill_fe_values (const Mapping &mapping, +fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const + FEValuesData &data) const { compute_fill(mapping, cell, invalid_face_number, invalid_face_number, - quadrature, mapping_data, fe_data, data); + quadrature, cell_similarity, mapping_data, fe_data, data); } @@ -688,8 +689,8 @@ fill_fe_face_values (const Mapping &mapping, typename Mapping::InternalDataBase &fe_data, FEValuesData &data) const { - compute_fill (mapping, cell, face_no, invalid_face_number, - quadrature, mapping_data, fe_data, data); + compute_fill (mapping, cell, face_no, invalid_face_number, quadrature, + CellSimilarity::no_similarity, mapping_data, fe_data, data); } @@ -698,17 +699,17 @@ fill_fe_face_values (const Mapping &mapping, template void FESystem:: -fill_fe_subface_values (const Mapping &mapping, +fill_fe_subface_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const + FEValuesData &data) const { - compute_fill (mapping, cell, face_no, sub_no, - quadrature, mapping_data, fe_data, data); + compute_fill (mapping, cell, face_no, sub_no, quadrature, + CellSimilarity::no_similarity, mapping_data, fe_data, data); } @@ -717,14 +718,15 @@ template template void FESystem:: -compute_fill (const Mapping &mapping, +compute_fill (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + FEValuesData &data) const { const unsigned int n_q_points = quadrature.size(); @@ -732,7 +734,8 @@ compute_fill (const Mapping &mapping, // data for this class. fails with // an exception if that is not // possible - InternalData & fe_data = dynamic_cast (fedata); + Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); + InternalData & fe_data = static_cast (fedata); // Either dim_1==dim // (fill_fe_values) or dim_1==dim-1 @@ -800,16 +803,18 @@ compute_fill (const Mapping &mapping, if (face_no==invalid_face_number) { Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim)); + Assert (dynamic_cast *>(quadrature_base_pointer) != 0, + ExcInternalError()); cell_quadrature - = dynamic_cast *>(quadrature_base_pointer); - Assert (cell_quadrature != 0, ExcInternalError()); + = static_cast *>(quadrature_base_pointer); } else { Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1)); + Assert (dynamic_cast *>(quadrature_base_pointer) != 0, + ExcInternalError()); face_quadrature - = dynamic_cast *>(quadrature_base_pointer); - Assert (face_quadrature != 0, ExcInternalError()); + = static_cast *>(quadrature_base_pointer); } // let base elements update the @@ -852,8 +857,8 @@ compute_fill (const Mapping &mapping, // therefore use // base_fe_data.update_flags. if (face_no==invalid_face_number) - base_fe.fill_fe_values(mapping, cell, - *cell_quadrature, mapping_data, base_fe_data, base_data); + base_fe.fill_fe_values(mapping, cell, *cell_quadrature, cell_similarity, + mapping_data, base_fe_data, base_data); else if (sub_no==invalid_face_number) base_fe.fill_fe_face_values(mapping, cell, face_no, *face_quadrature, mapping_data, base_fe_data, base_data); @@ -892,14 +897,21 @@ compute_fill (const Mapping &mapping, // base element const UpdateFlags base_flags(dim_1==dim ? base_fe_data.current_update_flags() : - base_fe_data.update_flags); - for (unsigned int system_index=0; system_indexdofs_per_cell; - ++system_index) - if (this->system_to_base_table[system_index].first.first == base_no) - { - const unsigned int - base_index = this->system_to_base_table[system_index].second; - Assert (base_indexdofs_per_cell; + ++system_index) + if (this->system_to_base_table[system_index].first.first == base_no) + { + const unsigned int + base_index = this->system_to_base_table[system_index].second; + Assert (base_index &mapping, // this one value, and // to which index to // put - unsigned int out_index = 0; - for (unsigned int i=0; in_nonzero_components(i); - unsigned int in_index = 0; - for (unsigned int i=0; in_nonzero_components(i); + unsigned int in_index = 0; + for (unsigned int i=0; in_nonzero_components(system_index) == - base_fe.n_nonzero_components(base_index), - ExcInternalError()); - for (unsigned int s=0; sn_nonzero_components(system_index); ++s) - { - if (base_flags & update_values) - for (unsigned int q=0; qn_nonzero_components(system_index) == + base_fe.n_nonzero_components(base_index), + ExcInternalError()); + for (unsigned int s=0; sn_nonzero_components(system_index); ++s) + { + if (base_flags & update_values) + for (unsigned int q=0; q &mapping, // should not // have computed // them! - Assert (!(base_flags & update_hessians), - ExcInternalError()); - }; - }; + Assert (!(base_flags & update_hessians), + ExcInternalError()); + }; + }; }; if (fe_data.is_first_cell()) @@ -981,7 +993,7 @@ compute_fill (const Mapping &mapping, } } - if (fe_data.compute_hessians) + if (fe_data.compute_hessians && cell_similarity != CellSimilarity::translation) { unsigned int offset = 0; if (face_no != invalid_face_number) diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index f70ba8aa1c..0cb2bd6f11 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -1690,6 +1690,7 @@ FEValuesData::initialize (const unsigned int n_quadrature_p } + /*------------------------------- FEValuesBase ---------------------------*/ @@ -3238,7 +3239,6 @@ FEValuesBase::compute_update_flags (const UpdateFlags update_flags } - /*------------------------------- FEValues -------------------------------*/ @@ -3252,9 +3252,9 @@ const unsigned int FEValues::integral_dimension; template FEValues::FEValues (const Mapping &mapping, - const FiniteElement &fe, - const Quadrature &q, - const UpdateFlags update_flags) + const FiniteElement &fe, + const Quadrature &q, + const UpdateFlags update_flags) : FEValuesBase (q.size(), fe.dofs_per_cell, @@ -3270,15 +3270,15 @@ FEValues::FEValues (const Mapping &mapping, template FEValues::FEValues (const FiniteElement &fe, - const Quadrature &q, - const UpdateFlags update_flags) + const Quadrature &q, + const UpdateFlags update_flags) : FEValuesBase (q.size(), - fe.dofs_per_cell, - update_default, + fe.dofs_per_cell, + update_default, StaticMappingQ1::mapping, - fe), - quadrature (q) + fe), + quadrature (q) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); initialize (update_flags); @@ -3330,9 +3330,11 @@ FEValues::reinit (const typename DoFHandler::cell_it typedef FEValuesBase FEVB; Assert (static_cast&>(*this->fe) == - static_cast&>(cell->get_fe()), + static_cast&>(cell->get_fe()), typename FEVB::ExcFEDontMatch()); + check_cell_similarity(cell); + // set new cell. auto_ptr will take // care that old object gets // destroyed and also that this @@ -3360,11 +3362,13 @@ FEValues::reinit (const typename hp::DoFHandler::cel // passed to the constructor and // used by the DoFHandler used by // this cell, are the same - typedef FEValuesBase FEVB; + typedef FEValuesBase FEVB; Assert (static_cast&>(*this->fe) == static_cast&>(cell->get_fe()), typename FEVB::ExcFEDontMatch()); + check_cell_similarity(cell); + // set new cell. auto_ptr will take // care that old object gets // destroyed and also that this @@ -3420,6 +3424,8 @@ FEValues::reinit (const typename MGDoFHandler::cell_ // static_cast&>(cell->get_fe()), // typename FEValuesBase::ExcFEDontMatch()); + check_cell_similarity(cell); + // set new cell. auto_ptr will take // care that old object gets // destroyed and also that this @@ -3442,8 +3448,9 @@ FEValues::reinit (const typename MGDoFHandler::cell_ template void FEValues::reinit (const typename Triangulation::cell_iterator &cell) { - // no FE in this cell, so no check + // no FE in this cell, so no assertion // necessary here + check_cell_similarity(cell); // set new cell. auto_ptr will take // care that old object gets @@ -3462,11 +3469,59 @@ void FEValues::reinit (const typename Triangulation: +template +inline +void +FEValues::check_cell_similarity (const typename Triangulation::cell_iterator &cell) +{ + // case that there has not been any cell + // before + if (&*this->present_cell == 0) + { + cell_similarity = CellSimilarity::no_similarity; + return; + } + + const typename Triangulation::cell_iterator & present_cell = + *this->present_cell; + + // test for translation + { + // otherwise, go through the vertices and + // check... + bool is_translation = true; + const Point dist = cell->vertex(0) - present_cell->vertex(0); + for (unsigned int i=1; i::vertices_per_cell; ++i) + { + Point dist_new = cell->vertex(i) - present_cell->vertex(i) - dist; + if (dist_new.norm_square() > 1e-28) + { + is_translation = false; + break; + } + } + + cell_similarity = (is_translation == true + ? + CellSimilarity::translation + : + CellSimilarity::no_similarity); + return; + } + + // TODO: here, one could implement other + // checks for similarity, e.g. for + // children of a parallelogram. +} + + + template void FEValues::do_reinit () { this->get_mapping().fill_fe_values(*this->present_cell, quadrature, + this->cell_similarity, *this->mapping_data, this->quadrature_points, this->JxW_values, @@ -3478,6 +3533,7 @@ void FEValues::do_reinit () this->get_fe().fill_fe_values(this->get_mapping(), *this->present_cell, quadrature, + this->cell_similarity, *this->mapping_data, *this->fe_data, *this); @@ -3502,11 +3558,11 @@ FEValues::memory_consumption () const template FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, - const unsigned int dofs_per_cell, - const UpdateFlags, - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature& quadrature) + const unsigned int dofs_per_cell, + const UpdateFlags, + const Mapping &mapping, + const FiniteElement &fe, + const Quadrature& quadrature) : FEValuesBase (n_q_points, dofs_per_cell, diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index f908b2bb8b..ec021ded6c 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -79,7 +79,7 @@ MappingCartesian::update_each (const UpdateFlags in) const template typename Mapping::InternalDataBase * MappingCartesian::get_data (const UpdateFlags update_flags, - const Quadrature &q) const + const Quadrature &q) const { InternalData* data = new InternalData (q); @@ -95,7 +95,7 @@ MappingCartesian::get_data (const UpdateFlags update_flags, template typename Mapping::InternalDataBase * MappingCartesian::get_face_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const + const Quadrature& quadrature) const { InternalData* data = new InternalData (QProjector::project_to_all_faces(quadrature)); @@ -112,7 +112,7 @@ MappingCartesian::get_face_data (const UpdateFlags update_flags, template typename Mapping::InternalDataBase * MappingCartesian::get_subface_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const + const Quadrature &quadrature) const { InternalData* data = new InternalData (QProjector::project_to_all_subfaces(quadrature)); @@ -130,11 +130,12 @@ MappingCartesian::get_subface_data (const UpdateFlags update_flag template void MappingCartesian::compute_fill (const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - InternalData &data, - std::vector > &quadrature_points, - std::vector > &normal_vectors) const + const unsigned int face_no, + const unsigned int sub_no, + const enum CellSimilarity::Similarity cell_similarity, + InternalData &data, + std::vector > &quadrature_points, + std::vector > &normal_vectors) const { const UpdateFlags update_flags(data.current_update_flags()); @@ -299,6 +300,7 @@ void MappingCartesian:: fill_fe_values (const typename Triangulation::cell_iterator& cell, const Quadrature& q, + const enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase& mapping_data, std::vector >& quadrature_points, std::vector& JxW_values, @@ -311,11 +313,12 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, // data for this class. fails with // an exception if that is not // possible - InternalData &data = dynamic_cast (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); + InternalData &data = static_cast (mapping_data); std::vector > dummy; - compute_fill (cell, invalid_face_number, invalid_face_number, + compute_fill (cell, invalid_face_number, invalid_face_number, cell_similarity, data, quadrature_points, dummy); @@ -326,40 +329,44 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, // local lengths in each coordinate // direction if (data.current_update_flags() & (update_JxW_values | update_volume_elements)) - { - double J = data.length[0]; - for (unsigned int d=1;d(); - for (unsigned int j=0; j(); + for (unsigned int j=0; j(); + if (cell_similarity != CellSimilarity::translation) + 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(); + for (unsigned int j=0; j::fill_fe_face_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &data = dynamic_cast (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, + ExcInternalError()); + InternalData &data = static_cast (mapping_data); - compute_fill (cell, face_no, invalid_face_number, + compute_fill (cell, face_no, invalid_face_number, + CellSimilarity::no_similarity, data, quadrature_points, normal_vectors); @@ -471,9 +481,10 @@ MappingCartesian::fill_fe_subface_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &data = dynamic_cast (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); + InternalData &data = static_cast (mapping_data); - compute_fill (cell, face_no, sub_no, + compute_fill (cell, face_no, sub_no, CellSimilarity::no_similarity, data, quadrature_points, normal_vectors); @@ -535,7 +546,7 @@ MappingCartesian::transform ( AssertDimension (input.size(), output.size()); Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); - const InternalData &data = dynamic_cast (mapping_data); + const InternalData &data = static_cast (mapping_data); switch (mapping_type) { diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 658eb90f32..fc6c6336b5 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -198,7 +198,7 @@ MappingQ<1>::compute_shapes_virtual (const std::vector > &unit_points, template void MappingQ::compute_shapes_virtual (const std::vector > &unit_points, - typename MappingQ1::InternalData &data) const + typename MappingQ1::InternalData &data) const { const unsigned int n_points=unit_points.size(); std::vector values; @@ -250,7 +250,7 @@ MappingQ::compute_shapes_virtual (const std::vector > & template typename Mapping::InternalDataBase * MappingQ::get_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const + const Quadrature &quadrature) const { InternalData *data = new InternalData(n_shape_functions); this->compute_data (update_flags, quadrature, @@ -266,7 +266,7 @@ MappingQ::get_data (const UpdateFlags update_flags, template typename Mapping::InternalDataBase * MappingQ::get_face_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const + const Quadrature& quadrature) const { InternalData *data = new InternalData(n_shape_functions); const Quadrature q (QProjector::project_to_all_faces(quadrature)); @@ -284,7 +284,7 @@ MappingQ::get_face_data (const UpdateFlags update_flags, template typename Mapping::InternalDataBase * MappingQ::get_subface_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const + const Quadrature& quadrature) const { InternalData *data = new InternalData(n_shape_functions); const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); @@ -298,12 +298,16 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, } - + // Note that the CellSimilarity flag is + // modifyable, since MappingQ can need to + // recalculate data even when cells are + // similar. template void MappingQ::fill_fe_values ( const typename Triangulation::cell_iterator &cell, const Quadrature &q, + enum CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, std::vector > &quadrature_points, std::vector &JxW_values, @@ -316,7 +320,8 @@ MappingQ::fill_fe_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &data = dynamic_cast (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); + InternalData &data = static_cast (mapping_data); // check whether this cell needs // the full mapping or can be @@ -326,16 +331,30 @@ MappingQ::fill_fe_values ( data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use - // this or the other data object - // for the mapping + // depending on this result, use this or + // the other data object for the + // mapping. furthermore, we need to + // ensure that the flag indicating + // whether we can use some similarity has + // to be modified - for a general + // MappingQ, the data needs to be + // recomputed anyway since then the + // mapping changes the data. this needs + // to be known also for later operations, + // so modify the variable here. TODO: + // still need to check how this will work + // when the previous cell disabled this + // flag. typename MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) p_data=&data.mapping_q1_data; else - p_data=&data; + { + p_data=&data; + cell_similarity = CellSimilarity::no_similarity; + } - MappingQ1::fill_fe_values(cell, q, *p_data, + MappingQ1::fill_fe_values(cell, q, cell_similarity, *p_data, quadrature_points, JxW_values, jacobians, jacobian_grads, inverse_jacobians, cell_normal_vectors); @@ -359,7 +378,9 @@ MappingQ::fill_fe_face_values ( // data for this class. fails with // an exception if that is not // possible - InternalData &data = dynamic_cast (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, + ExcInternalError()); + InternalData &data = static_cast (mapping_data); // check whether this cell needs // the full mapping or can be @@ -419,7 +440,9 @@ MappingQ::fill_fe_subface_values (const typename Triangulation (mapping_data); + Assert (dynamic_cast (&mapping_data) != 0, + ExcInternalError()); + InternalData &data = static_cast (mapping_data); // check whether this cell needs // the full mapping or can be diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 6ac2276306..d7b3aa9f00 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -73,7 +73,7 @@ MappingQ1::MappingQ1 () template void MappingQ1::compute_shapes (const std::vector > &unit_points, - InternalData &data) const + InternalData &data) const { // choose either the function implemented // in this class, or whatever a virtual @@ -454,9 +454,9 @@ MappingQ1::update_each (const UpdateFlags in) const template void MappingQ1::compute_data (const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const { const unsigned int n_q_points = q.size(); @@ -504,9 +504,9 @@ MappingQ1::get_data (const UpdateFlags update_flags, template void MappingQ1::compute_face_data (const UpdateFlags update_flags, - const Quadrature& q, - const unsigned int n_original_q_points, - InternalData& data) const + const Quadrature& q, + const unsigned int n_original_q_points, + InternalData& data) const { compute_data (update_flags, q, n_original_q_points, data); @@ -572,7 +572,7 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, template typename Mapping::InternalDataBase * MappingQ1::get_face_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const + const Quadrature &quadrature) const { InternalData* data = new InternalData(n_shape_functions); compute_face_data (update_flags, @@ -606,33 +606,13 @@ void MappingQ1::compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int n_q_points, const DataSetDescriptor data_set, + const enum CellSimilarity::Similarity cell_similarity, InternalData &data, std::vector > &quadrature_points) const { const UpdateFlags update_flags(data.current_update_flags()); - if (update_flags & update_quadrature_points) - { - Assert (quadrature_points.size() == n_q_points, - ExcDimensionMismatch(quadrature_points.size(), n_q_points)); - std::fill(quadrature_points.begin(), quadrature_points.end(), - Point()); - } - - if (update_flags & update_contravariant_transformation) - { - Assert (data.contravariant.size() == n_q_points, - ExcDimensionMismatch(data.contravariant.size(), n_q_points)); - std::fill(data.contravariant.begin(), data.contravariant.end(), - Tensor<2,spacedim>()); - } - - if (update_flags & update_covariant_transformation) - { - Assert (data.covariant.size() == n_q_points, - ExcDimensionMismatch(data.covariant.size(), n_q_points)); - } // if necessary, recompute the // support points of the // transformation of this cell @@ -655,61 +635,89 @@ MappingQ1::compute_fill (const typename Triangulation()); + + for (unsigned int point=0; point()); + for (unsigned int point=0; point &data_derv = data.derivative(point+data_set, k); - const Tensor<1,spacedim> &supp_pts = data.mapping_support_points[k]; + const Tensor<1,dim> &data_derv = data.derivative(point+data_set, k); + const Tensor<1,spacedim> &supp_pts = data.mapping_support_points[k]; - for (unsigned int i=0; i contravariant_matrix(spacedim,dim); - FullMatrix covariant_matrix(dim,spacedim); + FullMatrix contravariant_matrix(spacedim,dim); + FullMatrix covariant_matrix(dim,spacedim); - for (unsigned int point=0; point void MappingQ1::fill_fe_values ( const typename Triangulation::cell_iterator &cell, - const Quadrature &q, - typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, - std::vector &JxW_values, - std::vector > &jacobians, - std::vector > &jacobian_grads, - std::vector > &inverse_jacobians, - std::vector > &cell_normal_vectors) const + const Quadrature &q, + const enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors) const { // ensure that the following cast is really correct: Assert (dynamic_cast(&mapping_data) != 0, @@ -746,14 +755,13 @@ MappingQ1::fill_fe_values ( const unsigned int n_q_points=q.size(); - compute_fill (cell, n_q_points, DataSetDescriptor::cell (), + compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity, data, quadrature_points); const UpdateFlags update_flags(data.current_update_flags()); const std::vector &weights=q.get_weights(); - // Multiply quadrature weights by // Jacobian determinants or by norm // of the cross product of the @@ -769,6 +777,7 @@ MappingQ1::fill_fe_values ( Assert( !(update_flags & update_cell_normal_vectors ) || (cell_normal_vectors.size() == n_q_points), ExcDimensionMismatch(cell_normal_vectors.size(), n_q_points)); + if (cell_similarity != CellSimilarity::translation) for (unsigned int point=0; point::fill_fe_values ( { Assert (jacobians.size() == n_q_points, ExcDimensionMismatch(jacobians.size(), n_q_points)); + if (cell_similarity != CellSimilarity::translation) for (unsigned int point=0; point::fill_fe_values ( { Assert (jacobian_grads.size() == n_q_points, ExcDimensionMismatch(jacobian_grads.size(), n_q_points)); - std::fill(jacobian_grads.begin(), - jacobian_grads.end(), - Tensor<3,spacedim>()); - for (unsigned int point=0; point()); + + for (unsigned int point=0; point::fill_fe_values ( { Assert (inverse_jacobians.size() == n_q_points, ExcDimensionMismatch(inverse_jacobians.size(), n_q_points)); - for (unsigned int point=0; point::compute_fill_face ( std::vector > &boundary_forms, std::vector > &normal_vectors) const { - compute_fill (cell, n_q_points, data_set, data, quadrature_points); + compute_fill (cell, n_q_points, data_set, CellSimilarity::no_similarity, + data, quadrature_points); const UpdateFlags update_flags(data.current_update_flags()); diff --git a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc index 7f4f1d5ac1..287a809e10 100644 --- a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc @@ -27,9 +27,9 @@ template MappingQ1Eulerian:: MappingQ1Eulerian (const EulerVectorType &euler_transform_vectors, const DoFHandler &shiftmap_dof_handler) - : - euler_transform_vectors(euler_transform_vectors), - shiftmap_dof_handler(&shiftmap_dof_handler) + : + euler_transform_vectors(euler_transform_vectors), + shiftmap_dof_handler(&shiftmap_dof_handler) {} @@ -114,6 +114,33 @@ MappingQ1Eulerian::clone () const } + +template +void +MappingQ1Eulerian::fill_fe_values ( + const typename Triangulation::cell_iterator &cell, + const Quadrature &q, + enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors) const +{ + // disable any previously detected + // similarity and hand on to the + // respective function of the base class. + cell_similarity = CellSimilarity::no_similarity; + MappingQ1::fill_fe_values (cell, q, cell_similarity, mapping_data, + quadrature_points, JxW_values, jacobians, + jacobian_grads, inverse_jacobians, + cell_normal_vectors); +} + + + // explicit instantiation template class MappingQ1Eulerian >; #ifdef DEAL_II_USE_PETSC diff --git a/deal.II/deal.II/source/fe/mapping_q_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc index f12842a2cf..cfa9ce1930 100644 --- a/deal.II/deal.II/source/fe/mapping_q_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc @@ -182,6 +182,32 @@ compute_mapping_support_points +template +void +MappingQEulerian::fill_fe_values ( + const typename Triangulation::cell_iterator &cell, + const Quadrature &q, + enum CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors) const +{ + // disable any previously detected + // similarity and hand on to the respective + // function of the base class. + cell_similarity = CellSimilarity::no_similarity; + MappingQ::fill_fe_values (cell, q, cell_similarity, mapping_data, + quadrature_points, JxW_values, jacobians, + jacobian_grads, inverse_jacobians, + cell_normal_vectors); +} + + + // explicit instantiation template class MappingQEulerian >; #ifdef DEAL_II_USE_PETSC -- 2.39.5