From 14f71faea435c55442887baa594ec2cc1c1bbc50 Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 9 Aug 2011 09:38:41 +0000 Subject: [PATCH] fix a bug in FESystem:\ the functions get_data, get_face_data and get_subface_data did not honor if a base element required a separate implementation git-svn-id: https://svn.dealii.org/trunk@24036 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/fe_system.h | 33 ++++- deal.II/source/fe/fe_system.cc | 197 +++++++++++++++++++++---- 2 files changed, 193 insertions(+), 37 deletions(-) diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index f560ed9d43..44dfd2e48d 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -678,16 +678,21 @@ class FESystem : public FiniteElement */ virtual FiniteElement * clone() const; - /** - * Prepare internal data - * structures and fill in values - * independent of the cell. - */ virtual typename Mapping::InternalDataBase* get_data (const UpdateFlags update_flags, const Mapping &mapping, const Quadrature &quadrature) const ; + virtual typename Mapping::InternalDataBase* + get_face_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const ; + + virtual typename Mapping::InternalDataBase* + get_subface_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const ; + /** * Implementation of the same * function in @@ -1262,6 +1267,22 @@ class FESystem : public FiniteElement }; }; +template <> +typename Mapping<1,1>::InternalDataBase* +FESystem<1,1>::get_face_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const; + +template <> +typename Mapping<1,1>::InternalDataBase* +FESystem<1,1>::get_subface_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const; + +template <> +typename Mapping<1,2>::InternalDataBase* +FESystem<1,2>::get_face_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const; + +template <> +typename Mapping<1,2>::InternalDataBase* +FESystem<1,2>::get_subface_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const; + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 6ad4d931c4..d4fdd7881c 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -685,7 +685,6 @@ FESystem::get_data (const UpdateFlags flags_, const Quadrature &quadrature) const { UpdateFlags flags = flags_; - InternalData* data = new InternalData(n_base_elements()); data->update_once = update_once (flags); @@ -693,7 +692,6 @@ FESystem::get_data (const UpdateFlags 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 @@ -708,7 +706,6 @@ FESystem::get_data (const UpdateFlags flags_, 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_no::get_data (const UpdateFlags flags_, return data; } +// The following function is a clone of get_data, with the exception +// that get_face_data of the base elements is called. + +template +typename Mapping::InternalDataBase * +FESystem::get_face_data ( + const UpdateFlags flags_, + const Mapping &mapping, + const Quadrature &quadrature) const +{ + UpdateFlags flags = flags_; + InternalData* data = new InternalData(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; + data->compute_hessians = flags & update_hessians; + 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_no::InternalDataBase *base_fe_data_base = + base_element(base_no).get_face_data(sub_flags, mapping, quadrature); + + typename FiniteElement::InternalDataBase *base_fe_data = + dynamic_cast::InternalDataBase *> + (base_fe_data_base); + Assert (base_fe_data != 0, ExcInternalError()); + + 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()); + + FEValuesData *base_data = new FEValuesData(); + data->set_fe_values_data(base_no, base_data); + } + data->update_flags = data->update_once | + data->update_each; + return data; +} + + + +// The following function is a clone of get_data, with the exception +// that get_subface_data of the base elements is called. + +template +typename Mapping::InternalDataBase * +FESystem::get_subface_data ( + const UpdateFlags flags_, + const Mapping &mapping, + const Quadrature &quadrature) const +{ + UpdateFlags flags = flags_; + InternalData* data = new InternalData(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; + data->compute_hessians = flags & update_hessians; + 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_no::InternalDataBase *base_fe_data_base = + base_element(base_no).get_subface_data(sub_flags, mapping, quadrature); + + typename FiniteElement::InternalDataBase *base_fe_data = + dynamic_cast::InternalDataBase *> + (base_fe_data_base); + Assert (base_fe_data != 0, ExcInternalError()); + + 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()); + + FEValuesData *base_data = new FEValuesData(); + data->set_fe_values_data(base_no, base_data); + } + data->update_flags = data->update_once | + data->update_each; + return data; +} + + +template <> +typename Mapping<1,1>::InternalDataBase * +FESystem<1,1>::get_face_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + +template <> +typename Mapping<1,1>::InternalDataBase * +FESystem<1,1>::get_subface_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + +template <> +typename Mapping<1,2>::InternalDataBase * +FESystem<1,2>::get_face_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + +template <> +typename Mapping<1,2>::InternalDataBase * +FESystem<1,2>::get_subface_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} template void -FESystem:: -fill_fe_values (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fe_data, - FEValuesData &data, - CellSimilarity::Similarity &cell_similarity) const +FESystem::fill_fe_values ( + const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fe_data, + FEValuesData &data, + CellSimilarity::Similarity &cell_similarity) const { compute_fill(mapping, cell, invalid_face_number, invalid_face_number, quadrature, cell_similarity, mapping_data, fe_data, data); @@ -803,15 +937,15 @@ fill_fe_face_values (const Mapping &mapping, template void -FESystem:: -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, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const +FESystem::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, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fe_data, + FEValuesData &data) const { compute_fill (mapping, cell, face_no, sub_no, quadrature, CellSimilarity::none, mapping_data, fe_data, data); @@ -822,16 +956,16 @@ fill_fe_subface_values (const Mapping &mappin template template void -FESystem:: -compute_fill (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, - CellSimilarity::Similarity cell_similarity, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const +FESystem::compute_fill ( + const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { const unsigned int n_q_points = quadrature.size(); @@ -1019,6 +1153,7 @@ compute_fill (const Mapping &mapping, // those shape functions // that belong to a given // base element +//TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style const UpdateFlags base_flags(dim_1==dim ? base_fe_data.current_update_flags() : base_fe_data.update_flags); @@ -1072,7 +1207,7 @@ compute_fill (const Mapping &mapping, for (unsigned int q=0; q