From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Wed, 7 Mar 2001 19:10:19 +0000 (+0000) Subject: Fix several bugs where by copy-paste some places were modified some X-Git-Tag: v8.0.0~19609 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5f370a40a77f6b3cec288c7e5d39e96e3af4bb6b;p=dealii.git 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 --- 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<deal_II_dimension> mapping_q1; template <int dim> void FEValuesData<dim>::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<dim>::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<n_shapes;++i) shape_2nd_derivatives[i].resize(n_quadrature_points); } @@ -94,6 +96,7 @@ FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points, } + template <int dim> FEValuesBase<dim>::~FEValuesBase () { @@ -113,6 +116,7 @@ FEValuesBase<dim>::~FEValuesBase () } + template <int dim> template <class InputVector, typename number> void FEValuesBase<dim>::get_function_values (const InputVector &fe_function, @@ -196,6 +200,7 @@ FEValuesBase<dim>::get_shape_grads () const }; + template <int dim> const typename std::vector<typename std::vector<Tensor<2,dim> > > & FEValuesBase<dim>::get_shape_2nd_derivatives () const @@ -205,6 +210,7 @@ FEValuesBase<dim>::get_shape_2nd_derivatives () const }; + template <int dim> const typename std::vector<Point<dim> > & FEValuesBase<dim>::get_quadrature_points () const @@ -214,6 +220,7 @@ FEValuesBase<dim>::get_quadrature_points () const }; + template <int dim> const std::vector<double> & FEValuesBase<dim>::get_JxW_values () const @@ -223,6 +230,7 @@ FEValuesBase<dim>::get_JxW_values () const } + template <int dim> template <class InputVector> void FEValuesBase<dim>::get_function_grads (const InputVector &fe_function, @@ -296,6 +304,7 @@ void FEValuesBase<dim>::get_function_grads (const InputVector &fe_ }; + template <int dim> template <class InputVector> void FEValuesBase<dim>::get_function_2nd_derivatives (const InputVector &fe_function, @@ -330,6 +339,7 @@ void FEValuesBase<dim>::get_function_2nd_derivatives (const InputVector &fe }; + template <int dim> template <class InputVector> void @@ -398,6 +408,7 @@ template <int dim> unsigned int FEValuesBase<dim>::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<dim>::memory_consumption () const /*------------------------------- FEValues -------------------------------*/ + +// TODO: unify the two constructor bodies by calling a single function that does all the work template <int dim> FEValues<dim>::FEValues (const Mapping<dim> &mapping, const FiniteElement<dim> &fe, @@ -425,56 +438,115 @@ FEValues<dim>::FEValues (const Mapping<dim> &mapping, fe), quadrature (q) { + // you can't compute normal vectors + // on cells, only on faces Assert ((update_flags & update_normal_vectors) == false, typename FEValuesBase<dim>::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<dim>::initialize(n_quadrature_points, dofs_per_cell, allflags); }; + template <int dim> FEValues<dim>::FEValues (const FiniteElement<dim> &fe, const Quadrature<dim> &q, const UpdateFlags update_flags) - : - FEValuesBase<dim> (q.n_quadrature_points, - fe.dofs_per_cell, - 1, - update_flags, - mapping_q1, - fe), - quadrature (q) + : + FEValuesBase<dim> (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<dim>::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<dim>::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); } + template <int dim> void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell) { @@ -508,6 +580,7 @@ template <int dim> unsigned int FEValues<dim>::memory_consumption () const { +//TODO: check whether this is up to date return FEValuesBase<dim>::memory_consumption (); }; @@ -523,17 +596,18 @@ FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points, const Mapping<dim> &mapping, const FiniteElement<dim> &fe, const Quadrature<dim-1>& quadrature) - : - FEValuesBase<dim> (n_q_points, - dofs_per_cell, - n_faces_or_subfaces, - update_flags, - mapping, - fe), - quadrature(quadrature) + : + FEValuesBase<dim> (n_q_points, + dofs_per_cell, + n_faces_or_subfaces, + update_flags, + mapping, + fe), + quadrature(quadrature) {}; + template <int dim> const std::vector<Point<dim> > & FEFaceValuesBase<dim>::get_normal_vectors () const @@ -544,6 +618,7 @@ FEFaceValuesBase<dim>::get_normal_vectors () const }; + template <int dim> const std::vector<Tensor<1,dim> > & FEFaceValuesBase<dim>::get_boundary_forms () const @@ -557,6 +632,7 @@ FEFaceValuesBase<dim>::get_boundary_forms () const /*------------------------------- FEFaceValues -------------------------------*/ +// TODO: unify the two constructor bodies by calling a single function that does all the work template <int dim> FEFaceValues<dim>::FEFaceValues (const Mapping<dim> &mapping, const FiniteElement<dim> &fe, @@ -570,17 +646,42 @@ FEFaceValues<dim>::FEFaceValues (const Mapping<dim> &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<dim>::initialize(n_quadrature_points, dofs_per_cell, allflags); @@ -600,24 +701,51 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &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<dim>::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); }; template <int dim> -void FEFaceValues<dim>::reinit (const typename - DoFHandler<dim>::cell_iterator &cell, +void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell, const unsigned int face_no) { // assert that the finite elements @@ -628,8 +756,8 @@ void FEFaceValues<dim>::reinit (const typename static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()), typename FEValuesBase<dim>::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<dim>::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 <int dim> FESubfaceValues<dim>::FESubfaceValues (const Mapping<dim> &mapping, const FiniteElement<dim> &fe, @@ -659,22 +789,51 @@ FESubfaceValues<dim>::FESubfaceValues (const Mapping<dim> &mapping, : FEFaceValuesBase<dim> (quadrature.n_quadrature_points, fe.dofs_per_cell, - GeometryInfo<dim>::faces_per_cell * GeometryInfo<dim>::subfaces_per_face, + GeometryInfo<dim>::faces_per_cell * + GeometryInfo<dim>::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<dim>::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); } @@ -686,22 +845,51 @@ FESubfaceValues<dim>::FESubfaceValues (const FiniteElement<dim> &fe, : FEFaceValuesBase<dim> (quadrature.n_quadrature_points, fe.dofs_per_cell, - GeometryInfo<dim>::faces_per_cell * GeometryInfo<dim>::subfaces_per_face, + GeometryInfo<dim>::faces_per_cell * + GeometryInfo<dim>::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<dim>::initialize(n_quadrature_points, dofs_per_cell, - flags); + allflags); }