From 5f370a40a77f6b3cec288c7e5d39e96e3af4bb6b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 7 Mar 2001 19:10:19 +0000 Subject: [PATCH] Fix several bugs where by copy-paste some places were modified some not. Many annotations, many TODOs. Major code merge is needed in several places. git-svn-id: https://svn.dealii.org/trunk@4162 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_values.cc | 342 +++++++++++++++++++------ 1 file changed, 265 insertions(+), 77 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 40b810436d..eab5e8a4ca 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -31,13 +31,15 @@ using namespace std; #endif + +// TODO: replace this by non-global object static MappingQ1 mapping_q1; template void FEValuesData::initialize (const unsigned int n_quadrature_points, - unsigned int n_shapes, - const UpdateFlags flags) + const unsigned int n_shapes, + const UpdateFlags flags) { if (flags & update_values) shape_values.reinit(n_shapes, n_quadrature_points); @@ -51,7 +53,7 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, if (flags & update_second_derivatives) { - shape_2nd_derivatives.resize(n_shapes); + shape_2nd_derivatives.resize(n_shapes); for (unsigned int i=0;i::FEValuesBase (const unsigned int n_q_points, } + template FEValuesBase::~FEValuesBase () { @@ -113,6 +116,7 @@ FEValuesBase::~FEValuesBase () } + template template void FEValuesBase::get_function_values (const InputVector &fe_function, @@ -196,6 +200,7 @@ FEValuesBase::get_shape_grads () const }; + template const typename std::vector > > & FEValuesBase::get_shape_2nd_derivatives () const @@ -205,6 +210,7 @@ FEValuesBase::get_shape_2nd_derivatives () const }; + template const typename std::vector > & FEValuesBase::get_quadrature_points () const @@ -214,6 +220,7 @@ FEValuesBase::get_quadrature_points () const }; + template const std::vector & FEValuesBase::get_JxW_values () const @@ -223,6 +230,7 @@ FEValuesBase::get_JxW_values () const } + template template void FEValuesBase::get_function_grads (const InputVector &fe_function, @@ -296,6 +304,7 @@ void FEValuesBase::get_function_grads (const InputVector &fe_ }; + template template void FEValuesBase::get_function_2nd_derivatives (const InputVector &fe_function, @@ -330,6 +339,7 @@ void FEValuesBase::get_function_2nd_derivatives (const InputVector &fe }; + template template void @@ -398,6 +408,7 @@ template unsigned int FEValuesBase::memory_consumption () const { +//TODO: check whether this is up to date return (MemoryConsumption::memory_consumption (shape_values) + MemoryConsumption::memory_consumption (shape_gradients) + MemoryConsumption::memory_consumption (shape_2nd_derivatives) + @@ -411,6 +422,8 @@ FEValuesBase::memory_consumption () const /*------------------------------- FEValues -------------------------------*/ + +// TODO: unify the two constructor bodies by calling a single function that does all the work template FEValues::FEValues (const Mapping &mapping, const FiniteElement &fe, @@ -425,56 +438,115 @@ FEValues::FEValues (const Mapping &mapping, fe), quadrature (q) { + // you can't compute normal vectors + // on cells, only on faces Assert ((update_flags & update_normal_vectors) == false, typename FEValuesBase::ExcInvalidUpdateFlag()); - UpdateFlags flags = mapping.update_once (update_flags); - flags |= mapping.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping.update_once (update_flags) | + mapping.update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit mapping_data = mapping.get_data(flags, quadrature); fe_data = fe.get_data(flags, mapping, quadrature); - UpdateFlags allflags = mapping_data->update_once | mapping_data->update_each; - allflags |= fe_data->update_once | fe_data->update_each; - + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, allflags); }; + template FEValues::FEValues (const FiniteElement &fe, const Quadrature &q, const UpdateFlags update_flags) - : - FEValuesBase (q.n_quadrature_points, - fe.dofs_per_cell, - 1, - update_flags, - mapping_q1, - fe), - quadrature (q) + : + FEValuesBase (q.n_quadrature_points, + fe.dofs_per_cell, + 1, + update_flags, + mapping_q1, + fe), + quadrature (q) { + // you can't compute normal vectors + // on cells, only on faces Assert ((update_flags & update_normal_vectors) == false, FEValuesBase::ExcInvalidUpdateFlag()); - UpdateFlags flags = mapping_q1.update_once (update_flags); - flags |= mapping_q1.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - - mapping_data = mapping_q1.get_data(flags, quadrature); - fe_data = fe.get_data(flags, mapping_q1, quadrature); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping->update_once (update_flags) | + mapping->update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit + mapping_data = mapping->get_data(flags, quadrature); + fe_data = fe.get_data(flags, *mapping, quadrature); + + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); } + template void FEValues::reinit (const typename DoFHandler::cell_iterator &cell) { @@ -508,6 +580,7 @@ template unsigned int FEValues::memory_consumption () const { +//TODO: check whether this is up to date return FEValuesBase::memory_consumption (); }; @@ -523,17 +596,18 @@ FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, const Mapping &mapping, const FiniteElement &fe, const Quadrature& quadrature) - : - FEValuesBase (n_q_points, - dofs_per_cell, - n_faces_or_subfaces, - update_flags, - mapping, - fe), - quadrature(quadrature) + : + FEValuesBase (n_q_points, + dofs_per_cell, + n_faces_or_subfaces, + update_flags, + mapping, + fe), + quadrature(quadrature) {}; + template const std::vector > & FEFaceValuesBase::get_normal_vectors () const @@ -544,6 +618,7 @@ FEFaceValuesBase::get_normal_vectors () const }; + template const std::vector > & FEFaceValuesBase::get_boundary_forms () const @@ -557,6 +632,7 @@ FEFaceValuesBase::get_boundary_forms () const /*------------------------------- FEFaceValues -------------------------------*/ +// TODO: unify the two constructor bodies by calling a single function that does all the work template FEFaceValues::FEFaceValues (const Mapping &mapping, const FiniteElement &fe, @@ -570,17 +646,42 @@ FEFaceValues::FEFaceValues (const Mapping &mapping, mapping, fe, quadrature) { - UpdateFlags flags = mapping.update_once (update_flags); - flags |= mapping.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping.update_once (update_flags) | + mapping.update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit mapping_data = mapping.get_face_data(flags, quadrature); fe_data = fe.get_face_data(flags, mapping, quadrature); - UpdateFlags allflags = mapping_data->update_once | mapping_data->update_each; - allflags |= fe_data->update_once | fe_data->update_each; - + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, allflags); @@ -600,24 +701,51 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, mapping_q1, fe, quadrature) { - UpdateFlags flags = mapping_q1.update_once (update_flags); - flags |= mapping_q1.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - - mapping_data = mapping_q1.get_face_data(flags, quadrature); - fe_data = fe.get_face_data(flags, mapping_q1, quadrature); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping->update_once (update_flags) | + mapping->update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit + mapping_data = mapping->get_face_data(flags, quadrature); + fe_data = fe.get_face_data(flags, *mapping, quadrature); + + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); }; template -void FEFaceValues::reinit (const typename - DoFHandler::cell_iterator &cell, +void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &cell, const unsigned int face_no) { // assert that the finite elements @@ -628,8 +756,8 @@ void FEFaceValues::reinit (const typename static_cast&>(cell->get_dof_handler().get_fe()), typename FEValuesBase::ExcFEDontMatch()); - present_cell = cell; - present_face = cell->face(face_no); + present_cell = cell; + present_face = cell->face(face_no); get_mapping().fill_fe_face_values(cell, face_no, quadrature, @@ -651,6 +779,8 @@ void FEFaceValues::reinit (const typename /*------------------------------- FESubFaceValues -------------------------------*/ +// TODO: unify the two constructor bodies by calling a single function that does all the work +// TODO: FESubfaceValues constructors are exactly the same as FEFaceValues constructors with exception of calling fe_sub_face_data instead of fe_face_data. Merge! template FESubfaceValues::FESubfaceValues (const Mapping &mapping, const FiniteElement &fe, @@ -659,22 +789,51 @@ FESubfaceValues::FESubfaceValues (const Mapping &mapping, : FEFaceValuesBase (quadrature.n_quadrature_points, fe.dofs_per_cell, - GeometryInfo::faces_per_cell * GeometryInfo::subfaces_per_face, + GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face, update_flags, mapping, fe, quadrature) { - UpdateFlags flags = mapping.update_once (update_flags); - flags |= mapping.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - - mapping_data = mapping.get_subface_data(flags, quadrature); - fe_data = fe.get_subface_data(flags, mapping, quadrature); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping.update_once (update_flags) | + mapping.update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit + mapping_data = mapping.get_sub_face_data(flags, quadrature); + fe_data = fe.get_sub_face_data(flags, mapping, quadrature); + + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); } @@ -686,22 +845,51 @@ FESubfaceValues::FESubfaceValues (const FiniteElement &fe, : FEFaceValuesBase (quadrature.n_quadrature_points, fe.dofs_per_cell, - GeometryInfo::faces_per_cell * GeometryInfo::subfaces_per_face, + GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face, update_flags, mapping_q1, fe, quadrature) { - UpdateFlags flags = mapping_q1.update_once (update_flags); - flags |= mapping_q1.update_each (update_flags); - flags |= fe.update_once (update_flags); - flags |= fe.update_each (update_flags); - - mapping_data = mapping_q1.get_subface_data(flags, quadrature); - fe_data = fe.get_subface_data(flags, mapping_q1, quadrature); - + // first find out which objects + // need to be recomputed on each + // cell we visit. this we have to + // ask the finite element + const UpdateFlags flags = UpdateFlags(mapping->update_once (update_flags) | + mapping->update_each (update_flags) | + fe.update_once (update_flags) | + fe.update_each (update_flags) ); + + // then get objects into which the + // FE and the Mapping can store + // intermediate data used across + // calls to reinit + mapping_data = mapping->get_sub_face_data(flags, quadrature); + fe_data = fe.get_sub_face_data(flags, *mapping, quadrature); + + // maybe some of the flags passed + // to the get_data functions + // resulted in other flags to be + // set, so collect them again. + // + // example: a FE might need to have + // the covariant transform of the + // mapping to compute its shape + // values +//TODO: why can't the flags be obtained in one cycle?, e.g. by calling +// flags = fe.update_once(update_flags) +// flags = mapping.update_once(flags) +// if the mapping should depend on the fe to compute some things, then +// we would need to iterate, and the one cycle here isn't sufficient!? + const UpdateFlags allflags = UpdateFlags(mapping_data->update_once | + mapping_data->update_each | + fe_data->update_once | + fe_data->update_each ); + + // set up objects within this class FEValuesData::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); } -- 2.39.5