From: Wolfgang Bangerth Date: Wed, 7 Mar 2001 15:08:38 +0000 (+0000) Subject: Further clean-ups. X-Git-Tag: v8.0.0~19623 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73664b2d5db9f2a4da2a4e90209a4feddf8cc05d;p=dealii.git Further clean-ups. git-svn-id: https://svn.dealii.org/trunk@4148 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index cf32082d61..28e4038147 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -48,6 +48,7 @@ template class MatrixCreator; template class FiniteElement : public FiniteElementBase { + private: /** * Copy constructor prohibited. */ @@ -67,21 +68,6 @@ class FiniteElement : public FiniteElementBase */ virtual ~FiniteElement (); - /** - * Compute flags for initial - * update only. - * @see{FEValuesBase} - */ - virtual UpdateFlags update_once (UpdateFlags flags) const = 0; - - /** - * Compute flags for update on - * each cell. - * @see{FEValuesBase} - */ - virtual UpdateFlags update_each (UpdateFlags flags) const = 0; - - /** * Return the support points of the * trial functions on the unit cell. @@ -116,16 +102,73 @@ class FiniteElement : public FiniteElementBase * will be zero. This is the standard behavior, * if the function is not overloaded. */ - virtual void get_unit_face_support_points (std::vector > &) const; + virtual void get_unit_face_support_points (std::vector > &points) const; - friend class FEValues; - friend class FEFaceValues; - friend class FESubfaceValues; - friend class FESystem; - friend class MatrixCreator; - friend class VectorTools; + /** + * Number of base elements in a mixed + * discretization. This function returns + * 1 for primitive elements. + */ + virtual unsigned int n_base_elements () const; + + /** + * Access to base element + * objects. By default, + * #base_element(0)# is #this#. + * This function is overloaded by + * system elements to allow + * access to the different + * components of mixed + * discretizations. + */ + virtual const FiniteElement & base_element (const unsigned int index) const; + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + * + * This function is made virtual, + * since finite element objects + * are usually accessed through + * pointers to their base class, + * rather than the class itself. + */ + virtual unsigned int memory_consumption () const; + + /** + * Exception + */ + DeclException0 (ExcBoundaryFaceUsed); + /** + * Exception + */ + DeclException0 (ExcJacobiDeterminantHasWrongSign); + /** + * Exception + */ + DeclException1 (ExcComputationNotUseful, + int, + << "The computation you required from this function is not " + << "feasible or not probable in the present dimension (" + << arg1 << ") because it would be prohibitively expensive."); + protected: + + /** + * Compute flags for initial + * update only. + * @see{FEValuesBase} + */ + virtual UpdateFlags update_once (UpdateFlags flags) const = 0; + + /** + * Compute flags for update on + * each cell. + * @see{FEValuesBase} + */ + virtual UpdateFlags update_each (UpdateFlags flags) const = 0; + /** * @p{clone} function instead of * a copy constructor. @@ -228,58 +271,14 @@ class FiniteElement : public FiniteElementBase typename Mapping::InternalDataBase &fe_internal, FEValuesData &data) const = 0; - - public: - - /** - * Number of base elements in a mixed - * discretization. This function returns - * 1 for simple elements. - */ - virtual unsigned int n_base_elements () const; - - /** - * Access to base element - * objects. By default, - * #base_element(0)# is #this#. - * This function is overloaded by - * system elements to allow - * access to the different - * components of mixed - * discretizations. - */ - virtual const FiniteElement& base_element (const unsigned int index) const; - - /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. - * - * This function is made virtual, - * since finite element objects - * are usually accessed through - * pointers to their base class, - * rather than the class itself. - */ - virtual unsigned int memory_consumption () const; - /** - * Exception + * Declare some other classes as + * friends of this class. */ - DeclException0 (ExcBoundaryFaceUsed); - /** - * Exception - */ - DeclException0 (ExcJacobiDeterminantHasWrongSign); - /** - * Exception - */ - DeclException1 (ExcComputationNotUseful, - int, - << "The computation you required from this function is not " - << "feasible or not probable in the present dimension (" - << arg1 << ") because it would be prohibitively expensive."); - + friend class FEValues; + friend class FEFaceValues; + friend class FESubfaceValues; + friend class FESystem; }; diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 1b0ce3ea8f..1eb6f2edda 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -335,12 +335,14 @@ class FiniteElementBase : public Subscriptor, { public: /** - * Initialize @ref{FEValues} - * pointers for given element. + * Initialize some pointers + * used in the computation of + * second derivatives by finite + * differencing of gradients. */ - void initialize (const FiniteElement* element, - const Mapping& mapping, - const Quadrature& quadrature); + void initialize_2nd (const FiniteElement *element, + const Mapping &mapping, + const Quadrature &quadrature); /** * Destructor. Needed to avoid @@ -377,14 +379,16 @@ class FiniteElementBase : public Subscriptor, const std::vector &restriction_is_additive_flags); /** - * Compute second differences. + * Compute second derivatives by + * finite differences of + * gradients. */ - void compute_2nd (const Mapping &mapping, - const DoFHandler::cell_iterator &cell, - const unsigned int offset, + void compute_2nd (const Mapping &mapping, + const DoFHandler::cell_iterator &cell, + const unsigned int offset, typename Mapping::InternalDataBase &mapping_internal, - InternalDataBase& fe_internal, - FEValuesData& data) const; + InternalDataBase &fe_internal, + FEValuesData &data) const; /** * Projection from a fine grid diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 087ea67dff..295ebc4962 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -466,10 +466,11 @@ class FESystem : public FiniteElement ~InternalData(); /** - * Flag for computation of - * second derivatives. + * Flag indicating whether + * second derivatives shall + * be computed. */ - bool second_flag; + bool compute_second_derivatives; /** * Gives write-access to the diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 7718537308..01eda10290 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -34,35 +34,51 @@ using namespace std; template void FiniteElementBase:: -InternalDataBase::initialize (const FiniteElement* element, - const Mapping& mapping, - const Quadrature& quadrature) +InternalDataBase::initialize_2nd (const FiniteElement *element, + const Mapping &mapping, + const Quadrature &quadrature) { - // We compute difference - // quotients of gradients - UpdateFlags diff_flags = update_gradients; - - // We will need shifted - // quadrature formulae + // if we shall compute second + // derivatives, then we do so by + // finite differencing the + // gradients. that we do by + // evaluating the gradients of + // shape values at points shifted + // star-like a little in each + // coordinate direction around each + // quadrature point. + // + // therefore generate 2*dim (the + // number of evaluation points) + // FEValues objects with slightly + // shifted positions std::vector > diff_points (quadrature.n_quadrature_points); - std::vector diff_weights (quadrature.n_quadrature_points, 0); + std::vector diff_weights (quadrature.n_quadrature_points, 0); - // The star has 2dim points differences.resize(2*dim); - for (unsigned int d=0;d shift; +//TODO: unify the places where the finite differencing step length is used shift (d) = 1.e-6; - for (unsigned int i=0;i plus_quad (diff_points, diff_weights); - differences[d] = - new FEValues (mapping, *element, plus_quad, diff_flags); - for (unsigned int i=0;i plus_quad (diff_points, diff_weights); + differences[d] = new FEValues (mapping, *element, + plus_quad, update_gradients); + + // now same in minus-direction + for (unsigned int i=0; i minus_quad (diff_points, diff_weights); - differences[d+dim] = - new FEValues (mapping, *element, minus_quad, diff_flags); + const Quadrature minus_quad (diff_points, diff_weights); + differences[d+dim] = new FEValues (mapping, *element, + minus_quad, update_gradients); } } @@ -72,9 +88,15 @@ InternalDataBase::initialize (const FiniteElement* element, template FiniteElementBase::InternalDataBase::~InternalDataBase () { - for (unsigned int i=0;i::InternalDataBase::~InternalDataBase () template FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags) : - FiniteElementData (fe_data), - system_to_component_table(dofs_per_cell), - face_system_to_component_table(dofs_per_face), - component_to_system_table(components, std::vector(dofs_per_cell)), - face_component_to_system_table(components, std::vector(dofs_per_face)), - component_to_base_table(components), - restriction_is_additive_flags(restriction_is_additive_flags) + const std::vector &restriction_is_additive_flags) + : + FiniteElementData (fe_data), + system_to_component_table(dofs_per_cell), + face_system_to_component_table(dofs_per_face), + component_to_system_table(components, std::vector(dofs_per_cell)), + face_component_to_system_table(components, std::vector(dofs_per_face)), + component_to_base_table(components), + restriction_is_additive_flags(restriction_is_additive_flags) { Assert(restriction_is_additive_flags.size()==fe_data.components, ExcDimensionMismatch(restriction_is_additive_flags.size(),fe_data.components)); @@ -100,13 +123,16 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data prolongation[i].reinit (dofs_per_cell, dofs_per_cell); }; + // first set sizes of some + // matrices. they will be filled by + // derived classes, or re-set to + // zero size of not defined switch (dim) { case 1: Assert ((interface_constraints.m() == 0) && (interface_constraints.n() == 0), ExcInternalError()); - break; case 2: @@ -142,6 +168,7 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data }; + template const FullMatrix & FiniteElementBase::restrict (const unsigned int child) const @@ -153,6 +180,7 @@ FiniteElementBase::restrict (const unsigned int child) const }; + template const FullMatrix & FiniteElementBase::prolongate (const unsigned int child) const @@ -164,6 +192,7 @@ FiniteElementBase::prolongate (const unsigned int child) const }; + template const FullMatrix & FiniteElementBase::constraints () const @@ -180,6 +209,7 @@ FiniteElementBase::constraints () const }; + template bool FiniteElementBase::operator == (const FiniteElementBase &f) const { @@ -211,52 +241,76 @@ FiniteElementBase::memory_consumption () const }; + template void FiniteElementBase:: -compute_2nd (const Mapping &mapping, +compute_2nd (const Mapping &mapping, const DoFHandler::cell_iterator &cell, - const unsigned int offset, - Mapping::InternalDataBase &mapping_internal, - InternalDataBase& fe_internal, - FEValuesData& data) const + const unsigned int offset, + Mapping::InternalDataBase &mapping_internal, + InternalDataBase &fe_internal, + FEValuesData &data) const { // Number of quadrature points - const unsigned int n = data.shape_2nd_derivatives[0].size(); - - for (unsigned int d=0;dreinit(cell); fe_internal.differences[d+dim]->reinit(cell); } - std::vector > > diff_quot (dim, std::vector >(n)); - std::vector > diff_quot2 (n); - // Loop over shape functions + // collection of difference + // quotients of gradients in each + // direction (first index) and at + // all q-points (second index) + std::vector > > diff_quot (dim, std::vector >(n_q_points)); + std::vector > diff_quot2 (n_q_points); + + // for all shape functions at all + // quadrature points and difference + // quotients in all directions: for (unsigned int shape=0; shape& right - = fe_internal.differences[d1]->shape_grad(shape, k); + = fe_internal.differences[d1]->shape_grad(shape, q); const Tensor<1,dim>& left - = fe_internal.differences[d1+dim]->shape_grad(shape, k); - for (unsigned int d=0;dshape_grad(shape, q); + +//TODO: unify the places where the finite differencing step length is used + // compute the second + // derivative from a + // symmetric difference + // approximation + for (unsigned int d=0; d void FiniteElement::get_unit_support_points (std::vector > &points) const { + // default implementation points.resize(0); } + template void FiniteElement::get_unit_face_support_points (std::vector > &points) const { + // default implementation points.resize(0); } @@ -298,8 +355,8 @@ FiniteElement::get_unit_face_support_points (std::vector > &po template Mapping::InternalDataBase* -FiniteElement::get_face_data (const UpdateFlags flags, - const Mapping& mapping, +FiniteElement::get_face_data (const UpdateFlags flags, + const Mapping &mapping, const Quadrature &quadrature) const { QProjector q(quadrature, false); @@ -307,25 +364,30 @@ FiniteElement::get_face_data (const UpdateFlags flags, } + template Mapping::InternalDataBase* -FiniteElement::get_subface_data (const UpdateFlags flags, - const Mapping& mapping, +FiniteElement::get_subface_data (const UpdateFlags flags, + const Mapping &mapping, const Quadrature &quadrature) const { QProjector q(quadrature, true); return get_data (flags, mapping, q); - } + template unsigned int FiniteElement::n_base_elements() const { + // default implementation return 1; } + + + template unsigned int FiniteElement::memory_consumption () const diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 9e9c75693b..6e8055717e 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -394,10 +394,12 @@ FE_DGQ::get_data (const UpdateFlags update_flags, std::vector >(quadrature.n_quadrature_points)); } + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that if (flags & update_second_derivatives) - { - data->initialize (this, mapping, quadrature); - } + data->initialize_2nd (this, mapping, quadrature); if (flags & (update_values | update_gradients)) diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 29cb80fd3a..0bc632c5ba 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -989,8 +989,12 @@ FE_Q::get_data (const UpdateFlags update_flags, std::vector >(n_q_points)); } + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that if (flags & update_second_derivatives) - data->initialize (this, mapping, quadrature); + data->initialize_2nd (this, mapping, quadrature); if (flags & (update_values | update_gradients)) for (unsigned int i=0; i::update_each (UpdateFlags flags) const template Mapping::InternalDataBase* -FESystem::get_data (UpdateFlags flags, - const Mapping& mapping, - const Quadrature& quadrature) const +FESystem::get_data ( UpdateFlags flags, + const Mapping &mapping, + const Quadrature &quadrature) const { InternalData* data = new InternalData(n_base_elements()); - data->second_flag = flags & update_second_derivatives; - - // Make sure that this object - // computes 2nd derivatives itself - if (data->second_flag) + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that + data->compute_second_derivatives = flags & update_second_derivatives; + if (data->compute_second_derivatives) { + // delete + // update_second_derivatives + // from flags list flags = UpdateFlags (flags ^ update_second_derivatives); - data->initialize (this, mapping, quadrature); + data->initialize_2nd (this, mapping, quadrature); } @@ -177,13 +181,17 @@ FESystem::get_data (UpdateFlags flags, { typename Mapping::InternalDataBase *base_fe_data_base = base_element(base_no).get_data(flags, mapping, quadrature); + FiniteElementBase::InternalDataBase *base_fe_data = dynamic_cast::InternalDataBase *> (base_fe_data_base); data->set_fe_data(base_no, base_fe_data); - data->update_once|=base_fe_data->update_once; - data->update_each|=base_fe_data->update_each; + + // collect requirements of base + // elements + data->update_once |= base_fe_data->update_once; + data->update_each |= base_fe_data->update_each; // The FEValuesData @p{data} // given to the @@ -210,7 +218,7 @@ FESystem::get_data (UpdateFlags flags, // @p{data}, similar to the // storing of the // @p{base_fe_data}s. - FEValuesData *base_data=new FEValuesData(); + FEValuesData *base_data = new FEValuesData(); data->set_fe_values_data(base_no, base_data); } data->update_flags=data->update_once | data->update_each; @@ -451,7 +459,7 @@ FESystem::compute_fill (const Mapping &mapping, } } } - if (fe_data.second_flag) + if (fe_data.compute_second_derivatives) { unsigned int offset = 0; if (face_no != static_cast (-1)) @@ -1116,14 +1124,22 @@ FESystem::InternalData::InternalData(const unsigned int n_base_elements): template FESystem::InternalData::~InternalData() { + // delete pointers and set them to + // zero to avoid inadvertant use for (unsigned int i=0; i