From a099c12aefb0b4735e88e8fa55c0d1b55f18e367 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 16 Sep 2002 13:29:37 +0000 Subject: [PATCH] Change the meaning of the restriction_is_additive flag: previously it was indexed on the vector component. With non-primitive elements, this does no longer work, so index it on the shape function number. git-svn-id: https://svn.dealii.org/trunk@6428 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_base.h | 84 ++++++++--- deal.II/deal.II/include/fe/fe_system.h | 28 +++- deal.II/deal.II/source/fe/fe.cc | 49 +++---- deal.II/deal.II/source/fe/fe_dgp.cc | 27 +++- deal.II/deal.II/source/fe/fe_dgq.cc | 9 +- deal.II/deal.II/source/fe/fe_q.cc | 3 +- deal.II/deal.II/source/fe/fe_system.cc | 189 +++++++++++++++++++++---- 7 files changed, 311 insertions(+), 78 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 6d2312c7f5..e30daa0172 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -409,11 +409,17 @@ class FiniteElementBase : public Subscriptor, /** * Construct an object of this - * type. You have to set some + * type. You have to set some * member variables, for example * some matrices, explicitly * after calling this base class' - * constructor. + * constructor. For this see the + * existing finite element + * classes. For the second and + * third parameter of this + * constructor, see the document + * of the respective member + * variables. */ FiniteElementBase (const FiniteElementData &fe_data, const std::vector &restriction_is_additive_flags, @@ -816,11 +822,16 @@ class FiniteElementBase : public Subscriptor, component_to_base (unsigned int component) const; /** - * Access the @p{restriction_is_additive_flag} - * field. See there for more information on - * its contents. + * Access the + * @p{restriction_is_additive_flag} + * field. See there for more + * information on its contents. + * + * The index must be between zero + * and the number of shape + * functions of this element. */ - bool restriction_is_additive (const unsigned int component) const; + bool restriction_is_additive (const unsigned int index) const; /** * Return the support points of @@ -1187,18 +1198,20 @@ class FiniteElementBase : public Subscriptor, * * This array has valid values * also in the case of - * vector-value + * vector-valued * (i.e. non-primitive) shape * functions, in contrast to the * @p{system_to_component_table}. */ - std::vector,unsigned int> > system_to_base_table; + std::vector,unsigned int> > + system_to_base_table; /** * Likewise for the indices on * faces. */ - std::vector,unsigned int> > face_system_to_base_table; + std::vector,unsigned int> > + face_system_to_base_table; /** * Map between component and @@ -1278,9 +1291,48 @@ class FiniteElementBase : public Subscriptor, * build the complete matrix. The * flag should be @p{true}. * - * There is one flag per - * component in vector valued - * elements. + * For examples of use of these + * flags, see the places in the + * library where it is queried. + * + * There is one flag per shape + * function, indicating whether + * it belongs to the class of + * shape functions that are + * additive in the restriction or + * not. + * + * Note that in previous versions + * of the library, there was one + * flag per vector component of + * the element. This is based on + * the fact that all the shape + * functions that belong to the + * same vector component must + * necessarily behave in the same + * way, to make things + * reasonable. However, the + * problem is that it is + * sometimes impossible to query + * this flag in the vector-valued + * case: this used to be done + * with the + * @p{system_to_component_index} + * function that returns which + * vector component a shape + * function is associated + * with. The point is that since + * we now support shape functions + * that are associated with more + * than one vector component (for + * example the shape functions of + * Raviart-Thomas, or Nedelec + * elements), that function can + * no more be used, so it can be + * difficult to find out which + * for vector component we would + * like to query the + * restriction-is-additive flags. */ const std::vector restriction_is_additive_flags; @@ -1558,11 +1610,11 @@ FiniteElementBase::component_to_base (unsigned int index) const template inline bool -FiniteElementBase::restriction_is_additive (const unsigned int component) const +FiniteElementBase::restriction_is_additive (const unsigned int index) const { - Assert(componentn_components(), - ExcIndexRange(component, 0, this->n_components())); - return restriction_is_additive_flags[component]; + Assert(index < this->dofs_per_cell, + ExcIndexRange(index, 0, this->dofs_per_cell)); + return restriction_is_additive_flags[index]; } diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 04f164a91f..1247420c46 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -525,12 +525,15 @@ class FESystem : public FiniteElement /** - * Helper function used in the constructor: - * takes a @p{FiniteElement} object - * and returns an boolean vector including - * the @p{restriction_is_additive_flags} of - * the mixed element consisting of @p{N} - * elements of the sub-element @p{fe}. + * Helper function used in the + * constructor: takes a + * @p{FiniteElement} object and + * returns an boolean vector + * including the + * @p{restriction_is_additive_flags} + * of the mixed element + * consisting of @p{N} elements + * of the sub-element @p{fe}. */ static std::vector compute_restriction_is_additive_flags (const FiniteElement &fe, @@ -558,6 +561,19 @@ class FESystem : public FiniteElement const FiniteElement &fe3, const unsigned int N3); + /** + * Compute the named flags for a + * list of finite elements with + * multiplicities given in the + * second argument. This function + * is called from all the above + * functions. + */ + static std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities); + + /** * Compute the non-zero vector * components of a composed diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 6002b08e73..d0b1c0d7b9 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -108,34 +108,35 @@ FiniteElementBase::InternalDataBase::~InternalDataBase () template -FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags, - const std::vector > &nonzero_components) +FiniteElementBase:: +FiniteElementBase (const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components) : FiniteElementData (fe_data), - system_to_component_table(this->dofs_per_cell), - face_system_to_component_table(this->dofs_per_face), - system_to_base_table(this->dofs_per_cell), - face_system_to_base_table(this->dofs_per_face), - component_to_system_table(this->components, - std::vector(this->dofs_per_cell)), - face_component_to_system_table(this->components, - std::vector(this->dofs_per_face)), - component_to_base_table (this->components, - std::make_pair(0U, 0U)), - restriction_is_additive_flags(restriction_is_additive_flags), - nonzero_components (nonzero_components), - n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)), - cached_primitivity (std::find_if (n_nonzero_components_table.begin(), - n_nonzero_components_table.end(), - std::bind2nd(std::not_equal_to(), - 1U)) - == - n_nonzero_components_table.end()) + system_to_component_table (this->dofs_per_cell), + face_system_to_component_table(this->dofs_per_face), + system_to_base_table(this->dofs_per_cell), + face_system_to_base_table(this->dofs_per_face), + component_to_system_table(this->components, + std::vector(this->dofs_per_cell)), + face_component_to_system_table(this->components, + std::vector(this->dofs_per_face)), + component_to_base_table (this->components, + std::make_pair(0U, 0U)), + restriction_is_additive_flags(restriction_is_additive_flags), + nonzero_components (nonzero_components), + n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)), + cached_primitivity (std::find_if (n_nonzero_components_table.begin(), + n_nonzero_components_table.end(), + std::bind2nd(std::not_equal_to(), + 1U)) + == + n_nonzero_components_table.end()) { - Assert (restriction_is_additive_flags.size()==fe_data.components, + Assert (restriction_is_additive_flags.size() == this->dofs_per_cell, ExcDimensionMismatch(restriction_is_additive_flags.size(), - fe_data.components)); + this->dofs_per_cell)); Assert (nonzero_components.size() == this->dofs_per_cell, ExcInternalError()); for (unsigned int i=0; i FE_DGP::FE_DGP (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(1,true), + std::vector(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell,true), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - polynomial_space (Legendre::generate_complete_basis(degree)) + degree(degree), + polynomial_space (Legendre::generate_complete_basis(degree)) { // if defined, copy over matrices // from precomputed arrays @@ -41,7 +41,15 @@ FE_DGP::FE_DGP (const unsigned int degree) else for (unsigned int i=0; i::children_per_cell;++i) this->prolongation[i].reinit(0,0); - + + // restriction can be defined + // through projection for + // discontinuous elements, but is + // presently not implemented for DGP + // elements. + // + // if it were, then the following + // snippet would be the right code // if ((degree < Matrices::n_projection_matrices) && // (Matrices::projection_matrices[degree] != 0)) // { @@ -50,7 +58,16 @@ FE_DGP::FE_DGP (const unsigned int degree) // else // // matrix undefined, set size to zero // for (unsigned int i=0;i::children_per_cell;++i) -// restriction[i].reinit(0); +// restriction[i].reinit(0, 0); + // since not implemented, set to + // "empty" + for (unsigned int i=0;i::children_per_cell;++i) + restriction[i].reinit(0, 0); + + // note further, that these + // elements have neither support + // nor face-support points, so + // leave these fields empty }; diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index a638f981b0..6ddd18b930 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -28,7 +28,8 @@ template FE_DGQ::FE_DGQ (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(1,true), + std::vector(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + true), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), degree(degree), @@ -47,6 +48,12 @@ FE_DGQ::FE_DGQ (const unsigned int degree) // from precomputed arrays and // generate all other matrices by // permutations + // + // (note: the matrix is defined if + // something was entered into the + // respective table, and what was + // entered is not a NULL pointer -- + // this would allow for "holes") if ((degree < Matrices::n_embedding_matrices) && (Matrices::embedding[degree] != 0)) { diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 192fc8f0ed..e20e49d2e3 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -27,7 +27,8 @@ template FE_Q::FE_Q (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector (1,false), + std::vector (FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + false), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), degree(degree), diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 69b3ba4735..9b8e526df9 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -1680,11 +1680,13 @@ std::vector FESystem::compute_restriction_is_additive_flags (const FiniteElement &fe, const unsigned int n_elements) { - std::vector tmp; - for (unsigned int i=0; i*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe); + multiplicities.push_back (n_elements); + + return compute_restriction_is_additive_flags (fe_list, multiplicities); }; @@ -1696,14 +1698,16 @@ FESystem::compute_restriction_is_additive_flags (const FiniteElement & const FiniteElement &fe2, const unsigned int N2) { - std::vector tmp; - for (unsigned int i=0; i*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe1); + multiplicities.push_back (N1); + + fe_list.push_back (&fe2); + multiplicities.push_back (N2); + + return compute_restriction_is_additive_flags (fe_list, multiplicities); }; @@ -1717,17 +1721,151 @@ FESystem::compute_restriction_is_additive_flags (const FiniteElement & const FiniteElement &fe3, const unsigned int N3) { - std::vector tmp; - for (unsigned int i=0; i*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe1); + multiplicities.push_back (N1); + + fe_list.push_back (&fe2); + multiplicities.push_back (N2); + + fe_list.push_back (&fe3); + multiplicities.push_back (N3); + + return compute_restriction_is_additive_flags (fe_list, multiplicities); +}; + + + +template +std::vector +FESystem:: +compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities) +{ + Assert (fes.size() == multiplicities.size(), ExcInternalError()); + + // first count the number of dofs + // and components that will emerge + // from the given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; idofs_per_cell * multiplicities[i]; + + // generate the array that will + // hold the output + std::vector retval (n_shape_functions, false); + + // finally go through all the shape + // functions of the base elements, + // and copy their flags. this + // somehow copies the code in + // build_cell_table, which is not + // nice as it uses too much + // implicit knowledge about the + // layout of the individual bases + // in the composed FE, but there + // seems no way around... + // + // for each shape function, copy + // the flags from the base element + // to this one, taking into account + // multiplicities, and other + // complications + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + for(unsigned int base=0; basedofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 2. Lines + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + for (unsigned int base=0; basedofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 3. Quads + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + for (unsigned int base=0; basedofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + for(unsigned int base=0; basedofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + Assert (total_index == n_shape_functions, ExcInternalError()); + + return retval; }; @@ -1970,13 +2108,14 @@ compute_nonzero_components (const std::vector*> &fes, } } - Assert (total_index == retval.size(), ExcInternalError()); + Assert (total_index == n_shape_functions, ExcInternalError()); return retval; }; + template const FiniteElement & FESystem::base_element (const unsigned int index) const -- 2.39.5