From: guido Date: Wed, 3 Jul 2002 15:30:05 +0000 (+0000) Subject: fe sublib ISOified for icc -Xc X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=68c79089e65744c8bd59d8b49e47cb7650249dd6;p=dealii-svn.git fe sublib ISOified for icc -Xc git-svn-id: https://svn.dealii.org/trunk@6216 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h index d6bcd24230..5393e2aa27 100644 --- a/deal.II/deal.II/include/fe/fe_dgp.h +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -199,6 +199,53 @@ class FE_DGP : public FiniteElement */ virtual unsigned int memory_consumption () const; + + /** + * Declare a nested class which + * will hold static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_dgp_[123]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Pointers to the embedding + * matrices, one for each + * polynomial degree starting + * from constant elements + */ + static const double * const embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree plus one + * for which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As @p{embedding} but for + * projection matrices. + */ + static const double * const projection_matrices[][GeometryInfo::children_per_cell]; + + /** + * As + * @p{n_embedding_matrices} + * but for projection + * matrices. + */ + static const unsigned int n_projection_matrices; + }; + protected: /** @@ -264,53 +311,6 @@ class FE_DGP : public FiniteElement FEValuesData& data) const ; private: - - /** - * Declare a nested class which - * will hold static definitions of - * various matrices such as - * constraint and embedding - * matrices. The definition of - * the various static fields are - * in the files @p{fe_dgp_[123]d.cc} - * in the source directory. - */ - struct Matrices - { - /** - * Pointers to the embedding - * matrices, one for each - * polynomial degree starting - * from constant elements - */ - static const double * const embedding[][GeometryInfo::children_per_cell]; - - /** - * Number of elements (first - * index) the above field - * has. Equals the highest - * polynomial degree plus one - * for which the embedding - * matrices have been - * computed. - */ - static const unsigned int n_embedding_matrices; - - /** - * As @p{embedding} but for - * projection matrices. - */ - static const double * const projection_matrices[][GeometryInfo::children_per_cell]; - - /** - * As - * @p{n_embedding_matrices} - * but for projection - * matrices. - */ - static const unsigned int n_projection_matrices; - }; - /** * Only for internal use. Its diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index 131b0a92dd..84afc997bb 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -199,6 +199,54 @@ class FE_DGQ : public FiniteElement */ virtual unsigned int memory_consumption () const; + + /** + * Declare a nested class which + * will hold static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_dgq_[123]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Pointers to the embedding + * matrices, one for each + * polynomial degree starting + * from constant elements + */ + static const double * const embedding[]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree plus one + * for which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As @p{embedding} but for + * projection matrices. + */ + static const double * const projection_matrices[]; + + /** + * As + * @p{n_embedding_matrices} + * but for projection + * matrices. + */ + static const unsigned int n_projection_matrices; + }; + + protected: /** @@ -264,54 +312,6 @@ class FE_DGQ : public FiniteElement FEValuesData& data) const ; private: - - /** - * Declare a nested class which - * will hold static definitions of - * various matrices such as - * constraint and embedding - * matrices. The definition of - * the various static fields are - * in the files @p{fe_dgq_[123]d.cc} - * in the source directory. - */ - struct Matrices - { - /** - * Pointers to the embedding - * matrices, one for each - * polynomial degree starting - * from constant elements - */ - static const double * const embedding[]; - - /** - * Number of elements (first - * index) the above field - * has. Equals the highest - * polynomial degree plus one - * for which the embedding - * matrices have been - * computed. - */ - static const unsigned int n_embedding_matrices; - - /** - * As @p{embedding} but for - * projection matrices. - */ - static const double * const projection_matrices[]; - - /** - * As - * @p{n_embedding_matrices} - * but for projection - * matrices. - */ - static const unsigned int n_projection_matrices; - }; - - /** * Only for internal use. Its * full name is diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index d5e8f390e4..6c5717e2a4 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -191,6 +191,74 @@ class FE_Nedelec : public FiniteElement */ virtual unsigned int memory_consumption () const; + + /** + * Declare a nested class which + * will hold static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_nedelec_[23]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Embedding matrices. For + * each element type (the + * first index) there are as + * many embedding matrices as + * there are children per + * cell. The first index + * starts with linear + * elements and goes up in + * polynomial degree. The + * array may grow in the + * future with the number of + * elements for which these + * matrices have been + * computed. If for some + * element, the matrices have + * not been computed then you + * may use the element + * nevertheless but can not + * access the respective + * fields. + */ + static const double * const + embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree for + * which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As the + * @p{embedding_matrices} + * field, but for the + * interface constraints. One + * for each element for which + * it has been computed. + */ + static const double * const constraint_matrices[]; + + /** + * Like + * @p{n_embedding_matrices}, + * but for the number of + * interface constraint + * matrices. + */ + static const unsigned int n_constraint_matrices; + }; /** * Exception */ @@ -260,74 +328,6 @@ class FE_Nedelec : public FiniteElement FEValuesData& data) const ; private: - - /** - * Declare a nested class which - * will hold static definitions of - * various matrices such as - * constraint and embedding - * matrices. The definition of - * the various static fields are - * in the files @p{fe_nedelec_[23]d.cc} - * in the source directory. - */ - struct Matrices - { - /** - * Embedding matrices. For - * each element type (the - * first index) there are as - * many embedding matrices as - * there are children per - * cell. The first index - * starts with linear - * elements and goes up in - * polynomial degree. The - * array may grow in the - * future with the number of - * elements for which these - * matrices have been - * computed. If for some - * element, the matrices have - * not been computed then you - * may use the element - * nevertheless but can not - * access the respective - * fields. - */ - static const double * const - embedding[][GeometryInfo::children_per_cell]; - - /** - * Number of elements (first - * index) the above field - * has. Equals the highest - * polynomial degree for - * which the embedding - * matrices have been - * computed. - */ - static const unsigned int n_embedding_matrices; - - /** - * As the - * @p{embedding_matrices} - * field, but for the - * interface constraints. One - * for each element for which - * it has been computed. - */ - static const double * const constraint_matrices[]; - - /** - * Like - * @p{n_embedding_matrices}, - * but for the number of - * interface constraint - * matrices. - */ - static const unsigned int n_constraint_matrices; - }; /** * Only for internal use. Its diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 61abd087d4..c20d11b08b 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -403,71 +403,6 @@ class FE_Q : public FiniteElement */ virtual unsigned int memory_consumption () const; - protected: - /** - * @p{clone} function instead of - * a copy constructor. - * - * This function is needed by the - * constructors of @p{FESystem}. - */ - virtual FiniteElement * clone() const; - - /** - * Prepare internal data - * structures and fill in values - * independent of the cell. - */ - virtual - typename Mapping::InternalDataBase * - get_data (const UpdateFlags, - const Mapping& mapping, - const Quadrature& quadrature) const ; - - /** - * Implementation of the same - * function in - * @ref{FiniteElement}. - */ - virtual void - fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_internal, - typename Mapping::InternalDataBase &fe_internal, - FEValuesData& data) const; - - /** - * Implementation of the same - * function in - * @ref{FiniteElement}. - */ - virtual void - fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_internal, - typename Mapping::InternalDataBase &fe_internal, - FEValuesData& data) const ; - - /** - * Implementation of the same - * function in - * @ref{FiniteElement}. - */ - virtual void - fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_internal, - typename Mapping::InternalDataBase &fe_internal, - FEValuesData& data) const ; - - private: - /** * Declare a nested class which * will hold static definitions of @@ -535,6 +470,71 @@ class FE_Q : public FiniteElement */ static const unsigned int n_constraint_matrices; }; + + protected: + /** + * @p{clone} function instead of + * a copy constructor. + * + * This function is needed by the + * constructors of @p{FESystem}. + */ + virtual FiniteElement * clone() const; + + /** + * Prepare internal data + * structures and fill in values + * independent of the cell. + */ + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + private: /** * Only for internal use. Its diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index c1d9f86f12..8608ab103f 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -105,47 +105,6 @@ class MappingQ : public MappingQ1 */ unsigned int get_degree () const; - protected: - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values) const ; - - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values, - typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; - - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values, - typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; - /** * Storage for internal data of * Q_degree transformation. @@ -200,6 +159,47 @@ class MappingQ : public MappingQ1 typename MappingQ1::InternalData mapping_q1_data; }; + protected: + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_values (const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values) const ; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + typename std::vector > &exterior_form, + typename std::vector > &normal_vectors) const ; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + typename std::vector > &exterior_form, + typename std::vector > &normal_vectors) const ; + /** * For @p{dim=2,3}. Append the * support points of all shape diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 2ceba94765..e3ff5feacf 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -99,48 +99,6 @@ class MappingQ1 : public Mapping */ virtual UpdateFlags update_each (const UpdateFlags) const; - protected: - - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values) const ; - - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values, - typename std::vector > &boundary_form, - typename std::vector > &normal_vectors) const ; - - /** - * Implementation of the interface in - * @ref{Mapping}. - */ - virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature& quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values, - typename std::vector > &boundary_form, - typename std::vector > &normal_vectors) const ; - /** * Storage for internal data of * d-linear transformation. @@ -283,6 +241,48 @@ class MappingQ1 : public Mapping unsigned int n_shape_functions; }; + protected: + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_values (const typename DoFHandler::cell_iterator &cell, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values) const ; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + typename std::vector > &boundary_form, + typename std::vector > &normal_vectors) const ; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + typename std::vector > &boundary_form, + typename std::vector > &normal_vectors) const ; + /** * Compute shape values and/or * derivatives. diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index b74c9c74ce..78fd5c6fe2 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -113,31 +113,34 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data const std::vector > &nonzero_components) : FiniteElementData (fe_data), - system_to_component_table(dofs_per_cell), - face_system_to_component_table(dofs_per_face), - system_to_base_table(dofs_per_cell), - face_system_to_base_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, 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, ExcDimensionMismatch(restriction_is_additive_flags.size(), fe_data.components)); - Assert (nonzero_components.size() == dofs_per_cell, + Assert (nonzero_components.size() == this->dofs_per_cell, ExcInternalError()); for (unsigned int i=0; in_components(), ExcInternalError()); Assert (std::count (nonzero_components[i].begin(), nonzero_components[i].end(), @@ -146,14 +149,14 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data ExcInternalError()); Assert (n_nonzero_components_table[i] >= 1, ExcInternalError()); - Assert (n_nonzero_components_table[i] <= n_components(), + Assert (n_nonzero_components_table[i] <= this->n_components(), ExcInternalError()); }; for (unsigned int i=0; i::children_per_cell; ++i) { - restriction[i].reinit (dofs_per_cell, dofs_per_cell); - prolongation[i].reinit (dofs_per_cell, dofs_per_cell); + restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell); + prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell); }; // first set sizes of some @@ -163,38 +166,39 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data switch (dim) { case 1: - Assert ((interface_constraints.m() == 0) && - (interface_constraints.n() == 0), - ExcInternalError()); - break; + Assert ((interface_constraints.m() == 0) && + (interface_constraints.n() == 0), + ExcInternalError()); + break; case 2: - interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line, - dofs_per_face); - break; + interface_constraints.reinit (this->dofs_per_vertex + +2*this->dofs_per_line, + this->dofs_per_face); + break; case 3: - interface_constraints.reinit (5*dofs_per_vertex + - 12*dofs_per_line + - 4*dofs_per_quad, - dofs_per_face); - break; + interface_constraints.reinit (5*this->dofs_per_vertex + + 12*this->dofs_per_line + + 4*this->dofs_per_quad, + this->dofs_per_face); + break; default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); }; // this is the default way, if there is only // one component; if there are several, then // the constructor of the derived class needs // to fill these arrays - for (unsigned int j=0 ; jdofs_per_cell ; ++j) { system_to_component_table[j] = std::pair(0,j); system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); component_to_system_table[0][j] = j; } - for (unsigned int j=0 ; jdofs_per_face ; ++j) { face_system_to_component_table[j] = std::pair(0,j); face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); @@ -299,7 +303,7 @@ template const FullMatrix & FiniteElementBase::constraints () const { - Assert ((dofs_per_face == 0) || (interface_constraints.m() != 0), + Assert ((this->dofs_per_face == 0) || (interface_constraints.m() != 0), ExcConstraintsVoid()); if (dim==1) @@ -331,7 +335,7 @@ FiniteElementBase::get_unit_support_points () const // there are as many as there are // degrees of freedom Assert ((unit_support_points.size() == 0) || - (unit_support_points.size() == dofs_per_cell), + (unit_support_points.size() == this->dofs_per_cell), ExcInternalError()); return unit_support_points; }; @@ -356,7 +360,7 @@ FiniteElementBase::get_unit_face_support_points () const // there are as many as there are // degrees of freedom on a face Assert ((unit_face_support_points.size() == 0) || - (unit_face_support_points.size() == dofs_per_face), + (unit_face_support_points.size() == this->dofs_per_face), ExcInternalError()); return unit_face_support_points; }; @@ -411,7 +415,7 @@ compute_2nd (const Mapping &mapping, Assert ((fe_internal.update_each | fe_internal.update_once) & update_second_derivatives, ExcInternalError()); - Assert (data.shape_2nd_derivatives.size() == dofs_per_cell, + Assert (data.shape_2nd_derivatives.size() == this->dofs_per_cell, ExcInternalError()); // Number of quadrature points const unsigned int n_q_points = data.shape_2nd_derivatives[0].size(); @@ -435,7 +439,7 @@ compute_2nd (const Mapping &mapping, // for all shape functions at all // quadrature points and difference // quotients in all directions: - for (unsigned int shape=0; shapedofs_per_cell; ++shape) { for (unsigned int d1=0; d1::FE_DGP (const unsigned int degree) std::vector(1,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 if ((degree < Matrices::n_embedding_matrices) && (Matrices::embedding[degree][0] != 0)) for (unsigned int c=0; c::children_per_cell; ++c) - prolongation[c].fill (Matrices::embedding[degree][c]); + this->prolongation[c].fill (Matrices::embedding[degree][c]); else for (unsigned int i=0; i::children_per_cell;++i) - prolongation[i].reinit(0,0); + this->prolongation[i].reinit(0,0); // // same as above: copy over matrix // // from predefined values and @@ -72,7 +72,7 @@ double FE_DGP::shape_value (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_value(i, p); } @@ -84,7 +84,7 @@ FE_DGP::shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_value(i, p); } @@ -96,7 +96,7 @@ Tensor<1,dim> FE_DGP::shape_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_grad(i, p); } @@ -107,7 +107,7 @@ FE_DGP::shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad(i, p); } @@ -119,7 +119,7 @@ Tensor<2,dim> FE_DGP::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_grad_grad(i, p); } @@ -131,7 +131,7 @@ FE_DGP::shape_grad_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad_grad(i, p); } @@ -219,15 +219,15 @@ FE_DGP::get_data (const UpdateFlags update_flags, // allocate memory if (flags & update_values) { - values.resize (dofs_per_cell); - data->shape_values.resize(dofs_per_cell, + values.resize (this->dofs_per_cell); + data->shape_values.resize(this->dofs_per_cell, std::vector(n_q_points)); } if (flags & update_gradients) { - grads.resize (dofs_per_cell); - data->shape_gradients.resize(dofs_per_cell, + grads.resize (this->dofs_per_cell); + data->shape_gradients.resize(this->dofs_per_cell, std::vector >(n_q_points)); } @@ -250,7 +250,7 @@ FE_DGP::get_data (const UpdateFlags update_flags, { polynomial_space.compute(quadrature.point(i), values, grads, grad_grads); - for (unsigned int k=0; kdofs_per_cell; ++k) { if (flags & update_values) data->shape_values[k][i] = values[k]; @@ -284,7 +284,7 @@ FE_DGP::fill_fe_values (const Mapping &mapping, const UpdateFlags flags(fe_data.current_update_flags()); - for (unsigned int k=0; kdofs_per_cell; ++k) { if (flags & update_values) for (unsigned int i=0;i::fill_fe_face_values (const Mapping &mapping, const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::fill_fe_subface_values (const Mapping &mappi const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::FE_DGQ (const unsigned int degree) std::vector(1,true), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - polynomial_space (LagrangeEquidistant::generate_complete_basis(degree)) + degree(degree), + polynomial_space (LagrangeEquidistant::generate_complete_basis(degree)) { // generate permutation/rotation // index sets to generate some @@ -50,45 +50,45 @@ FE_DGQ::FE_DGQ (const unsigned int degree) if ((degree < Matrices::n_embedding_matrices) && (Matrices::embedding[degree] != 0)) { - prolongation[0].fill (Matrices::embedding[degree]); + this->prolongation[0].fill (Matrices::embedding[degree]); switch (dim) { case 1: - prolongation[1].fill_permutation (prolongation[0], - right, right); - break; + this->prolongation[1].fill_permutation (this->prolongation[0], + right, right); + break; case 2: - prolongation[1].fill_permutation (prolongation[0], - right, right); - prolongation[2].fill_permutation (prolongation[1], - right, right); - prolongation[3].fill_permutation (prolongation[2], - right, right); - break; + this->prolongation[1].fill_permutation (this->prolongation[0], + right, right); + this->prolongation[2].fill_permutation (this->prolongation[1], + right, right); + this->prolongation[3].fill_permutation (this->prolongation[2], + right, right); + break; case 3: - prolongation[1].fill_permutation (prolongation[0], - right, right); - prolongation[5].fill_permutation (prolongation[1], - right, right); - prolongation[4].fill_permutation (prolongation[5], - right, right); - prolongation[7].fill_permutation (prolongation[4], - top, top); - prolongation[3].fill_permutation (prolongation[7], - top, top); - prolongation[6].fill_permutation (prolongation[5], - top, top); - prolongation[2].fill_permutation (prolongation[6], - top, top); - break; + this->prolongation[1].fill_permutation (this->prolongation[0], + right, right); + this->prolongation[5].fill_permutation (this->prolongation[1], + right, right); + this->prolongation[4].fill_permutation (this->prolongation[5], + right, right); + this->prolongation[7].fill_permutation (this->prolongation[4], + top, top); + this->prolongation[3].fill_permutation (this->prolongation[7], + top, top); + this->prolongation[6].fill_permutation (this->prolongation[5], + top, top); + this->prolongation[2].fill_permutation (this->prolongation[6], + top, top); + break; default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } } else // matrix undefined, set size to zero for (unsigned int i=0;i::children_per_cell;++i) - prolongation[i].reinit(0); + this->prolongation[i].reinit(0); // same as above: copy over matrix // from predefined values and @@ -96,45 +96,45 @@ FE_DGQ::FE_DGQ (const unsigned int degree) if ((degree < Matrices::n_projection_matrices) && (Matrices::projection_matrices[degree] != 0)) { - restriction[0].fill (Matrices::projection_matrices[degree]); + this->restriction[0].fill (Matrices::projection_matrices[degree]); switch (dim) { case 1: - restriction[1].fill_permutation (restriction[0], - right, right); - break; + this->restriction[1].fill_permutation (this->restriction[0], + right, right); + break; case 2: - restriction[1].fill_permutation (restriction[0], - right, right); - restriction[2].fill_permutation (restriction[1], - right, right); - restriction[3].fill_permutation (restriction[2], - right, right); - break; + this->restriction[1].fill_permutation (this->restriction[0], + right, right); + this->restriction[2].fill_permutation (this->restriction[1], + right, right); + this->restriction[3].fill_permutation (this->restriction[2], + right, right); + break; case 3: - restriction[1].fill_permutation (restriction[0], - right, right); - restriction[5].fill_permutation (restriction[1], - right, right); - restriction[4].fill_permutation (restriction[5], - right, right); - restriction[7].fill_permutation (restriction[4], - top, top); - restriction[3].fill_permutation (restriction[7], - top, top); - restriction[6].fill_permutation (restriction[5], - top, top); - restriction[2].fill_permutation (restriction[6], - top, top); - break; + this->restriction[1].fill_permutation (this->restriction[0], + right, right); + this->restriction[5].fill_permutation (this->restriction[1], + right, right); + this->restriction[4].fill_permutation (this->restriction[5], + right, right); + this->restriction[7].fill_permutation (this->restriction[4], + top, top); + this->restriction[3].fill_permutation (this->restriction[7], + top, top); + this->restriction[6].fill_permutation (this->restriction[5], + top, top); + this->restriction[2].fill_permutation (this->restriction[6], + top, top); + break; default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } } else // matrix undefined, set size to zero for (unsigned int i=0;i::children_per_cell;++i) - restriction[i].reinit(0); + this->restriction[i].reinit(0); // finally fill in support points @@ -142,9 +142,9 @@ FE_DGQ::FE_DGQ (const unsigned int degree) { // constant elements, take // midpoint - unit_support_points.resize(1); + this->unit_support_points.resize(1); for (unsigned int i=0; iunit_support_points[0](i) = 0.5; } else { @@ -153,7 +153,7 @@ FE_DGQ::FE_DGQ (const unsigned int degree) for (unsigned int i=1; iunit_support_points.resize(n); const double step = 1./degree; Point p; @@ -169,7 +169,7 @@ FE_DGQ::FE_DGQ (const unsigned int degree) if (dim>2) p(2) = iz * step; - unit_support_points[k++] = p; + this->unit_support_points[k++] = p; }; }; @@ -193,7 +193,7 @@ double FE_DGQ::shape_value (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_value(i, p); } @@ -205,7 +205,7 @@ FE_DGQ::shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_value(i, p); } @@ -217,7 +217,7 @@ Tensor<1,dim> FE_DGQ::shape_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_grad(i, p); } @@ -228,7 +228,7 @@ FE_DGQ::shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad(i, p); } @@ -240,7 +240,7 @@ Tensor<2,dim> FE_DGQ::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); return polynomial_space.compute_grad_grad(i, p); } @@ -252,7 +252,7 @@ FE_DGQ::shape_grad_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad_grad(i, p); } @@ -331,51 +331,51 @@ FE_DGQ::rotate_indices (std::vector &numbers, // Rotate xy-plane // counter-clockwise case 'z': - for (unsigned int iz=0;iz<((dim>2) ? n:1);++iz) - for (unsigned int j=0;j2) ? n:1);++iz) + for (unsigned int j=0;j2) ? n:1);++iz) - for (unsigned int iy=0;iy2) ? n:1);++iz) + for (unsigned int iy=0;iy2, - typename FiniteElementData::ExcSpaceDimensionMismatch (dim,3)); - for (unsigned int iz=0;iz2, + typename FiniteElementData::ExcSpaceDimensionMismatch (dim,3)); + for (unsigned int iz=0;iz2, - typename FiniteElementData::ExcSpaceDimensionMismatch (dim,3)); - for (unsigned int iz=0;iz2, + typename FiniteElementData::ExcSpaceDimensionMismatch (dim,3)); + for (unsigned int iz=0;iz::get_data (const UpdateFlags update_flags, // allocate memory if (flags & update_values) { - values.resize (dofs_per_cell); - data->shape_values.resize(dofs_per_cell, + values.resize (this->dofs_per_cell); + data->shape_values.resize(this->dofs_per_cell, std::vector(n_q_points)); } if (flags & update_gradients) { - grads.resize (dofs_per_cell); - data->shape_gradients.resize(dofs_per_cell, + grads.resize (this->dofs_per_cell); + data->shape_gradients.resize(this->dofs_per_cell, std::vector >(n_q_points)); } @@ -450,7 +450,7 @@ FE_DGQ::get_data (const UpdateFlags update_flags, { polynomial_space.compute(quadrature.point(i), values, grads, grad_grads); - for (unsigned int k=0; kdofs_per_cell; ++k) { if (flags & update_values) data->shape_values[k][i] = values[k]; @@ -484,7 +484,7 @@ FE_DGQ::fill_fe_values (const Mapping &mapping, const UpdateFlags flags(fe_data.current_update_flags()); - for (unsigned int k=0; kdofs_per_cell; ++k) { if (flags & update_values) for (unsigned int i=0;i::fill_fe_face_values (const Mapping &mapping, const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::fill_fe_subface_values (const Mapping &mappi const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const { - Assert (shape_index < dofs_per_cell, - ExcIndexRange (shape_index, 0, dofs_per_cell)); + Assert (shape_index < this->dofs_per_cell, + ExcIndexRange (shape_index, 0, this->dofs_per_cell)); Assert (face_index < GeometryInfo::faces_per_cell, ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); @@ -630,57 +630,57 @@ FE_DGQ::has_support_on_face (const unsigned int shape_index, switch (dim) { case 1: - { - // in 1d, things are simple. since - // there is only one degree of - // freedom per vertex in this - // class, the first is on vertex 0 - // (==face 0 in some sense), the - // second on face 1: - return (((shape_index == 0) && (face_index == 0)) || - ((shape_index == 1) && (face_index == 1))); - }; + { + // in 1d, things are simple. since + // there is only one degree of + // freedom per vertex in this + // class, the first is on vertex 0 + // (==face 0 in some sense), the + // second on face 1: + return (((shape_index == 0) && (face_index == 0)) || + ((shape_index == 1) && (face_index == 1))); + }; case 2: - { - if (face_index==0 && shape_index < n) - return true; - if (face_index==1 && (shape_index % n) == degree) - return true; - if (face_index==2 && shape_index >= dofs_per_cell-n) - return true; - if (face_index==3 && (shape_index % n) == 0) - return true; - return false; - }; + { + if (face_index==0 && shape_index < n) + return true; + if (face_index==1 && (shape_index % n) == degree) + return true; + if (face_index==2 && shape_index >= this->dofs_per_cell-n) + return true; + if (face_index==3 && (shape_index % n) == 0) + return true; + return false; + }; case 3: - { - const unsigned int in2 = shape_index % n2; + { + const unsigned int in2 = shape_index % n2; - // y=0 - if (face_index==0 && in2 < n ) - return true; - // y=1 - if (face_index==1 && in2 >= n2-n) - return true; - // z=0 - if (face_index==2 && shape_index < n2) - return true; - // x=1 - if (face_index==3 && (shape_index % n) == n-1) - return true; - // z=1 - if (face_index==4 && shape_index >= dofs_per_cell - n2) - return true; - // x=0 - if (face_index==5 && (shape_index % n) == 0) - return true; - return false; - }; + // y=0 + if (face_index==0 && in2 < n ) + return true; + // y=1 + if (face_index==1 && in2 >= n2-n) + return true; + // z=0 + if (face_index==2 && shape_index < n2) + return true; + // x=1 + if (face_index==3 && (shape_index % n) == n-1) + return true; + // z=1 + if (face_index==4 && shape_index >= this->dofs_per_cell - n2) + return true; + // x=0 + if (face_index==5 && (shape_index % n) == 0) + return true; + return false; + }; default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } return true; } diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 0645bf2104..c8c6115012 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -33,7 +33,7 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) std::vector (dim,false), std::vector >(FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell, std::vector(dim,true))), - degree(degree) + degree(degree) { Assert (dim >= 2, ExcNotUsefulInThisDimension()); @@ -41,19 +41,19 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) // are defined. otherwise set them // to invalid size if (degreeinterface_constraints.fill (Matrices::constraint_matrices[degree-1]); else - interface_constraints.reinit(0,0); + this->interface_constraints.reinit(0,0); // next copy over embedding // matrices if they are defined if ((degree < Matrices::n_embedding_matrices+1) && (Matrices::embedding[degree-1][0] != 0)) for (unsigned int c=0; c::children_per_cell; ++c) - prolongation[c].fill (Matrices::embedding[degree-1][c]); + this->prolongation[c].fill (Matrices::embedding[degree-1][c]); else for (unsigned int i=0; i::children_per_cell;++i) - prolongation[i].reinit(0,0); + this->prolongation[i].reinit(0,0); // then fill restriction // matrices. they are hardcoded for @@ -61,125 +61,125 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) switch (dim) { case 2: // 2d - { - switch (degree) - { - case 1: - { - // DoF on bottom line - // of coarse cell will - // be mean value of - // bottom DoFs on the - // two adjacent child - // cells - restriction[0](0,0) = 0.5; - restriction[1](0,0) = 0.5; - // same for other DoFs - restriction[1](1,1) = 0.5; - restriction[2](1,1) = 0.5; - - restriction[2](2,2) = 0.5; - restriction[3](2,2) = 0.5; - - restriction[3](3,3) = 0.5; - restriction[0](3,3) = 0.5; - - break; - }; + { + switch (degree) + { + case 1: + { + // DoF on bottom line + // of coarse cell will + // be mean value of + // bottom DoFs on the + // two adjacent child + // cells + this->restriction[0](0,0) = 0.5; + this->restriction[1](0,0) = 0.5; + // same for other DoFs + this->restriction[1](1,1) = 0.5; + this->restriction[2](1,1) = 0.5; + + this->restriction[2](2,2) = 0.5; + this->restriction[3](2,2) = 0.5; + + this->restriction[3](3,3) = 0.5; + this->restriction[0](3,3) = 0.5; + + break; + }; - default: - { - // in case we don't - // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - restriction[i].reinit(0,0); - }; - }; + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + }; - break; - }; + break; + }; case 3: // 3d - { - switch (degree) - { - case 1: - { - // same principle as in - // 2d - restriction[0](0,0) = 0.5; - restriction[1](0,0) = 0.5; + { + switch (degree) + { + case 1: + { + // same principle as in + // 2d + this->restriction[0](0,0) = 0.5; + this->restriction[1](0,0) = 0.5; - restriction[1](1,1) = 0.5; - restriction[2](1,1) = 0.5; + this->restriction[1](1,1) = 0.5; + this->restriction[2](1,1) = 0.5; - restriction[2](2,2) = 0.5; - restriction[3](2,2) = 0.5; + this->restriction[2](2,2) = 0.5; + this->restriction[3](2,2) = 0.5; - restriction[3](3,3) = 0.5; - restriction[0](3,3) = 0.5; + this->restriction[3](3,3) = 0.5; + this->restriction[0](3,3) = 0.5; - restriction[4](4,4) = 0.5; - restriction[5](4,4) = 0.5; + this->restriction[4](4,4) = 0.5; + this->restriction[5](4,4) = 0.5; - restriction[5](5,5) = 0.5; - restriction[6](5,5) = 0.5; + this->restriction[5](5,5) = 0.5; + this->restriction[6](5,5) = 0.5; - restriction[6](6,6) = 0.5; - restriction[7](6,6) = 0.5; + this->restriction[6](6,6) = 0.5; + this->restriction[7](6,6) = 0.5; - restriction[7](7,7) = 0.5; - restriction[4](7,7) = 0.5; + this->restriction[7](7,7) = 0.5; + this->restriction[4](7,7) = 0.5; - restriction[1](8,8) = 0.5; - restriction[5](8,8) = 0.5; + this->restriction[1](8,8) = 0.5; + this->restriction[5](8,8) = 0.5; - restriction[2](9,9) = 0.5; - restriction[6](9,9) = 0.5; + this->restriction[2](9,9) = 0.5; + this->restriction[6](9,9) = 0.5; - restriction[3](10,10) = 0.5; - restriction[7](10,10) = 0.5; + this->restriction[3](10,10) = 0.5; + this->restriction[7](10,10) = 0.5; - restriction[0](11,11) = 0.5; - restriction[5](11,11) = 0.5; + this->restriction[0](11,11) = 0.5; + this->restriction[5](11,11) = 0.5; - break; - }; + break; + }; - default: - { - // in case we don't - // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - restriction[i].reinit(0,0); - }; - }; + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + }; - break; - }; + break; + }; default: - Assert (false,ExcNotImplemented()); + Assert (false,ExcNotImplemented()); } // finally fill in support points @@ -205,126 +205,126 @@ FE_Nedelec::shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); switch (dim) { case 2: // 2D - { - switch (degree) - { - // first order Nedelec - // elements - case 1: + { + switch (degree) + { + // first order Nedelec + // elements + case 1: + { + switch (i) { - switch (i) - { - // (1-y, 0) - case 0: return (component == 0 ? 1-p(1) : 0); - // (0,x) - case 1: return (component == 0 ? 0 : p(0)); - // (y, 0) - case 2: return (component == 0 ? p(1) : 0); - // (0, 1-x) - case 3: return (component == 0 ? 0 : 1-p(0)); - - // there are - // only four - // shape - // functions!? - default: - Assert (false, ExcInternalError()); - return 0; - }; + // (1-y, 0) + case 0: return (component == 0 ? 1-p(1) : 0); + // (0,x) + case 1: return (component == 0 ? 0 : p(0)); + // (y, 0) + case 2: return (component == 0 ? p(1) : 0); + // (0, 1-x) + case 3: return (component == 0 ? 0 : 1-p(0)); + + // there are + // only four + // shape + // functions!? + default: + Assert (false, ExcInternalError()); + return 0; }; + }; - // no other degrees - // implemented - default: - Assert (false, ExcNotImplemented()); - }; - }; + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; case 3: // 3D - { - switch (degree) - { - // first order Nedelec - // elements - case 1: - { - // note that the - // degrees of freedom - // on opposite faces - // have a common vector - // direction, so simply - // that a little. these - // directions are: - // - // for lines 0, 2, 4, 6: - // (1,0,0) - // for lines 1, 3, 5, 7: - // (0,0,1) - // for lines 8, 9, 10, 11: - // (0,1,0) - // - // thus, sort out all - // those cases where - // the component is - // zero anyway, and - // only otherwise - // compute the - // spatially dependent - // part which is then - // also the return - // value - if (((i<8) && (((i%2==0) && (component!=0)) || - ((i%2==1) && (component!=2)))) || - ((i>=8) && (component != 1))) - return 0; - - // now we know that the - // only non-zero - // component is - // requested: + { + switch (degree) + { + // first order Nedelec + // elements + case 1: + { + // note that the + // degrees of freedom + // on opposite faces + // have a common vector + // direction, so simply + // that a little. these + // directions are: + // + // for lines 0, 2, 4, 6: + // (1,0,0) + // for lines 1, 3, 5, 7: + // (0,0,1) + // for lines 8, 9, 10, 11: + // (0,1,0) + // + // thus, sort out all + // those cases where + // the component is + // zero anyway, and + // only otherwise + // compute the + // spatially dependent + // part which is then + // also the return + // value + if (((i<8) && (((i%2==0) && (component!=0)) || + ((i%2==1) && (component!=2)))) || + ((i>=8) && (component != 1))) + return 0; + + // now we know that the + // only non-zero + // component is + // requested: //[Anna]:checked: OK - const double x = p(0), - y = p(1), - z = p(2); - switch (i) - { - case 0: return (1-y)*(1-z); - case 2: return (1-y)*z; - case 1: return x*(1-y); - case 3: return (1-x)*(1-y); - - case 4: return y*(1-z); - case 6: return y*z; - case 5: return x*y; - case 7: return (1-x)*y; + const double x = p(0), + y = p(1), + z = p(2); + switch (i) + { + case 0: return (1-y)*(1-z); + case 2: return (1-y)*z; + case 1: return x*(1-y); + case 3: return (1-x)*(1-y); + + case 4: return y*(1-z); + case 6: return y*z; + case 5: return x*y; + case 7: return (1-x)*y; - case 8: return (1-x)*(1-z); - case 9: return x*(1-z); - case 10: return x*z; - case 11: return (1-x)*z; - default: - Assert (false, ExcInternalError()); - return 0; - }; + case 8: return (1-x)*(1-z); + case 9: return x*(1-z); + case 10: return x*z; + case 11: return (1-x)*z; + default: + Assert (false, ExcInternalError()); + return 0; }; + }; - // no other degrees - // implemented - default: - Assert (false, ExcNotImplemented()); - }; - }; + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; - // presently no other space - // dimension implemented + // presently no other space + // dimension implemented default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); }; return 0; @@ -338,54 +338,54 @@ FE_Nedelec::shape_grad_component (const unsigned int i, const Point &/*p*/, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); switch (dim) { case 2: // 2D - { - switch (degree) - { - // first order Nedelec - // elements - case 1: - { - // on the unit cell, - // the gradients of - // these shape - // functions are - // constant, so we pack - // them into a table - // for simpler lookup - // - // the format is: first - // index=shape function - // number; second - // index=vector - // component, thrid - // index=component - // within gradient + { + switch (degree) + { + // first order Nedelec + // elements + case 1: + { + // on the unit cell, + // the gradients of + // these shape + // functions are + // constant, so we pack + // them into a table + // for simpler lookup + // + // the format is: first + // index=shape function + // number; second + // index=vector + // component, thrid + // index=component + // within gradient //[Anna]: checked:OK - static const double unit_gradients[4][2][2] - = { { {0.,-1.}, {0.,0.} }, - { {0.,0.}, {1.,0.} }, - { {0.,+1.}, {0.,0.} }, - { {0.,0.}, {-1.,0.} } }; - return Tensor<1,dim>(unit_gradients[i][component]); - }; - - // no other degrees - // implemented - default: - Assert (false, ExcNotImplemented()); - }; - }; + static const double unit_gradients[4][2][2] + = { { {0.,-1.}, {0.,0.} }, + { {0.,0.}, {1.,0.} }, + { {0.,+1.}, {0.,0.} }, + { {0.,0.}, {-1.,0.} } }; + return Tensor<1,dim>(unit_gradients[i][component]); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; - // presently no other space - // dimension implemented + // presently no other space + // dimension implemented default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); }; return Tensor<1,dim>(); @@ -399,56 +399,56 @@ FE_Nedelec::shape_grad_grad_component (const unsigned int i, const Point &/*p*/, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); switch (dim) { case 2: // 2D - { - switch (degree) - { - // first order Nedelec - // elements. their second - // derivatives on the - // unit cell are zero - case 1: - { - return Tensor<2,dim>(); - }; - - // no other degrees - // implemented - default: - Assert (false, ExcNotImplemented()); - }; - }; + { + switch (degree) + { + // first order Nedelec + // elements. their second + // derivatives on the + // unit cell are zero + case 1: + { + return Tensor<2,dim>(); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; case 3: // 3D - { - switch (degree) - { - // first order Nedelec - // elements. their second - // derivatives on the - // unit cell are zero - case 1: - { - return Tensor<2,dim>(); - }; - - // no other degrees - // implemented - default: - Assert (false, ExcNotImplemented()); - }; - }; + { + switch (degree) + { + // first order Nedelec + // elements. their second + // derivatives on the + // unit cell are zero + case 1: + { + return Tensor<2,dim>(); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; - // presently no other space - // dimension implemented + // presently no other space + // dimension implemented default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); }; return Tensor<2,dim>(); @@ -470,7 +470,7 @@ void FE_Nedelec::initialize_unit_support_points () // all degrees of freedom are on // edges, and their order is the // same as the edges themselves - unit_support_points.resize(GeometryInfo::lines_per_cell * degree); + this->unit_support_points.resize(GeometryInfo::lines_per_cell * degree); unsigned int index = 0; for (unsigned int line=0; line::lines_per_cell; ++line) { @@ -486,7 +486,7 @@ void FE_Nedelec::initialize_unit_support_points () // between the vertices of each // line for (unsigned int d=0; dunit_support_points[index] = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1); }; }; @@ -510,7 +510,7 @@ void FE_Nedelec::initialize_unit_face_support_points () // is this correct? all DoFs on lines, none on faces or bubbles? // do this the same as above, but // for one dimension less - unit_face_support_points.resize(GeometryInfo::lines_per_cell * degree); + this->unit_face_support_points.resize(GeometryInfo::lines_per_cell * degree); unsigned int index = 0; for (unsigned int line=0; line::lines_per_cell; ++line) { @@ -526,7 +526,7 @@ void FE_Nedelec::initialize_unit_face_support_points () // between the vertices of each // line for (unsigned int d=0; dunit_face_support_points[index] = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1); }; }; @@ -823,8 +823,8 @@ bool FE_Nedelec::has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const { - Assert (shape_index < dofs_per_cell, - ExcIndexRange (shape_index, 0, dofs_per_cell)); + Assert (shape_index < this->dofs_per_cell, + ExcIndexRange (shape_index, 0, this->dofs_per_cell)); Assert (face_index < GeometryInfo::faces_per_cell, ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 097849cb68..6e8f6bb684 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -30,11 +30,11 @@ FE_Q::FE_Q (const unsigned int degree) std::vector (1,false), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - renumber(dofs_per_cell, 0), - renumber_inverse(dofs_per_cell, 0), - face_renumber(dofs_per_face, 0), - polynomial_space(LagrangeEquidistant::generate_complete_basis(degree)) + degree(degree), + renumber(this->dofs_per_cell, 0), + renumber_inverse(this->dofs_per_cell, 0), + face_renumber(this->dofs_per_face, 0), + polynomial_space(LagrangeEquidistant::generate_complete_basis(degree)) { // do some internal book-keeping on // cells and faces. if in 1d, the @@ -44,26 +44,26 @@ FE_Q::FE_Q (const unsigned int degree) // build inverse of renumbering // vector - for (unsigned int i=0; idofs_per_cell; ++i) renumber_inverse[renumber[i]]=i; // copy constraint matrices if they // are defined. otherwise set them // to invalid size if ((dim>1) && (degreeinterface_constraints.fill (Matrices::constraint_matrices[degree-1]); else - interface_constraints.reinit(0,0); + this->interface_constraints.reinit(0,0); // next copy over embedding // matrices if they are defined if ((degree < Matrices::n_embedding_matrices+1) && (Matrices::embedding[degree-1][0] != 0)) for (unsigned int c=0; c::children_per_cell; ++c) - prolongation[c].fill (Matrices::embedding[degree-1][c]); + this->prolongation[c].fill (Matrices::embedding[degree-1][c]); else for (unsigned int i=0; i::children_per_cell;++i) - prolongation[i].reinit(0,0); + this->prolongation[i].reinit(0,0); // then fill restriction // matrices. they are hardcoded for @@ -79,545 +79,545 @@ FE_Q::FE_Q (const unsigned int degree) switch (dim) { case 1: // 1d - { - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](2,1) = 1; - restriction[1](1,1) = 1; - restriction[1](1,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[1](1,1) = 1; - restriction[1](3,2) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[0](3,1) = 1; - restriction[1](1,1) = 1; - restriction[1](3,0) = 1; - restriction[1](4,3) = 1; - break; + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](2,1) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](1,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](2,3) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](3,2) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](2,3) = 1; + this->restriction[0](3,1) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](3,0) = 1; + this->restriction[1](4,3) = 1; + break; - default: - { - // in case we don't - // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - restriction[i].reinit(0,0); - }; - } - break; - }; + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; case 2: // 2d - { - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](4,1) = 1; - restriction[0](7,3) = 1; - restriction[0](8,2) = 1; - restriction[1](1,1) = 1; - restriction[1](4,0) = 1; - restriction[1](5,2) = 1; - restriction[1](8,3) = 1; - restriction[2](2,2) = 1; - restriction[2](5,1) = 1; - restriction[2](6,3) = 1; - restriction[2](8,0) = 1; - restriction[3](3,3) = 1; - restriction[3](6,2) = 1; - restriction[3](7,0) = 1; - restriction[3](8,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](10,11) = 1; - restriction[0](12,15) = 1; - restriction[1](1,1) = 1; - restriction[1](5,4) = 1; - restriction[1](6,7) = 1; - restriction[1](13,14) = 1; - restriction[2](2,2) = 1; - restriction[2](7,6) = 1; - restriction[2](9,8) = 1; - restriction[2](15,12) = 1; - restriction[3](3,3) = 1; - restriction[3](8,9) = 1; - restriction[3](11,10) = 1; - restriction[3](14,13) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](5,1) = 1; - restriction[0](13,14) = 1; - restriction[0](14,3) = 1; - restriction[0](16,20) = 1; - restriction[0](17,8) = 1; - restriction[0](19,11) = 1; - restriction[0](20,2) = 1; - restriction[1](1,1) = 1; - restriction[1](5,0) = 1; - restriction[1](6,5) = 1; - restriction[1](7,8) = 1; - restriction[1](8,2) = 1; - restriction[1](17,14) = 1; - restriction[1](18,20) = 1; - restriction[1](20,3) = 1; - restriction[1](21,11) = 1; - restriction[2](2,2) = 1; - restriction[2](8,1) = 1; - restriction[2](9,8) = 1; - restriction[2](11,3) = 1; - restriction[2](12,11) = 1; - restriction[2](20,0) = 1; - restriction[2](21,5) = 1; - restriction[2](23,14) = 1; - restriction[2](24,20) = 1; - restriction[3](3,3) = 1; - restriction[3](10,11) = 1; - restriction[3](11,2) = 1; - restriction[3](14,0) = 1; - restriction[3](15,14) = 1; - restriction[3](19,5) = 1; - restriction[3](20,1) = 1; - restriction[3](22,20) = 1; - restriction[3](23,8) = 1; - break; - - default: - { - // in case we don't - // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - restriction[i].reinit(0,0); - }; - } - break; - }; + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + this->restriction[2](2,2) = 1; + this->restriction[3](3,3) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](4,1) = 1; + this->restriction[0](7,3) = 1; + this->restriction[0](8,2) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](4,0) = 1; + this->restriction[1](5,2) = 1; + this->restriction[1](8,3) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](5,1) = 1; + this->restriction[2](6,3) = 1; + this->restriction[2](8,0) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](6,2) = 1; + this->restriction[3](7,0) = 1; + this->restriction[3](8,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](4,5) = 1; + this->restriction[0](10,11) = 1; + this->restriction[0](12,15) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](5,4) = 1; + this->restriction[1](6,7) = 1; + this->restriction[1](13,14) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](7,6) = 1; + this->restriction[2](9,8) = 1; + this->restriction[2](15,12) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](8,9) = 1; + this->restriction[3](11,10) = 1; + this->restriction[3](14,13) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](4,5) = 1; + this->restriction[0](5,1) = 1; + this->restriction[0](13,14) = 1; + this->restriction[0](14,3) = 1; + this->restriction[0](16,20) = 1; + this->restriction[0](17,8) = 1; + this->restriction[0](19,11) = 1; + this->restriction[0](20,2) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](5,0) = 1; + this->restriction[1](6,5) = 1; + this->restriction[1](7,8) = 1; + this->restriction[1](8,2) = 1; + this->restriction[1](17,14) = 1; + this->restriction[1](18,20) = 1; + this->restriction[1](20,3) = 1; + this->restriction[1](21,11) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](8,1) = 1; + this->restriction[2](9,8) = 1; + this->restriction[2](11,3) = 1; + this->restriction[2](12,11) = 1; + this->restriction[2](20,0) = 1; + this->restriction[2](21,5) = 1; + this->restriction[2](23,14) = 1; + this->restriction[2](24,20) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](10,11) = 1; + this->restriction[3](11,2) = 1; + this->restriction[3](14,0) = 1; + this->restriction[3](15,14) = 1; + this->restriction[3](19,5) = 1; + this->restriction[3](20,1) = 1; + this->restriction[3](22,20) = 1; + this->restriction[3](23,8) = 1; + break; + + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; case 3: // 3d - { - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - restriction[4](4,4) = 1; - restriction[5](5,5) = 1; - restriction[6](6,6) = 1; - restriction[7](7,7) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](8,1) = 1; - restriction[0](11,3) = 1; - restriction[0](16,4) = 1; - restriction[0](20,2) = 1; - restriction[0](22,5) = 1; - restriction[0](25,7) = 1; - restriction[0](26,6) = 1; - restriction[1](1,1) = 1; - restriction[1](8,0) = 1; - restriction[1](9,2) = 1; - restriction[1](17,5) = 1; - restriction[1](20,3) = 1; - restriction[1](22,4) = 1; - restriction[1](23,6) = 1; - restriction[1](26,7) = 1; - restriction[2](2,2) = 1; - restriction[2](9,1) = 1; - restriction[2](10,3) = 1; - restriction[2](18,6) = 1; - restriction[2](20,0) = 1; - restriction[2](23,5) = 1; - restriction[2](24,7) = 1; - restriction[2](26,4) = 1; - restriction[3](3,3) = 1; - restriction[3](10,2) = 1; - restriction[3](11,0) = 1; - restriction[3](19,7) = 1; - restriction[3](20,1) = 1; - restriction[3](24,6) = 1; - restriction[3](25,4) = 1; - restriction[3](26,5) = 1; - restriction[4](4,4) = 1; - restriction[4](12,5) = 1; - restriction[4](15,7) = 1; - restriction[4](16,0) = 1; - restriction[4](21,6) = 1; - restriction[4](22,1) = 1; - restriction[4](25,3) = 1; - restriction[4](26,2) = 1; - restriction[5](5,5) = 1; - restriction[5](12,4) = 1; - restriction[5](13,6) = 1; - restriction[5](17,1) = 1; - restriction[5](21,7) = 1; - restriction[5](22,0) = 1; - restriction[5](23,2) = 1; - restriction[5](26,3) = 1; - restriction[6](6,6) = 1; - restriction[6](13,5) = 1; - restriction[6](14,7) = 1; - restriction[6](18,2) = 1; - restriction[6](21,4) = 1; - restriction[6](23,1) = 1; - restriction[6](24,3) = 1; - restriction[6](26,0) = 1; - restriction[7](7,7) = 1; - restriction[7](14,6) = 1; - restriction[7](15,4) = 1; - restriction[7](19,3) = 1; - restriction[7](21,5) = 1; - restriction[7](24,2) = 1; - restriction[7](25,0) = 1; - restriction[7](26,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](14,15) = 1; - restriction[0](24,25) = 1; - restriction[0](32,35) = 1; - restriction[0](40,43) = 1; - restriction[0](52,55) = 1; - restriction[0](56,63) = 1; - restriction[1](1,1) = 1; - restriction[1](9,8) = 1; - restriction[1](10,11) = 1; - restriction[1](26,27) = 1; - restriction[1](33,34) = 1; - restriction[1](41,42) = 1; - restriction[1](44,47) = 1; - restriction[1](57,62) = 1; - restriction[2](2,2) = 1; - restriction[2](11,10) = 1; - restriction[2](13,12) = 1; - restriction[2](28,29) = 1; - restriction[2](35,32) = 1; - restriction[2](46,45) = 1; - restriction[2](49,50) = 1; - restriction[2](61,58) = 1; - restriction[3](3,3) = 1; - restriction[3](12,13) = 1; - restriction[3](15,14) = 1; - restriction[3](30,31) = 1; - restriction[3](34,33) = 1; - restriction[3](48,51) = 1; - restriction[3](54,53) = 1; - restriction[3](60,59) = 1; - restriction[4](4,4) = 1; - restriction[4](16,17) = 1; - restriction[4](22,23) = 1; - restriction[4](25,24) = 1; - restriction[4](36,39) = 1; - restriction[4](42,41) = 1; - restriction[4](53,54) = 1; - restriction[4](58,61) = 1; - restriction[5](5,5) = 1; - restriction[5](17,16) = 1; - restriction[5](18,19) = 1; - restriction[5](27,26) = 1; - restriction[5](37,38) = 1; - restriction[5](43,40) = 1; - restriction[5](45,46) = 1; - restriction[5](59,60) = 1; - restriction[6](6,6) = 1; - restriction[6](19,18) = 1; - restriction[6](21,20) = 1; - restriction[6](29,28) = 1; - restriction[6](39,36) = 1; - restriction[6](47,44) = 1; - restriction[6](51,48) = 1; - restriction[6](63,56) = 1; - restriction[7](7,7) = 1; - restriction[7](20,21) = 1; - restriction[7](23,22) = 1; - restriction[7](31,30) = 1; - restriction[7](38,37) = 1; - restriction[7](50,49) = 1; - restriction[7](55,52) = 1; - restriction[7](62,57) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](9,1) = 1; - restriction[0](17,18) = 1; - restriction[0](18,3) = 1; - restriction[0](32,33) = 1; - restriction[0](33,4) = 1; - restriction[0](44,48) = 1; - restriction[0](45,12) = 1; - restriction[0](47,15) = 1; - restriction[0](48,2) = 1; - restriction[0](62,66) = 1; - restriction[0](63,36) = 1; - restriction[0](65,21) = 1; - restriction[0](66,5) = 1; - restriction[0](89,93) = 1; - restriction[0](90,30) = 1; - restriction[0](92,42) = 1; - restriction[0](93,7) = 1; - restriction[0](98,111) = 1; - restriction[0](99,75) = 1; - restriction[0](101,57) = 1; - restriction[0](102,24) = 1; - restriction[0](107,84) = 1; - restriction[0](108,39) = 1; - restriction[0](110,27) = 1; - restriction[0](111,6) = 1; - restriction[1](1,1) = 1; - restriction[1](9,0) = 1; - restriction[1](10,9) = 1; - restriction[1](11,12) = 1; - restriction[1](12,2) = 1; - restriction[1](35,36) = 1; - restriction[1](36,5) = 1; - restriction[1](45,18) = 1; - restriction[1](46,48) = 1; - restriction[1](48,3) = 1; - restriction[1](49,15) = 1; - restriction[1](63,33) = 1; - restriction[1](64,66) = 1; - restriction[1](66,4) = 1; - restriction[1](67,21) = 1; - restriction[1](71,75) = 1; - restriction[1](72,24) = 1; - restriction[1](74,39) = 1; - restriction[1](75,6) = 1; - restriction[1](99,93) = 1; - restriction[1](100,111) = 1; - restriction[1](102,30) = 1; - restriction[1](103,57) = 1; - restriction[1](108,42) = 1; - restriction[1](109,84) = 1; - restriction[1](111,7) = 1; - restriction[1](112,27) = 1; - restriction[2](2,2) = 1; - restriction[2](12,1) = 1; - restriction[2](13,12) = 1; - restriction[2](15,3) = 1; - restriction[2](16,15) = 1; - restriction[2](38,39) = 1; - restriction[2](39,6) = 1; - restriction[2](48,0) = 1; - restriction[2](49,9) = 1; - restriction[2](51,18) = 1; - restriction[2](52,48) = 1; - restriction[2](74,36) = 1; - restriction[2](75,5) = 1; - restriction[2](77,75) = 1; - restriction[2](78,24) = 1; - restriction[2](81,42) = 1; - restriction[2](82,84) = 1; - restriction[2](84,7) = 1; - restriction[2](85,27) = 1; - restriction[2](108,33) = 1; - restriction[2](109,66) = 1; - restriction[2](111,4) = 1; - restriction[2](112,21) = 1; - restriction[2](117,93) = 1; - restriction[2](118,111) = 1; - restriction[2](120,30) = 1; - restriction[2](121,57) = 1; - restriction[3](3,3) = 1; - restriction[3](14,15) = 1; - restriction[3](15,2) = 1; - restriction[3](18,0) = 1; - restriction[3](19,18) = 1; - restriction[3](41,42) = 1; - restriction[3](42,7) = 1; - restriction[3](47,9) = 1; - restriction[3](48,1) = 1; - restriction[3](50,48) = 1; - restriction[3](51,12) = 1; - restriction[3](80,84) = 1; - restriction[3](81,39) = 1; - restriction[3](83,27) = 1; - restriction[3](84,6) = 1; - restriction[3](92,33) = 1; - restriction[3](93,4) = 1; - restriction[3](95,93) = 1; - restriction[3](96,30) = 1; - restriction[3](107,66) = 1; - restriction[3](108,36) = 1; - restriction[3](110,21) = 1; - restriction[3](111,5) = 1; - restriction[3](116,111) = 1; - restriction[3](117,75) = 1; - restriction[3](119,57) = 1; - restriction[3](120,24) = 1; - restriction[4](4,4) = 1; - restriction[4](20,21) = 1; - restriction[4](21,5) = 1; - restriction[4](29,30) = 1; - restriction[4](30,7) = 1; - restriction[4](33,0) = 1; - restriction[4](34,33) = 1; - restriction[4](53,57) = 1; - restriction[4](54,24) = 1; - restriction[4](56,27) = 1; - restriction[4](57,6) = 1; - restriction[4](65,9) = 1; - restriction[4](66,1) = 1; - restriction[4](68,66) = 1; - restriction[4](69,36) = 1; - restriction[4](90,18) = 1; - restriction[4](91,93) = 1; - restriction[4](93,3) = 1; - restriction[4](94,42) = 1; - restriction[4](101,48) = 1; - restriction[4](102,12) = 1; - restriction[4](104,111) = 1; - restriction[4](105,75) = 1; - restriction[4](110,15) = 1; - restriction[4](111,2) = 1; - restriction[4](113,84) = 1; - restriction[4](114,39) = 1; - restriction[5](5,5) = 1; - restriction[5](21,4) = 1; - restriction[5](22,21) = 1; - restriction[5](23,24) = 1; - restriction[5](24,6) = 1; - restriction[5](36,1) = 1; - restriction[5](37,36) = 1; - restriction[5](54,30) = 1; - restriction[5](55,57) = 1; - restriction[5](57,7) = 1; - restriction[5](58,27) = 1; - restriction[5](66,0) = 1; - restriction[5](67,9) = 1; - restriction[5](69,33) = 1; - restriction[5](70,66) = 1; - restriction[5](72,12) = 1; - restriction[5](73,75) = 1; - restriction[5](75,2) = 1; - restriction[5](76,39) = 1; - restriction[5](102,18) = 1; - restriction[5](103,48) = 1; - restriction[5](105,93) = 1; - restriction[5](106,111) = 1; - restriction[5](111,3) = 1; - restriction[5](112,15) = 1; - restriction[5](114,42) = 1; - restriction[5](115,84) = 1; - restriction[6](6,6) = 1; - restriction[6](24,5) = 1; - restriction[6](25,24) = 1; - restriction[6](27,7) = 1; - restriction[6](28,27) = 1; - restriction[6](39,2) = 1; - restriction[6](40,39) = 1; - restriction[6](57,4) = 1; - restriction[6](58,21) = 1; - restriction[6](60,30) = 1; - restriction[6](61,57) = 1; - restriction[6](75,1) = 1; - restriction[6](76,36) = 1; - restriction[6](78,12) = 1; - restriction[6](79,75) = 1; - restriction[6](84,3) = 1; - restriction[6](85,15) = 1; - restriction[6](87,42) = 1; - restriction[6](88,84) = 1; - restriction[6](111,0) = 1; - restriction[6](112,9) = 1; - restriction[6](114,33) = 1; - restriction[6](115,66) = 1; - restriction[6](120,18) = 1; - restriction[6](121,48) = 1; - restriction[6](123,93) = 1; - restriction[6](124,111) = 1; - restriction[7](7,7) = 1; - restriction[7](26,27) = 1; - restriction[7](27,6) = 1; - restriction[7](30,4) = 1; - restriction[7](31,30) = 1; - restriction[7](42,3) = 1; - restriction[7](43,42) = 1; - restriction[7](56,21) = 1; - restriction[7](57,5) = 1; - restriction[7](59,57) = 1; - restriction[7](60,24) = 1; - restriction[7](83,15) = 1; - restriction[7](84,2) = 1; - restriction[7](86,84) = 1; - restriction[7](87,39) = 1; - restriction[7](93,0) = 1; - restriction[7](94,33) = 1; - restriction[7](96,18) = 1; - restriction[7](97,93) = 1; - restriction[7](110,9) = 1; - restriction[7](111,1) = 1; - restriction[7](113,66) = 1; - restriction[7](114,36) = 1; - restriction[7](119,48) = 1; - restriction[7](120,12) = 1; - restriction[7](122,111) = 1; - restriction[7](123,75) = 1; - break; - default: - { - // in case we don't - // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - restriction[i].reinit(0,0); - }; - } - break; - }; + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + this->restriction[2](2,2) = 1; + this->restriction[3](3,3) = 1; + this->restriction[4](4,4) = 1; + this->restriction[5](5,5) = 1; + this->restriction[6](6,6) = 1; + this->restriction[7](7,7) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](8,1) = 1; + this->restriction[0](11,3) = 1; + this->restriction[0](16,4) = 1; + this->restriction[0](20,2) = 1; + this->restriction[0](22,5) = 1; + this->restriction[0](25,7) = 1; + this->restriction[0](26,6) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](8,0) = 1; + this->restriction[1](9,2) = 1; + this->restriction[1](17,5) = 1; + this->restriction[1](20,3) = 1; + this->restriction[1](22,4) = 1; + this->restriction[1](23,6) = 1; + this->restriction[1](26,7) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](9,1) = 1; + this->restriction[2](10,3) = 1; + this->restriction[2](18,6) = 1; + this->restriction[2](20,0) = 1; + this->restriction[2](23,5) = 1; + this->restriction[2](24,7) = 1; + this->restriction[2](26,4) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](10,2) = 1; + this->restriction[3](11,0) = 1; + this->restriction[3](19,7) = 1; + this->restriction[3](20,1) = 1; + this->restriction[3](24,6) = 1; + this->restriction[3](25,4) = 1; + this->restriction[3](26,5) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](12,5) = 1; + this->restriction[4](15,7) = 1; + this->restriction[4](16,0) = 1; + this->restriction[4](21,6) = 1; + this->restriction[4](22,1) = 1; + this->restriction[4](25,3) = 1; + this->restriction[4](26,2) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](12,4) = 1; + this->restriction[5](13,6) = 1; + this->restriction[5](17,1) = 1; + this->restriction[5](21,7) = 1; + this->restriction[5](22,0) = 1; + this->restriction[5](23,2) = 1; + this->restriction[5](26,3) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](13,5) = 1; + this->restriction[6](14,7) = 1; + this->restriction[6](18,2) = 1; + this->restriction[6](21,4) = 1; + this->restriction[6](23,1) = 1; + this->restriction[6](24,3) = 1; + this->restriction[6](26,0) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](14,6) = 1; + this->restriction[7](15,4) = 1; + this->restriction[7](19,3) = 1; + this->restriction[7](21,5) = 1; + this->restriction[7](24,2) = 1; + this->restriction[7](25,0) = 1; + this->restriction[7](26,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](8,9) = 1; + this->restriction[0](14,15) = 1; + this->restriction[0](24,25) = 1; + this->restriction[0](32,35) = 1; + this->restriction[0](40,43) = 1; + this->restriction[0](52,55) = 1; + this->restriction[0](56,63) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](9,8) = 1; + this->restriction[1](10,11) = 1; + this->restriction[1](26,27) = 1; + this->restriction[1](33,34) = 1; + this->restriction[1](41,42) = 1; + this->restriction[1](44,47) = 1; + this->restriction[1](57,62) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](11,10) = 1; + this->restriction[2](13,12) = 1; + this->restriction[2](28,29) = 1; + this->restriction[2](35,32) = 1; + this->restriction[2](46,45) = 1; + this->restriction[2](49,50) = 1; + this->restriction[2](61,58) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](12,13) = 1; + this->restriction[3](15,14) = 1; + this->restriction[3](30,31) = 1; + this->restriction[3](34,33) = 1; + this->restriction[3](48,51) = 1; + this->restriction[3](54,53) = 1; + this->restriction[3](60,59) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](16,17) = 1; + this->restriction[4](22,23) = 1; + this->restriction[4](25,24) = 1; + this->restriction[4](36,39) = 1; + this->restriction[4](42,41) = 1; + this->restriction[4](53,54) = 1; + this->restriction[4](58,61) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](17,16) = 1; + this->restriction[5](18,19) = 1; + this->restriction[5](27,26) = 1; + this->restriction[5](37,38) = 1; + this->restriction[5](43,40) = 1; + this->restriction[5](45,46) = 1; + this->restriction[5](59,60) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](19,18) = 1; + this->restriction[6](21,20) = 1; + this->restriction[6](29,28) = 1; + this->restriction[6](39,36) = 1; + this->restriction[6](47,44) = 1; + this->restriction[6](51,48) = 1; + this->restriction[6](63,56) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](20,21) = 1; + this->restriction[7](23,22) = 1; + this->restriction[7](31,30) = 1; + this->restriction[7](38,37) = 1; + this->restriction[7](50,49) = 1; + this->restriction[7](55,52) = 1; + this->restriction[7](62,57) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](8,9) = 1; + this->restriction[0](9,1) = 1; + this->restriction[0](17,18) = 1; + this->restriction[0](18,3) = 1; + this->restriction[0](32,33) = 1; + this->restriction[0](33,4) = 1; + this->restriction[0](44,48) = 1; + this->restriction[0](45,12) = 1; + this->restriction[0](47,15) = 1; + this->restriction[0](48,2) = 1; + this->restriction[0](62,66) = 1; + this->restriction[0](63,36) = 1; + this->restriction[0](65,21) = 1; + this->restriction[0](66,5) = 1; + this->restriction[0](89,93) = 1; + this->restriction[0](90,30) = 1; + this->restriction[0](92,42) = 1; + this->restriction[0](93,7) = 1; + this->restriction[0](98,111) = 1; + this->restriction[0](99,75) = 1; + this->restriction[0](101,57) = 1; + this->restriction[0](102,24) = 1; + this->restriction[0](107,84) = 1; + this->restriction[0](108,39) = 1; + this->restriction[0](110,27) = 1; + this->restriction[0](111,6) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](9,0) = 1; + this->restriction[1](10,9) = 1; + this->restriction[1](11,12) = 1; + this->restriction[1](12,2) = 1; + this->restriction[1](35,36) = 1; + this->restriction[1](36,5) = 1; + this->restriction[1](45,18) = 1; + this->restriction[1](46,48) = 1; + this->restriction[1](48,3) = 1; + this->restriction[1](49,15) = 1; + this->restriction[1](63,33) = 1; + this->restriction[1](64,66) = 1; + this->restriction[1](66,4) = 1; + this->restriction[1](67,21) = 1; + this->restriction[1](71,75) = 1; + this->restriction[1](72,24) = 1; + this->restriction[1](74,39) = 1; + this->restriction[1](75,6) = 1; + this->restriction[1](99,93) = 1; + this->restriction[1](100,111) = 1; + this->restriction[1](102,30) = 1; + this->restriction[1](103,57) = 1; + this->restriction[1](108,42) = 1; + this->restriction[1](109,84) = 1; + this->restriction[1](111,7) = 1; + this->restriction[1](112,27) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](12,1) = 1; + this->restriction[2](13,12) = 1; + this->restriction[2](15,3) = 1; + this->restriction[2](16,15) = 1; + this->restriction[2](38,39) = 1; + this->restriction[2](39,6) = 1; + this->restriction[2](48,0) = 1; + this->restriction[2](49,9) = 1; + this->restriction[2](51,18) = 1; + this->restriction[2](52,48) = 1; + this->restriction[2](74,36) = 1; + this->restriction[2](75,5) = 1; + this->restriction[2](77,75) = 1; + this->restriction[2](78,24) = 1; + this->restriction[2](81,42) = 1; + this->restriction[2](82,84) = 1; + this->restriction[2](84,7) = 1; + this->restriction[2](85,27) = 1; + this->restriction[2](108,33) = 1; + this->restriction[2](109,66) = 1; + this->restriction[2](111,4) = 1; + this->restriction[2](112,21) = 1; + this->restriction[2](117,93) = 1; + this->restriction[2](118,111) = 1; + this->restriction[2](120,30) = 1; + this->restriction[2](121,57) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](14,15) = 1; + this->restriction[3](15,2) = 1; + this->restriction[3](18,0) = 1; + this->restriction[3](19,18) = 1; + this->restriction[3](41,42) = 1; + this->restriction[3](42,7) = 1; + this->restriction[3](47,9) = 1; + this->restriction[3](48,1) = 1; + this->restriction[3](50,48) = 1; + this->restriction[3](51,12) = 1; + this->restriction[3](80,84) = 1; + this->restriction[3](81,39) = 1; + this->restriction[3](83,27) = 1; + this->restriction[3](84,6) = 1; + this->restriction[3](92,33) = 1; + this->restriction[3](93,4) = 1; + this->restriction[3](95,93) = 1; + this->restriction[3](96,30) = 1; + this->restriction[3](107,66) = 1; + this->restriction[3](108,36) = 1; + this->restriction[3](110,21) = 1; + this->restriction[3](111,5) = 1; + this->restriction[3](116,111) = 1; + this->restriction[3](117,75) = 1; + this->restriction[3](119,57) = 1; + this->restriction[3](120,24) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](20,21) = 1; + this->restriction[4](21,5) = 1; + this->restriction[4](29,30) = 1; + this->restriction[4](30,7) = 1; + this->restriction[4](33,0) = 1; + this->restriction[4](34,33) = 1; + this->restriction[4](53,57) = 1; + this->restriction[4](54,24) = 1; + this->restriction[4](56,27) = 1; + this->restriction[4](57,6) = 1; + this->restriction[4](65,9) = 1; + this->restriction[4](66,1) = 1; + this->restriction[4](68,66) = 1; + this->restriction[4](69,36) = 1; + this->restriction[4](90,18) = 1; + this->restriction[4](91,93) = 1; + this->restriction[4](93,3) = 1; + this->restriction[4](94,42) = 1; + this->restriction[4](101,48) = 1; + this->restriction[4](102,12) = 1; + this->restriction[4](104,111) = 1; + this->restriction[4](105,75) = 1; + this->restriction[4](110,15) = 1; + this->restriction[4](111,2) = 1; + this->restriction[4](113,84) = 1; + this->restriction[4](114,39) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](21,4) = 1; + this->restriction[5](22,21) = 1; + this->restriction[5](23,24) = 1; + this->restriction[5](24,6) = 1; + this->restriction[5](36,1) = 1; + this->restriction[5](37,36) = 1; + this->restriction[5](54,30) = 1; + this->restriction[5](55,57) = 1; + this->restriction[5](57,7) = 1; + this->restriction[5](58,27) = 1; + this->restriction[5](66,0) = 1; + this->restriction[5](67,9) = 1; + this->restriction[5](69,33) = 1; + this->restriction[5](70,66) = 1; + this->restriction[5](72,12) = 1; + this->restriction[5](73,75) = 1; + this->restriction[5](75,2) = 1; + this->restriction[5](76,39) = 1; + this->restriction[5](102,18) = 1; + this->restriction[5](103,48) = 1; + this->restriction[5](105,93) = 1; + this->restriction[5](106,111) = 1; + this->restriction[5](111,3) = 1; + this->restriction[5](112,15) = 1; + this->restriction[5](114,42) = 1; + this->restriction[5](115,84) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](24,5) = 1; + this->restriction[6](25,24) = 1; + this->restriction[6](27,7) = 1; + this->restriction[6](28,27) = 1; + this->restriction[6](39,2) = 1; + this->restriction[6](40,39) = 1; + this->restriction[6](57,4) = 1; + this->restriction[6](58,21) = 1; + this->restriction[6](60,30) = 1; + this->restriction[6](61,57) = 1; + this->restriction[6](75,1) = 1; + this->restriction[6](76,36) = 1; + this->restriction[6](78,12) = 1; + this->restriction[6](79,75) = 1; + this->restriction[6](84,3) = 1; + this->restriction[6](85,15) = 1; + this->restriction[6](87,42) = 1; + this->restriction[6](88,84) = 1; + this->restriction[6](111,0) = 1; + this->restriction[6](112,9) = 1; + this->restriction[6](114,33) = 1; + this->restriction[6](115,66) = 1; + this->restriction[6](120,18) = 1; + this->restriction[6](121,48) = 1; + this->restriction[6](123,93) = 1; + this->restriction[6](124,111) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](26,27) = 1; + this->restriction[7](27,6) = 1; + this->restriction[7](30,4) = 1; + this->restriction[7](31,30) = 1; + this->restriction[7](42,3) = 1; + this->restriction[7](43,42) = 1; + this->restriction[7](56,21) = 1; + this->restriction[7](57,5) = 1; + this->restriction[7](59,57) = 1; + this->restriction[7](60,24) = 1; + this->restriction[7](83,15) = 1; + this->restriction[7](84,2) = 1; + this->restriction[7](86,84) = 1; + this->restriction[7](87,39) = 1; + this->restriction[7](93,0) = 1; + this->restriction[7](94,33) = 1; + this->restriction[7](96,18) = 1; + this->restriction[7](97,93) = 1; + this->restriction[7](110,9) = 1; + this->restriction[7](111,1) = 1; + this->restriction[7](113,66) = 1; + this->restriction[7](114,36) = 1; + this->restriction[7](119,48) = 1; + this->restriction[7](120,12) = 1; + this->restriction[7](122,111) = 1; + this->restriction[7](123,75) = 1; + break; + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; default: - Assert (false,ExcNotImplemented()); + Assert (false,ExcNotImplemented()); } // finally fill in support points @@ -642,7 +642,7 @@ double FE_Q::shape_value (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); return polynomial_space.compute_value(renumber_inverse[i], p); } @@ -653,7 +653,7 @@ FE_Q::shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_value(renumber_inverse[i], p); } @@ -665,7 +665,7 @@ Tensor<1,dim> FE_Q::shape_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); return polynomial_space.compute_grad(renumber_inverse[i], p); } @@ -677,7 +677,7 @@ FE_Q::shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad(renumber_inverse[i], p); } @@ -689,7 +689,7 @@ Tensor<2,dim> FE_Q::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); return polynomial_space.compute_grad_grad(renumber_inverse[i], p); } @@ -701,7 +701,7 @@ FE_Q::shape_grad_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); return polynomial_space.compute_grad_grad(renumber_inverse[i], p); } @@ -721,7 +721,7 @@ void FE_Q::initialize_unit_support_points () for (unsigned int i=1; iunit_support_points.resize(n); const double step = 1./degree; Point p; @@ -737,7 +737,7 @@ void FE_Q::initialize_unit_support_points () if (dim>2) p(2) = iz * step; - unit_support_points[renumber[k++]] = p; + this->unit_support_points[renumber[k++]] = p; }; }; @@ -763,7 +763,7 @@ void FE_Q::initialize_unit_face_support_points () for (unsigned int i=1; iunit_face_support_points.resize(n); const double step = 1./degree; Point p; @@ -779,7 +779,7 @@ void FE_Q::initialize_unit_face_support_points () if (codim>2) p(2) = iz * step; - unit_face_support_points[face_renumber[k++]] = p; + this->unit_face_support_points[face_renumber[k++]] = p; }; }; @@ -844,34 +844,34 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & switch (dim) { case 1: - { - const unsigned int values[GeometryInfo<1>::vertices_per_cell] - = { 0, degree }; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<1>::vertices_per_cell] + = { 0, degree }; + index = values[i]; + break; + }; case 2: - { - const unsigned int values[GeometryInfo<2>::vertices_per_cell] - = { 0, degree, n*degree+degree, n*degree }; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<2>::vertices_per_cell] + = { 0, degree, n*degree+degree, n*degree }; + index = values[i]; + break; + }; case 3: - { - const unsigned int values[GeometryInfo<3>::vertices_per_cell] - = { 0, degree, - n*n*degree + degree, n*n*degree, - n*degree, n*degree+degree, - n*n*degree + n*degree+degree, n*n*degree + n*degree}; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<3>::vertices_per_cell] + = { 0, degree, + n*n*degree + degree, n*n*degree, + n*degree, n*degree+degree, + n*n*degree + n*degree+degree, n*n*degree + n*degree}; + index = values[i]; + break; + }; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } renumber[index] = i; @@ -903,19 +903,19 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & case 100: case 200: case 202: case 300: case 302: case 304: case 306: - incr = 1; - break; - // lines in y-direction + incr = 1; + break; + // lines in y-direction case 201: case 203: case 308: case 309: case 310: case 311: - incr = n; - break; - // lines in z-direction + incr = n; + break; + // lines in z-direction case 301: case 303: case 305: case 307: - incr = n*n; - break; + incr = n*n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } switch (i+100*dim) { @@ -923,36 +923,36 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & case 100: case 200: case 203: case 300: case 303: case 308: - tensorstart = 0; - break; - // x=1 y=z=0 + tensorstart = 0; + break; + // x=1 y=z=0 case 201: case 301: case 309: - tensorstart = degree; - break; - // y=1 x=z=0 + tensorstart = degree; + break; + // y=1 x=z=0 case 202: case 304: case 307: - tensorstart = n*degree; - break; - // x=z=1 y=0 + tensorstart = n*degree; + break; + // x=z=1 y=0 case 310: - tensorstart = n*n*degree+degree; - break; - // z=1 x=y=0 + tensorstart = n*n*degree+degree; + break; + // z=1 x=y=0 case 302: case 311: - tensorstart = n*n*degree; - break; - // x=y=1 z=0 + tensorstart = n*n*degree; + break; + // x=y=1 z=0 case 305: - tensorstart = n*degree+degree; - break; - // y=z=1 x=0 + tensorstart = n*degree+degree; + break; + // y=z=1 x=0 case 306: - tensorstart = n*n*n-n; - break; + tensorstart = n*n*n-n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } for (unsigned int jx = 1; jx::lexicographic_to_hierarchic_numbering (const FiniteElementData & switch (i) { case 0: - tensorstart = 0; incx = 1; - if (dim==2) - incy = n; - else - incy = n*n; - break; + tensorstart = 0; incx = 1; + if (dim==2) + incy = n; + else + incy = n*n; + break; case 1: - tensorstart = n*degree; incx = 1; incy = n*n; - break; + tensorstart = n*degree; incx = 1; incy = n*n; + break; case 2: - tensorstart = 0; incx = 1; incy = n; - break; + tensorstart = 0; incx = 1; incy = n; + break; case 3: - tensorstart = degree; incx = n; incy = n*n; - break; + tensorstart = degree; incx = n; incy = n*n; + break; case 4: - tensorstart = n*n*degree; incx = 1; incy = n; - break; + tensorstart = n*n*degree; incx = 1; incy = n; + break; case 5: - tensorstart = 0; incx = n; incy = n*n; - break; + tensorstart = 0; incx = n; incy = n*n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } for (unsigned int jy = 1; jy::get_data (const UpdateFlags update_flags, // allocate memory if (flags & update_values) { - values.resize (dofs_per_cell); - data->shape_values.resize(dofs_per_cell, + values.resize (this->dofs_per_cell); + data->shape_values.resize(this->dofs_per_cell, std::vector(n_q_points)); } if (flags & update_gradients) { - grads.resize (dofs_per_cell); - data->shape_gradients.resize(dofs_per_cell, + grads.resize (this->dofs_per_cell); + data->shape_gradients.resize(this->dofs_per_cell, std::vector >(n_q_points)); } @@ -1142,11 +1142,11 @@ FE_Q::get_data (const UpdateFlags update_flags, values, grads, grad_grads); if (flags & update_values) - for (unsigned int k=0; kdofs_per_cell; ++k) data->shape_values[renumber[k]][i] = values[k]; if (flags & update_gradients) - for (unsigned int k=0; kdofs_per_cell; ++k) data->shape_gradients[renumber[k]][i] = grads[k]; } return data; @@ -1176,7 +1176,7 @@ FE_Q::fill_fe_values (const Mapping &mapping, const UpdateFlags flags(fe_data.current_update_flags()); - for (unsigned int k=0; kdofs_per_cell; ++k) { if (flags & update_values) for (unsigned int i=0; i::fill_fe_face_values (const Mapping &mapping, const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::fill_fe_subface_values (const Mapping &mapping const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0; kdofs_per_cell; ++k) { for (unsigned int i=0;i::has_support_on_face (const unsigned int shape_index_, const unsigned int face_index) const { - Assert (shape_index_ < dofs_per_cell, - ExcIndexRange (shape_index_, 0, dofs_per_cell)); + Assert (shape_index_ < this->dofs_per_cell, + ExcIndexRange (shape_index_, 0, this->dofs_per_cell)); Assert (face_index < GeometryInfo::faces_per_cell, ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); @@ -1333,11 +1333,11 @@ FE_Q::has_support_on_face (const unsigned int shape_index_, Assert (dim<=3, ExcNotImplemented()); const unsigned int cell_start = (dim==2) - ? first_quad_index - : first_hex_index; + ? this->first_quad_index + : this->first_hex_index; const unsigned int face_start = (dim==2) - ? first_line_index - : first_quad_index; + ? this->first_line_index + : this->first_quad_index; // Interior degrees of // freedom correspond to @@ -1352,13 +1352,13 @@ FE_Q::has_support_on_face (const unsigned int shape_index_, // equal to the face index. if (shape_index >= face_start) { - shape_index -= first_line_index; - shape_index /= dofs_per_face; + shape_index -= this->first_line_index; + shape_index /= this->dofs_per_face; return (shape_index == face_index); } // Only degrees of freedom on // a vertex are left. - shape_index /= dofs_per_vertex; + shape_index /= this->dofs_per_vertex; // Use a table to figure out // which face is neighbor to // which vertex. @@ -1372,9 +1372,9 @@ FE_Q::has_support_on_face (const unsigned int shape_index_, case 314: case 324: case 354: case 315: case 325: case 335: case 316: case 336: case 346: case 317: case 347: case 357: - return true; + return true; default: - return false; + return false; } return true; }; diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 9551529c03..69b3ba4735 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -140,7 +140,7 @@ FESystem::FESystem (const FiniteElement &fe, FiniteElement (multiply_dof_numbers(fe, n_elements), compute_restriction_is_additive_flags (fe, n_elements), compute_nonzero_components(fe, n_elements)), - base_elements(1) + base_elements(1) { base_elements[0] = ElementPair(fe.clone(), n_elements); base_elements[0].first->subscribe (); @@ -159,7 +159,7 @@ FESystem::FESystem (const FiniteElement &fe1, fe2, n2), compute_nonzero_components(fe1, n1, fe2, n2)), - base_elements(2) + base_elements(2) { base_elements[0] = ElementPair(fe1.clone(), n1); base_elements[0].first->subscribe (); @@ -186,7 +186,7 @@ FESystem::FESystem (const FiniteElement &fe1, compute_nonzero_components(fe1, n1, fe2, n2, fe3, n3)), - base_elements(3) + base_elements(3) { base_elements[0] = ElementPair(fe1.clone(), n1); base_elements[0].first->subscribe (); @@ -221,22 +221,22 @@ FESystem::clone() const switch (n_base_elements()) { case 1: - return new FESystem(base_element(0), - element_multiplicity(0)); + return new FESystem(base_element(0), + element_multiplicity(0)); case 2: - return new FESystem(base_element(0), - element_multiplicity(0), - base_element(1), - element_multiplicity(1)); + return new FESystem(base_element(0), + element_multiplicity(0), + base_element(1), + element_multiplicity(1)); case 3: - return new FESystem(base_element(0), - element_multiplicity(0), - base_element(1), - element_multiplicity(1), - base_element(2), - element_multiplicity(2)); + return new FESystem(base_element(0), + element_multiplicity(0), + base_element(1), + element_multiplicity(1), + base_element(2), + element_multiplicity(2)); default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } return 0; } @@ -248,12 +248,12 @@ double FESystem::shape_value (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (this->is_primitive(i), typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); - return (base_element(system_to_base_table[i].first.first) - .shape_value(system_to_base_table[i].second, p)); + return (base_element(this->system_to_base_table[i].first.first) + .shape_value(this->system_to_base_table[i].second, p)); } @@ -264,17 +264,17 @@ FESystem::shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (component < this->n_components(), + ExcIndexRange (component, 0, this->n_components())); // first find out to which of the // base elements this desired // component belongs, and which // component within this base // element it is - const unsigned int base = component_to_base(i).first; - const unsigned int component_in_base = component_to_base(i).second; + const unsigned int base = this->component_to_base(i).first; + const unsigned int component_in_base = this->component_to_base(i).second; // then get value from base // element. note that that will @@ -283,7 +283,7 @@ FESystem::shape_value_component (const unsigned int i, // primitive; thus, there is no // need to check this here return (base_element(base). - shape_value_component(system_to_base_table[i].second, + shape_value_component(this->system_to_base_table[i].second, p, component_in_base)); } @@ -295,12 +295,12 @@ Tensor<1,dim> FESystem::shape_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (this->is_primitive(i), typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); - return (base_element(system_to_base_table[i].first.first) - .shape_grad(system_to_base_table[i].second, p)); + return (base_element(this->system_to_base_table[i].first.first) + .shape_grad(this->system_to_base_table[i].second, p)); } @@ -311,17 +311,17 @@ FESystem::shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (component < this->n_components(), + ExcIndexRange (component, 0, this->n_components())); // first find out to which of the // base elements this desired // component belongs, and which // component within this base // element it is - const unsigned int base = component_to_base(i).first; - const unsigned int component_in_base = component_to_base(i).second; + const unsigned int base = this->component_to_base(i).first; + const unsigned int component_in_base = this->component_to_base(i).second; // then get value from base // element. note that that will @@ -330,7 +330,7 @@ FESystem::shape_grad_component (const unsigned int i, // primitive; thus, there is no // need to check this here return (base_element(base). - shape_grad_component(system_to_base_table[i].second, + shape_grad_component(this->system_to_base_table[i].second, p, component_in_base)); } @@ -342,12 +342,12 @@ Tensor<2,dim> FESystem::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (this->is_primitive(i), typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); - return (base_element(system_to_base_table[i].first.first) - .shape_grad_grad(system_to_base_table[i].second, p)); + return (base_element(this->system_to_base_table[i].first.first) + .shape_grad_grad(this->system_to_base_table[i].second, p)); } @@ -358,17 +358,17 @@ FESystem::shape_grad_grad_component (const unsigned int i, const Point &p, const unsigned int component) const { - Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (component < this->n_components(), + ExcIndexRange (component, 0, this->n_components())); // first find out to which of the // base elements this desired // component belongs, and which // component within this base // element it is - const unsigned int base = component_to_base(i).first; - const unsigned int component_in_base = component_to_base(i).second; + const unsigned int base = this->component_to_base(i).first; + const unsigned int component_in_base = this->component_to_base(i).second; // then get value from base // element. note that that will @@ -377,7 +377,7 @@ FESystem::shape_grad_grad_component (const unsigned int i, // primitive; thus, there is no // need to check this here return (base_element(base). - shape_grad_grad_component(system_to_base_table[i].second, + shape_grad_grad_component(this->system_to_base_table[i].second, p, component_in_base)); } @@ -705,7 +705,7 @@ FESystem::compute_fill (const Mapping &mapping, for (unsigned int point=0; pointcomponent_to_system_index(comp,k); if (base_flags & update_values) data.shape_values(system_index, point)= base_data.shape_values(k,point); @@ -771,8 +771,8 @@ FESystem::build_cell_tables() for (unsigned int base=0; base < n_base_elements(); ++base) for (unsigned int m = 0; m < element_multiplicity(base); ++m) for (unsigned int k=0; kcomponent_to_base_table[total_index++] = std::make_pair(base,k); + Assert (total_index == this->component_to_base_table.size(), ExcInternalError()); // Initialize index tables. @@ -806,7 +806,7 @@ FESystem::build_cell_tables() = (base_element(base).dofs_per_vertex*vertex_number + local_index); - system_to_base_table[total_index] + this->system_to_base_table[total_index] = std::make_pair (std::make_pair(base, m), index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -817,11 +817,11 @@ FESystem::build_cell_tables() = comp_start + comp_in_base; const unsigned int index_in_comp = base_element(base).system_to_component_index(index_in_base).second; - system_to_component_table[total_index] + this->system_to_component_table[total_index] = std::make_pair (comp, index_in_comp); } else - system_to_component_table[total_index] = non_primitive_index; + this->system_to_component_table[total_index] = non_primitive_index; } } @@ -844,7 +844,7 @@ FESystem::build_cell_tables() local_index + base_element(base).first_line_index); - system_to_base_table[total_index] + this->system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -853,13 +853,13 @@ FESystem::build_cell_tables() = base_element(base).system_to_component_index(index_in_base).first; const unsigned int comp = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); + const unsigned int index_in_comp + = base_element(base).system_to_component_index(index_in_base).second; + this->system_to_component_table[total_index] + = std::make_pair (comp, index_in_comp); } else - system_to_component_table[total_index] = non_primitive_index; + this->system_to_component_table[total_index] = non_primitive_index; } } @@ -882,7 +882,7 @@ FESystem::build_cell_tables() local_index + base_element(base).first_quad_index); - system_to_base_table[total_index] + this->system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -891,13 +891,13 @@ FESystem::build_cell_tables() = base_element(base).system_to_component_index(index_in_base).first; const unsigned int comp = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); + const unsigned int index_in_comp + = base_element(base).system_to_component_index(index_in_base).second; + this->system_to_component_table[total_index] + = std::make_pair (comp, index_in_comp); } else - system_to_component_table[total_index] = non_primitive_index; + this->system_to_component_table[total_index] = non_primitive_index; } } @@ -920,7 +920,7 @@ FESystem::build_cell_tables() local_index + base_element(base).first_hex_index); - system_to_base_table[total_index] + this->system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -929,13 +929,13 @@ FESystem::build_cell_tables() = base_element(base).system_to_component_index(index_in_base).first; const unsigned int comp = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); + const unsigned int index_in_comp + = base_element(base).system_to_component_index(index_in_base).second; + this->system_to_component_table[total_index] + = std::make_pair (comp, index_in_comp); } else - system_to_component_table[total_index] = non_primitive_index; + this->system_to_component_table[total_index] = non_primitive_index; } } @@ -943,24 +943,24 @@ FESystem::build_cell_tables() // components to linear // index. Fortunately, this is the // inverse of what we just did. - std::vector dofs_per_component (n_components(), 0); - for (unsigned int sys=0; sys dofs_per_component (this->n_components(), 0); + for (unsigned int sys=0; sysdofs_per_cell; ++sys) + ++dofs_per_component[this->system_to_component_index(sys).first]; + for (unsigned int component=0; componentn_components(); ++component) + this->component_to_system_table[component].resize(dofs_per_component[component]); // then go the reverse way to fill the array - for (unsigned int sys=0; sysdofs_per_cell; ++sys) { const unsigned int - comp = system_to_component_index(sys).first, - index_in_comp = system_to_component_index(sys).second; + comp = this->system_to_component_index(sys).first, + index_in_comp = this->system_to_component_index(sys).second; - Assert (comp < component_to_system_table.size(), + Assert (comp < this->component_to_system_table.size(), ExcInternalError()); - Assert (index_in_comp < component_to_system_table[comp].size(), + Assert (index_in_comp < this->component_to_system_table[comp].size(), ExcInternalError()); - component_to_system_table[comp][index_in_comp] = sys; + this->component_to_system_table[comp][index_in_comp] = sys; }; }; @@ -1029,7 +1029,7 @@ FESystem::build_face_tables() = (base_element(base).dofs_per_vertex*vertex_number + local_index); - face_system_to_base_table[total_index] + this->face_system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), face_index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -1040,11 +1040,11 @@ FESystem::build_face_tables() = comp_start + comp_in_base; const unsigned int face_index_in_comp = base_element(base).face_system_to_component_index(face_index_in_base).second; - face_system_to_component_table[total_index] + this->face_system_to_component_table[total_index] = std::make_pair (comp, face_index_in_comp); } else - face_system_to_component_table[total_index] = non_primitive_index; + this->face_system_to_component_table[total_index] = non_primitive_index; } } @@ -1075,7 +1075,7 @@ FESystem::build_face_tables() base_element(base).dofs_per_line * line_number + local_index); - face_system_to_base_table[total_index] + this->face_system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), face_index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -1086,11 +1086,11 @@ FESystem::build_face_tables() = comp_start + comp_in_base; const unsigned int face_index_in_comp = base_element(base).face_system_to_component_index(face_index_in_base).second; - face_system_to_component_table[total_index] + this->face_system_to_component_table[total_index] = std::make_pair (comp, face_index_in_comp); } else - face_system_to_component_table[total_index] = non_primitive_index; + this->face_system_to_component_table[total_index] = non_primitive_index; } } @@ -1121,7 +1121,7 @@ FESystem::build_face_tables() base_element(base).dofs_per_quad * quad_number + local_index); - face_system_to_base_table[total_index] + this->face_system_to_base_table[total_index] = std::make_pair (std::make_pair(base,m), face_index_in_base); if (base_element(base).is_primitive(index_in_base)) @@ -1132,39 +1132,39 @@ FESystem::build_face_tables() = comp_start + comp_in_base; const unsigned int face_index_in_comp = base_element(base).face_system_to_component_index(face_index_in_base).second; - face_system_to_component_table[total_index] + this->face_system_to_component_table[total_index] = std::make_pair (comp, face_index_in_comp); } else - face_system_to_component_table[total_index] = non_primitive_index; + this->face_system_to_component_table[total_index] = non_primitive_index; } } - Assert (total_index == dofs_per_face, ExcInternalError()); - Assert (total_index == face_system_to_component_table.size(), + Assert (total_index == this->dofs_per_face, ExcInternalError()); + Assert (total_index == this->face_system_to_component_table.size(), ExcInternalError()); - Assert (total_index == face_system_to_base_table.size(), + Assert (total_index == this->face_system_to_base_table.size(), ExcInternalError()); // finally, initialize reverse mapping - std::vector dofs_per_component (n_components(), 0); - for (unsigned int sys=0; sys dofs_per_component (this->n_components(), 0); + for (unsigned int sys=0; sysdofs_per_face; ++sys) + ++dofs_per_component[this->face_system_to_component_index(sys).first]; + for (unsigned int component=0; componentn_components(); ++component) + this->face_component_to_system_table[component].resize(dofs_per_component[component]); // then go the reverse way to fill // the array - for (unsigned int sys=0; sysdofs_per_face; ++sys) { const unsigned int - comp = face_system_to_component_index(sys).first, - index_in_comp = face_system_to_component_index(sys).second; + comp = this->face_system_to_component_index(sys).first, + index_in_comp = this->face_system_to_component_index(sys).second; - Assert (comp < face_component_to_system_table.size(), + Assert (comp < this->face_component_to_system_table.size(), ExcInternalError()); - Assert (index_in_comp < face_component_to_system_table[comp].size(), + Assert (index_in_comp < this->face_component_to_system_table[comp].size(), ExcInternalError()); - face_component_to_system_table[comp][index_in_comp] = sys; + this->face_component_to_system_table[comp][index_in_comp] = sys; }; }; @@ -1191,8 +1191,8 @@ void FESystem::build_interface_constraints () // function in conjunction with the // numbers // first_{line,quad,...}_index - for (unsigned int n=0; ninterface_constraints.n(); ++n) + for (unsigned int m=0; minterface_constraints.m(); ++m) { // for the pair (n,m) find // out which base element @@ -1210,7 +1210,7 @@ void FESystem::build_interface_constraints () // base element), second is // index within this instance const std::pair, unsigned int> n_index - = face_system_to_base_table[n]; + = this->face_system_to_base_table[n]; // likewise for the m // index. this is more @@ -1222,178 +1222,178 @@ void FESystem::build_interface_constraints () switch (dim) { case 1: - { - // we should never get here! - // (in 1d, the constraints matrix - // should be of size zero) - Assert (false, ExcInternalError()); - break; - }; + { + // we should never get here! + // (in 1d, the constraints matrix + // should be of size zero) + Assert (false, ExcInternalError()); + break; + }; case 2: - { - // the indices m=0..d_v-1 are - // from the center vertex. - // their order is the same - // as for the first vertex - // of the whole cell, so we - // can use the - // system_to_base_table - // variable (using the - // face_s_t_base_t function would - // yield the same) - if (m < dofs_per_vertex) - m_index = system_to_base_table[m]; - else - // then come the two sets of - // line indices + { + // the indices m=0..d_v-1 are + // from the center vertex. + // their order is the same + // as for the first vertex + // of the whole cell, so we + // can use the + // system_to_base_table + // variable (using the + // face_s_t_base_t function would + // yield the same) + if (m < this->dofs_per_vertex) + m_index = this->system_to_base_table[m]; + else + // then come the two sets of + // line indices + { + const unsigned int index_in_line + = (m-this->dofs_per_vertex) % this->dofs_per_line; + const unsigned int sub_line + = (m-this->dofs_per_vertex) / this->dofs_per_line; + Assert (sub_line < 2, ExcInternalError()); + + // from this + // information, try + // to get base + // element and + // instance of base + // element. we do + // so by + // constructing the + // corresponding + // face index of m + // in the present + // element, then + // use + // face_system_to_base_table + const unsigned int tmp1 = 2*this->dofs_per_vertex+index_in_line; + m_index.first = this->face_system_to_base_table[tmp1].first; + + // what we are + // still missing is + // the index of m + // within the base + // elements + // interface_constraints + // table + // + // here, the second + // value of + // face_system_to_base_table + // can help: it + // denotes the face + // index of that + // shape function + // within the base + // element. since + // we know that it + // is a line dof, + // we can construct + // the rest: tmp2 + // will denote the + // index of this + // shape function + // among the line + // shape functions: + Assert (this->face_system_to_base_table[tmp1].second >= + 2*base_element(m_index.first.first).dofs_per_vertex, + ExcInternalError()); + const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second - + 2*base_element(m_index.first.first).dofs_per_vertex; + Assert (tmp2 < base_element(m_index.first.first).dofs_per_line, + ExcInternalError()); + m_index.second = base_element(m_index.first.first).dofs_per_vertex + + base_element(m_index.first.first).dofs_per_line*sub_line + + tmp2; + }; + break; + }; + + case 3: + { + // same way as above, + // although a little + // more complicated... + + // the indices + // m=0..5*d_v-1 are + // from the center and + // the four subline + // vertices. their + // order is the same as + // for the first vertex + // of the whole cell, + // so we can use the + // simple arithmetic + if (m < 5*this->dofs_per_vertex) + m_index = this->system_to_base_table[m]; + else + // then come the 12 sets of + // line indices + if (m < 5*this->dofs_per_vertex + 12*this->dofs_per_line) { + // for the + // meaning of all + // this, see the + // 2d part const unsigned int index_in_line - = (m-dofs_per_vertex) % dofs_per_line; + = (m-5*this->dofs_per_vertex) % this->dofs_per_line; const unsigned int sub_line - = (m-dofs_per_vertex) / dofs_per_line; - Assert (sub_line < 2, ExcInternalError()); - - // from this - // information, try - // to get base - // element and - // instance of base - // element. we do - // so by - // constructing the - // corresponding - // face index of m - // in the present - // element, then - // use - // face_system_to_base_table - const unsigned int tmp1 = 2*dofs_per_vertex+index_in_line; - m_index.first = face_system_to_base_table[tmp1].first; - - // what we are - // still missing is - // the index of m - // within the base - // elements - // interface_constraints - // table - // - // here, the second - // value of - // face_system_to_base_table - // can help: it - // denotes the face - // index of that - // shape function - // within the base - // element. since - // we know that it - // is a line dof, - // we can construct - // the rest: tmp2 - // will denote the - // index of this - // shape function - // among the line - // shape functions: - Assert (face_system_to_base_table[tmp1].second >= - 2*base_element(m_index.first.first).dofs_per_vertex, + = (m-5*this->dofs_per_vertex) / this->dofs_per_line; + Assert (sub_line < 12, ExcInternalError()); + + const unsigned int tmp1 = 4*this->dofs_per_vertex+index_in_line; + m_index.first = this->face_system_to_base_table[tmp1].first; + + Assert (this->face_system_to_base_table[tmp1].second >= + 4*base_element(m_index.first.first).dofs_per_vertex, ExcInternalError()); - const unsigned int tmp2 = face_system_to_base_table[tmp1].second - - 2*base_element(m_index.first.first).dofs_per_vertex; + const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second - + 4*base_element(m_index.first.first).dofs_per_vertex; Assert (tmp2 < base_element(m_index.first.first).dofs_per_line, ExcInternalError()); - m_index.second = base_element(m_index.first.first).dofs_per_vertex + + m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex + base_element(m_index.first.first).dofs_per_line*sub_line + tmp2; - }; - break; - }; - - case 3: - { - // same way as above, - // although a little - // more complicated... - - // the indices - // m=0..5*d_v-1 are - // from the center and - // the four subline - // vertices. their - // order is the same as - // for the first vertex - // of the whole cell, - // so we can use the - // simple arithmetic - if (m < 5*dofs_per_vertex) - m_index = system_to_base_table[m]; + } else - // then come the 12 sets of - // line indices - if (m < 5*dofs_per_vertex + 12*dofs_per_line) - { - // for the - // meaning of all - // this, see the - // 2d part - const unsigned int index_in_line - = (m-5*dofs_per_vertex) % dofs_per_line; - const unsigned int sub_line - = (m-5*dofs_per_vertex) / dofs_per_line; - Assert (sub_line < 12, ExcInternalError()); - - const unsigned int tmp1 = 4*dofs_per_vertex+index_in_line; - m_index.first = face_system_to_base_table[tmp1].first; - - Assert (face_system_to_base_table[tmp1].second >= - 4*base_element(m_index.first.first).dofs_per_vertex, - ExcInternalError()); - const unsigned int tmp2 = face_system_to_base_table[tmp1].second - - 4*base_element(m_index.first.first).dofs_per_vertex; - Assert (tmp2 < base_element(m_index.first.first).dofs_per_line, - ExcInternalError()); - m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex + - base_element(m_index.first.first).dofs_per_line*sub_line + - tmp2; - } - else - // on one of the four sub-quads - { - // for the - // meaning of all - // this, see the - // 2d part - const unsigned int index_in_quad - = (m-5*dofs_per_vertex-12*dofs_per_line) % dofs_per_line; - const unsigned int sub_quad - = (m-5*dofs_per_vertex-12*dofs_per_line) / dofs_per_line; - Assert (sub_quad < 4, ExcInternalError()); - - const unsigned int tmp1 = 4*dofs_per_vertex+4*dofs_per_line+index_in_quad; - m_index.first = face_system_to_base_table[tmp1].first; - - Assert (face_system_to_base_table[tmp1].second >= - 4*base_element(m_index.first.first).dofs_per_vertex + - 4*base_element(m_index.first.first).dofs_per_line, - ExcInternalError()); - const unsigned int tmp2 = face_system_to_base_table[tmp1].second - - 4*base_element(m_index.first.first).dofs_per_vertex - - 4*base_element(m_index.first.first).dofs_per_line; - Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad, - ExcInternalError()); - m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex + - 12*base_element(m_index.first.first).dofs_per_line + - base_element(m_index.first.first).dofs_per_quad*sub_quad + - tmp2; - }; + // on one of the four sub-quads + { + // for the + // meaning of all + // this, see the + // 2d part + const unsigned int index_in_quad + = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) % this->dofs_per_line; + const unsigned int sub_quad + = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) / this->dofs_per_line; + Assert (sub_quad < 4, ExcInternalError()); + + const unsigned int tmp1 = 4*this->dofs_per_vertex+4*this->dofs_per_line+index_in_quad; + m_index.first = this->face_system_to_base_table[tmp1].first; + + Assert (this->face_system_to_base_table[tmp1].second >= + 4*base_element(m_index.first.first).dofs_per_vertex + + 4*base_element(m_index.first.first).dofs_per_line, + ExcInternalError()); + const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second - + 4*base_element(m_index.first.first).dofs_per_vertex - + 4*base_element(m_index.first.first).dofs_per_line; + Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad, + ExcInternalError()); + m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex + + 12*base_element(m_index.first.first).dofs_per_line + + base_element(m_index.first.first).dofs_per_quad*sub_quad + + tmp2; + }; - break; - }; + break; + }; default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); }; // now that we gathered all @@ -1405,7 +1405,7 @@ void FESystem::build_interface_constraints () // definitely will be no // coupling if (n_index.first == m_index.first) - interface_constraints(m,n) + this->interface_constraints(m,n) = (base_element(n_index.first.first).constraints()(m_index.second, n_index.second)); }; @@ -1437,10 +1437,10 @@ void FESystem::initialize () // element as well if (!do_restriction) for (unsigned int i=0;i::children_per_cell;++i) - restriction[i].reinit(0,0); + this->restriction[i].reinit(0,0); if (!do_prolongation) for (unsigned int i=0;i::children_per_cell;++i) - prolongation[i].reinit(0,0); + this->prolongation[i].reinit(0,0); // distribute the matrices of the // base finite elements to the @@ -1456,8 +1456,8 @@ void FESystem::initialize () // couple. only DoFs that belong to // the same instance of a base // element may couple - for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) { // first find out to which // base element indices i and @@ -1471,19 +1471,19 @@ void FESystem::initialize () // element, then they cannot // couple, so go on with the // next index - if (system_to_base_table[i].first != - system_to_base_table[j].first) + if (this->system_to_base_table[i].first != + this->system_to_base_table[j].first) continue; // so get the common base // element and the indices // therein: const unsigned int - base = system_to_base_table[i].first.first; + base = this->system_to_base_table[i].first.first; const unsigned int - base_index_i = system_to_base_table[i].second, - base_index_j = system_to_base_table[j].second; + base_index_i = this->system_to_base_table[i].second, + base_index_j = this->system_to_base_table[j].second; // if we are sure that DoFs i // and j may couple, then @@ -1492,11 +1492,11 @@ void FESystem::initialize () for (unsigned int child=0; child::children_per_cell; ++child) { if (do_restriction) - restriction[child] (i,j) + this->restriction[child] (i,j) = base_element(base).restrict(child)(base_index_i,base_index_j); if (do_prolongation) - prolongation[child] (i,j) + this->prolongation[child] (i,j) = base_element(base).prolongate(child)(base_index_i,base_index_j); }; }; @@ -1536,35 +1536,35 @@ void FESystem:: initialize_unit_support_points () { - // if one of the base elements - // has no support points, then - // it makes no sense to define - // support points for the - // composed element, so return - // an empty array to - // demonstrate that - // fact + // if one of the base elements + // has no support points, then + // it makes no sense to define + // support points for the + // composed element, so return + // an empty array to + // demonstrate that + // fact for (unsigned int base_el=0; base_elunit_support_points.resize(0); return; }; // generate unit support points // from unit support points of sub // elements - unit_support_points.resize(dofs_per_cell); + this->unit_support_points.resize(this->dofs_per_cell); - for (unsigned int i=0; idofs_per_cell; ++i) { const unsigned int - base = system_to_base_table[i].first.first, - base_index = system_to_base_table[i].second; + base = this->system_to_base_table[i].first.first, + base_index = this->system_to_base_table[i].second; Assert (baseunit_support_points[i] = base_element(base).unit_support_points[base_index]; }; } @@ -1599,7 +1599,7 @@ initialize_unit_face_support_points () for (unsigned int base_el=0; base_elunit_face_support_points.resize(0); return; }; @@ -1607,17 +1607,17 @@ initialize_unit_face_support_points () // generate unit face support points // from unit support points of sub // elements - unit_face_support_points.resize(dofs_per_face); + this->unit_face_support_points.resize(this->dofs_per_face); - for (unsigned int i=0; idofs_per_face; ++i) { - const unsigned int base_i = face_system_to_base_table[i].first.first; - const unsigned int index_in_base = face_system_to_base_table[i].second; + const unsigned int base_i = this->face_system_to_base_table[i].first.first; + const unsigned int index_in_base = this->face_system_to_base_table[i].second; Assert (index_in_base < base_element(base_i).unit_face_support_points.size(), ExcInternalError()); - unit_face_support_points[i] + this->unit_face_support_points[i] = base_element(base_i).unit_face_support_points[index_in_base]; }; } @@ -2002,8 +2002,8 @@ FESystem::has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const { const std::pair component - = system_to_component_index(shape_index); - const unsigned int base = component_to_base(component.first).first; + = this->system_to_component_index(shape_index); + const unsigned int base = this->component_to_base(component.first).first; return base_element(base).has_support_on_face(component.second, face_index); } @@ -2015,12 +2015,12 @@ template unsigned int FESystem::memory_consumption () const { - // neglect size of data stored in - // @p{base_elements} due to some - // problems with teh - // compiler. should be neglectable - // after all, considering the size - // of the data of the subelements + // neglect size of data stored in + // @p{base_elements} due to some + // problems with teh + // compiler. should be neglectable + // after all, considering the size + // of the data of the subelements unsigned int mem = (FiniteElement::memory_consumption () + sizeof (base_elements)); for (unsigned int i=0; i::initialize (const unsigned int n_quadrature_points, const FiniteElement &fe, const UpdateFlags flags) { - update_flags = flags; + this->update_flags = flags; // initialize the table mapping // from shape function number to @@ -48,11 +48,11 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, // the total number of non-zero // components accumulated over all // shape functions - shape_function_to_row_table.resize (fe.dofs_per_cell); + this->shape_function_to_row_table.resize (fe.dofs_per_cell); unsigned int row = 0; for (unsigned int i=0; ishape_function_to_row_table[i] = row; row += fe.n_nonzero_components (i); }; @@ -65,33 +65,33 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, // that we will need to their // correct size if (flags & update_values) - shape_values.reinit(n_nonzero_shape_components, n_quadrature_points); + this->shape_values.reinit(n_nonzero_shape_components, n_quadrature_points); if (flags & update_gradients) { - shape_gradients.resize(n_nonzero_shape_components); + this->shape_gradients.resize(n_nonzero_shape_components); for (unsigned int i=0; ishape_gradients[i].resize(n_quadrature_points); } if (flags & update_second_derivatives) { - shape_2nd_derivatives.resize(n_nonzero_shape_components); + this->shape_2nd_derivatives.resize(n_nonzero_shape_components); for (unsigned int i=0; ishape_2nd_derivatives[i].resize(n_quadrature_points); } if (flags & update_q_points) - quadrature_points.resize(n_quadrature_points); + this->quadrature_points.resize(n_quadrature_points); if (flags & update_JxW_values) - JxW_values.resize(n_quadrature_points); + this->JxW_values.resize(n_quadrature_points); if (flags & update_boundary_forms) - boundary_forms.resize(n_quadrature_points); + this->boundary_forms.resize(n_quadrature_points); if (flags & update_normal_vectors) - normal_vectors.resize(n_quadrature_points); + this->normal_vectors.resize(n_quadrature_points); } @@ -112,7 +112,7 @@ FEValuesBase::FEValuesBase (const unsigned int n_q_points, mapping_data(0), fe_data(0) { - update_flags = flags; + this->update_flags = flags; } @@ -147,7 +147,7 @@ template void FEValuesBase::get_function_values (const InputVector &fe_function, typename std::vector &values) const { - Assert (update_flags & update_values, ExcAccessToUninitializedField()); + Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); Assert (fe->n_components() == 1, ExcWrongNoOfComponents()); Assert (values.size() == n_quadrature_points, @@ -175,7 +175,7 @@ void FEValuesBase::get_function_values (const InputVector &fe_fu for (unsigned int point=0; pointshape_value(shape_func, point)); }; @@ -192,7 +192,7 @@ void FEValuesBase::get_function_values (const InputVector for (unsigned i=0;iupdate_flags & update_values, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -231,8 +231,8 @@ template const typename std::vector > & FEValuesBase::get_quadrature_points () const { - Assert (update_flags & update_q_points, ExcAccessToUninitializedField()); - return quadrature_points; + Assert (this->update_flags & update_q_points, ExcAccessToUninitializedField()); + return this->quadrature_points; }; @@ -241,8 +241,8 @@ template const std::vector & FEValuesBase::get_JxW_values () const { - Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField()); - return JxW_values; + Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField()); + return this->JxW_values; } @@ -254,7 +254,7 @@ FEValuesBase:: get_function_grads (const InputVector &fe_function, typename std::vector > &gradients) const { - Assert (update_flags & update_gradients, ExcAccessToUninitializedField()); + Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); Assert (fe->n_components() == 1, ExcWrongNoOfComponents()); @@ -283,7 +283,7 @@ get_function_grads (const InputVector &fe_function, for (unsigned int point=0; point tmp = shape_grad(shape_func,point); + Tensor<1,dim> tmp = this->shape_grad(shape_func,point); tmp *= dof_values(shape_func); gradients[point] += tmp; }; @@ -305,7 +305,7 @@ get_function_grads (const InputVector for (unsigned i=0; iupdate_flags & update_gradients, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -331,7 +331,7 @@ get_function_grads (const InputVector for (unsigned int shape_func=0; shape_funcis_primitive (shape_func)) { - Tensor<1,dim> tmp = shape_grad(shape_func,point); + Tensor<1,dim> tmp = this->shape_grad(shape_func,point); tmp *= dof_values(shape_func); gradients[point][fe->system_to_component_index(shape_func).first] += tmp; @@ -339,7 +339,7 @@ get_function_grads (const InputVector else for (unsigned int c=0; c tmp = shape_grad_component(shape_func,point,c); + Tensor<1,dim> tmp = this->shape_grad_component(shape_func,point,c); tmp *= dof_values(shape_func); gradients[point][c] += tmp; }; @@ -358,7 +358,7 @@ get_function_2nd_derivatives (const InputVector &fe_function, ExcWrongNoOfComponents()); Assert (second_derivatives.size() == n_quadrature_points, ExcWrongVectorSize(second_derivatives.size(), n_quadrature_points)); - Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField()); + Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -382,7 +382,7 @@ get_function_2nd_derivatives (const InputVector &fe_function, for (unsigned int point=0; point tmp = shape_2nd_derivative(shape_func,point); + Tensor<2,dim> tmp = this->shape_2nd_derivative(shape_func,point); tmp *= dof_values(shape_func); second_derivatives[point] += tmp; }; @@ -404,7 +404,7 @@ get_function_2nd_derivatives (const InputVector for (unsigned i=0;iupdate_flags & update_second_derivatives, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -434,7 +434,7 @@ get_function_2nd_derivatives (const InputVector else for (unsigned int c=0; c tmp = shape_2nd_derivative_component(shape_func,point,c); + Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c); tmp *= dof_values(shape_func); second_derivs[point][c] += tmp; }; @@ -446,10 +446,10 @@ template const Point & FEValuesBase::quadrature_point (const unsigned int i) const { - Assert (update_flags & update_q_points, ExcAccessToUninitializedField()); - Assert (iupdate_flags & update_q_points, ExcAccessToUninitializedField()); + Assert (iquadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size())); - return quadrature_points[i]; + return this->quadrature_points[i]; }; @@ -458,10 +458,10 @@ FEValuesBase::quadrature_point (const unsigned int i) const template double FEValuesBase::JxW (const unsigned int i) const { - Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField()); - Assert (iupdate_flags & update_JxW_values, ExcAccessToUninitializedField()); + Assert (iJxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size())); - return JxW_values[i]; + return this->JxW_values[i]; }; @@ -470,14 +470,14 @@ template unsigned int FEValuesBase::memory_consumption () const { - return (MemoryConsumption::memory_consumption (shape_values) + - MemoryConsumption::memory_consumption (shape_gradients) + - MemoryConsumption::memory_consumption (shape_2nd_derivatives) + - MemoryConsumption::memory_consumption (JxW_values) + - MemoryConsumption::memory_consumption (quadrature_points) + - MemoryConsumption::memory_consumption (normal_vectors) + - MemoryConsumption::memory_consumption (boundary_forms) + - sizeof(update_flags) + + return (MemoryConsumption::memory_consumption (this->shape_values) + + MemoryConsumption::memory_consumption (this->shape_gradients) + + MemoryConsumption::memory_consumption (this->shape_2nd_derivatives) + + MemoryConsumption::memory_consumption (this->JxW_values) + + MemoryConsumption::memory_consumption (this->quadrature_points) + + MemoryConsumption::memory_consumption (this->normal_vectors) + + MemoryConsumption::memory_consumption (this->boundary_forms) + + sizeof(this->update_flags) + MemoryConsumption::memory_consumption (present_cell) + MemoryConsumption::memory_consumption (n_quadrature_points) + MemoryConsumption::memory_consumption (dofs_per_cell) + @@ -487,7 +487,7 @@ FEValuesBase::memory_consumption () const MemoryConsumption::memory_consumption (*mapping_data) + MemoryConsumption::memory_consumption (fe_data) + MemoryConsumption::memory_consumption (*fe_data) + - MemoryConsumption::memory_consumption (shape_function_to_row_table)); + MemoryConsumption::memory_consumption (this->shape_function_to_row_table)); }; @@ -538,7 +538,7 @@ FEValues::FEValues (const Mapping &mapping, update_default, mapping, fe), - quadrature (q) + quadrature (q) { initialize (update_flags); }; @@ -553,9 +553,9 @@ FEValues::FEValues (const FiniteElement &fe, FEValuesBase (q.n_quadrature_points, fe.dofs_per_cell, update_default, - get_default_mapping(), + this->get_default_mapping(), fe), - quadrature (q) + quadrature (q) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); initialize (update_flags); @@ -572,17 +572,17 @@ FEValues::initialize (const UpdateFlags update_flags) Assert ((update_flags & update_normal_vectors) == false, typename FEValuesBase::ExcInvalidUpdateFlag()); - const UpdateFlags flags = compute_update_flags (update_flags); + const UpdateFlags flags = this->compute_update_flags (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); + this->mapping_data = this->mapping->get_data(flags, quadrature); + this->fe_data = this->fe->get_data(flags, *this->mapping, quadrature); // set up objects within this class - FEValuesData::initialize (n_quadrature_points, *fe, flags); + FEValuesData::initialize (this->n_quadrature_points, *this->fe, flags); }; @@ -594,24 +594,24 @@ void FEValues::reinit (const typename DoFHandler::cell_iterator &cell) // passed to the constructor and // used by the DoFHandler used by // this cell, are the same - Assert (static_cast&>(*fe) == + Assert (static_cast&>(*this->fe) == static_cast&>(cell->get_dof_handler().get_fe()), typename FEValuesBase::ExcFEDontMatch()); - present_cell = cell; + this->present_cell = cell; - get_mapping().fill_fe_values(cell, - quadrature, - *mapping_data, - quadrature_points, - JxW_values); + this->get_mapping().fill_fe_values(cell, + quadrature, + *this->mapping_data, + this->quadrature_points, + this->JxW_values); - get_fe().fill_fe_values(get_mapping(), - cell, - quadrature, - *mapping_data, - *fe_data, - *this); + this->get_fe().fill_fe_values(this->get_mapping(), + cell, + quadrature, + *this->mapping_data, + *this->fe_data, + *this); } @@ -641,7 +641,7 @@ FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, update_default, mapping, fe), - quadrature(quadrature) + quadrature(quadrature) {}; @@ -650,9 +650,9 @@ template const typename std::vector > & FEFaceValuesBase::get_normal_vectors () const { - Assert (update_flags & update_normal_vectors, + Assert (this->update_flags & update_normal_vectors, typename FEValuesBase::ExcAccessToUninitializedField()); - return normal_vectors; + return this->normal_vectors; }; @@ -661,9 +661,9 @@ template const typename std::vector > & FEFaceValuesBase::get_boundary_forms () const { - Assert (update_flags & update_boundary_forms, + Assert (this->update_flags & update_boundary_forms, FEValuesBase::ExcAccessToUninitializedField()); - return boundary_forms; + return this->boundary_forms; }; @@ -706,7 +706,7 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, FEFaceValuesBase (quadrature.n_quadrature_points, fe.dofs_per_cell, update_flags, - get_default_mapping(), + this->get_default_mapping(), fe, quadrature) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -719,17 +719,17 @@ template void FEFaceValues::initialize (const UpdateFlags update_flags) { - const UpdateFlags flags = compute_update_flags (update_flags); + const UpdateFlags flags = this->compute_update_flags (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); + this->mapping_data = this->mapping->get_face_data(flags, this->quadrature); + this->fe_data = this->fe->get_face_data(flags, *this->mapping, this->quadrature); // set up objects within this class - FEValuesData::initialize(n_quadrature_points, *fe, flags); + FEValuesData::initialize(this->n_quadrature_points, *this->fe, flags); }; @@ -742,27 +742,27 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c // passed to the constructor and // used by the DoFHandler used by // this cell, are the same - Assert (static_cast&>(*fe) == + Assert (static_cast&>(*this->fe) == static_cast&>(cell->get_dof_handler().get_fe()), typename FEValuesBase::ExcFEDontMatch()); - present_cell = cell; - present_face = cell->face(face_no); + this->present_cell = cell; + this->present_face = cell->face(face_no); - get_mapping().fill_fe_face_values(cell, face_no, - quadrature, - *mapping_data, - quadrature_points, - JxW_values, - boundary_forms, - normal_vectors); + this->get_mapping().fill_fe_face_values(cell, face_no, + this->quadrature, + *this->mapping_data, + this->quadrature_points, + this->JxW_values, + this->boundary_forms, + this->normal_vectors); - get_fe().fill_fe_face_values(get_mapping(), - cell, face_no, - quadrature, - *mapping_data, - *fe_data, - *this); + this->get_fe().fill_fe_face_values(this->get_mapping(), + cell, face_no, + this->quadrature, + *this->mapping_data, + *this->fe_data, + *this); }; @@ -794,7 +794,7 @@ FESubfaceValues::FESubfaceValues (const FiniteElement &fe, FEFaceValuesBase (quadrature.n_quadrature_points, fe.dofs_per_cell, update_flags, - get_default_mapping(), + this->get_default_mapping(), fe, quadrature) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -807,17 +807,19 @@ template void FESubfaceValues::initialize (const UpdateFlags update_flags) { - const UpdateFlags flags = compute_update_flags (update_flags); + const UpdateFlags flags = this->compute_update_flags (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_subface_data(flags, quadrature); - fe_data = fe->get_subface_data(flags, *mapping, quadrature); + this->mapping_data = this->mapping->get_subface_data(flags, this->quadrature); + this->fe_data = this->fe->get_subface_data(flags, + *this->mapping, + this->quadrature); // set up objects within this class - FEValuesData::initialize(n_quadrature_points, *fe, flags); + FEValuesData::initialize(this->n_quadrature_points, *this->fe, flags); }; @@ -831,7 +833,7 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator // passed to the constructor and // used by the DoFHandler used by // this cell, are the same - Assert (static_cast&>(*fe) == + Assert (static_cast&>(*this->fe) == static_cast&>(cell->get_dof_handler().get_fe()), typename FEValuesBase::ExcFEDontMatch()); Assert (face_no < GeometryInfo::faces_per_cell, @@ -839,23 +841,23 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator Assert (subface_no < GeometryInfo::subfaces_per_face, ExcIndexRange (subface_no, 0, GeometryInfo::subfaces_per_face)); - present_cell = cell; - present_face = cell->face(face_no); - - get_mapping().fill_fe_subface_values(cell, face_no, subface_no, - quadrature, - *mapping_data, - quadrature_points, - JxW_values, - boundary_forms, - normal_vectors); + this->present_cell = cell; + this->present_face = cell->face(face_no); + + this->get_mapping().fill_fe_subface_values(cell, face_no, subface_no, + this->quadrature, + *this->mapping_data, + this->quadrature_points, + this->JxW_values, + this->boundary_forms, + this->normal_vectors); - get_fe().fill_fe_subface_values(get_mapping(), - cell, face_no, subface_no, - quadrature, - *mapping_data, - *fe_data, - *this); + this->get_fe().fill_fe_subface_values(this->get_mapping(), + cell, face_no, subface_no, + this->quadrature, + *this->mapping_data, + *this->fe_data, + *this); }; diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index a8609cc495..b1c4e8facf 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -81,7 +81,7 @@ typename Mapping::InternalDataBase * MappingCartesian::get_data (const UpdateFlags update_flags, const Quadrature &q) const { - // Assert (flags & update_normal_vectors == 0, ExcNotImplemented()); + // Assert (flags & update_normal_vectors == 0, ExcNotImplemented()); InternalData* data = new InternalData (q); @@ -231,38 +231,38 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat { // 2D case 200: - n (1) = -1.; - break; + n (1) = -1.; + break; case 201: - n (0) = 1.; - break; + n (0) = 1.; + break; case 202: - n (1) = 1.; - break; + n (1) = 1.; + break; case 203: - n (0) = -1.; - break; - // 3D + n (0) = -1.; + break; + // 3D case 300: - n (1) = -1.; - break; + n (1) = -1.; + break; case 301: - n (1) = 1.; - break; + n (1) = 1.; + break; case 302: - n (2) = -1.; - break; + n (2) = -1.; + break; case 303: - n (0) = 1.; - break; + n (0) = 1.; + break; case 304: - n (2) = 1.; - break; + n (2) = 1.; + break; case 305: - n (0) = -1.; - break; + n (0) = -1.; + break; default: - Assert (false, ExcInternalError()); + Assert (false, ExcInternalError()); } // furthermore, all normal // vectors on a face are equal @@ -339,7 +339,7 @@ MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell // since the jacobian is diagonal double J = 1.; for (unsigned int d=0;dnormal_directions[face_no]/2)) J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) @@ -359,14 +359,14 @@ MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell template void MappingCartesian::fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &q, - typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, - std::vector &JxW_values, - typename std::vector > &boundary_forms, - typename std::vector > &normal_vectors) const + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + typename std::vector > &boundary_forms, + typename std::vector > &normal_vectors) const { // convert data object to internal // data for this class. fails with @@ -385,7 +385,7 @@ MappingCartesian::fill_fe_subface_values (const typename DoFHandler::c // since the jacobian is diagonal double J = 1.; for (unsigned int d=0;dnormal_directions[face_no]/2)) J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index b88510a027..3f5baa1d12 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -35,10 +35,10 @@ template MappingQ::InternalData::InternalData (const unsigned int n_shape_functions) : MappingQ1::InternalData(n_shape_functions), - use_mapping_q1_on_current_cell(false), - mapping_q1_data(1 << dim) + use_mapping_q1_on_current_cell(false), + mapping_q1_data(1 << dim) { - is_mapping_q1_data=false; + this->is_mapping_q1_data=false; } @@ -525,7 +525,7 @@ MappingQ<3>::set_laplace_on_hex_vector(std::vector > &lohvs) } } else - // not precomputed, then do so now + // not precomputed, then do so now compute_laplace_vector(lohvs); // the sum of weights of the points @@ -686,7 +686,7 @@ MappingQ::compute_mapping_support_points( a.resize(GeometryInfo::vertices_per_cell); for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(vertex_mapping[i]); + a[i] = cell->vertex(this->vertex_mapping[i]); } } @@ -706,27 +706,27 @@ MappingQ::compute_support_points_laplace(const typename Triangulation: switch (dim) { case 2: - // in 2d, add the - // points on the four - // bounding lines to - // the exterior (outer) - // points - add_line_support_points (cell, a); - apply_laplace_vector (laplace_on_quad_vector,a); - break; + // in 2d, add the + // points on the four + // bounding lines to + // the exterior (outer) + // points + add_line_support_points (cell, a); + apply_laplace_vector (laplace_on_quad_vector,a); + break; case 3: - // in 3d also add the - // points located on - // the boundary faces - add_line_support_points (cell, a); - add_quad_support_points (cell, a); - apply_laplace_vector (laplace_on_hex_vector, a); - break; + // in 3d also add the + // points located on + // the boundary faces + add_line_support_points (cell, a); + add_quad_support_points (cell, a); + apply_laplace_vector (laplace_on_hex_vector, a); + break; default: - Assert(false, ExcNotImplemented()); - break; + Assert(false, ExcNotImplemented()); + break; }; } @@ -762,50 +762,50 @@ MappingQ::compute_support_points_simple(const typename Triangulation:: switch (degree) { case 2: - { - a.push_back(middle); - break; - }; + { + a.push_back(middle); + break; + }; case 3: - { - // The four points in the - // cell are located at - // the midpoint between - // the middle point and - // the 4 vertices - for (unsigned int i=0; i::vertices_per_cell; ++i) - a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.); - break; - }; + { + // The four points in the + // cell are located at + // the midpoint between + // the middle point and + // the 4 vertices + for (unsigned int i=0; i::vertices_per_cell; ++i) + a.push_back(middle*2./3.+cell->vertex(this->vertex_mapping[i])/3.); + break; + }; case 4: - { - Assert(a.size()==16, ExcInternalError()); - a.insert(a.end(), 9, Point()); + { + Assert(a.size()==16, ExcInternalError()); + a.insert(a.end(), 9, Point()); - const unsigned int inner_map[8]= - { 0, 1, 2, 5, 8, 7, 6, 3 }; + const unsigned int inner_map[8]= + { 0, 1, 2, 5, 8, 7, 6, 3 }; - // The nine points in the - // cell are located at the - // midpoint between the - // middle point and (the 4 - // vertices and the face - // midpoints) + // The nine points in the + // cell are located at the + // midpoint between the + // middle point and (the 4 + // vertices and the face + // midpoints) - a[16+4]=middle; - for (unsigned int i=0, j=0; i::vertices_per_cell; ++i) - { - a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.; - a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.; - } - break; - }; + a[16+4]=middle; + for (unsigned int i=0, j=0; i::vertices_per_cell; ++i) + { + a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.; + a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.; + } + break; + }; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); }; }; @@ -1092,7 +1092,7 @@ MappingQ::transform_covariant ( else tensor = data->covariant.begin(); } - while (begin!=end) + while (begin!=end) { contract (*(begin++), *(src++), *(tensor++)); } @@ -1125,7 +1125,7 @@ MappingQ::transform_contravariant ( else tensor = data->contravariant.begin(); } - while (begin!=end) + while (begin!=end) { contract (*(begin++), *(tensor++), *(src++)); } @@ -1148,7 +1148,7 @@ Point MappingQ::transform_unit_to_real_cell ( Assert(mdata!=0, ExcInternalError()); mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells - || cell->has_boundary_lines()); + || cell->has_boundary_lines()); typename MappingQ1::InternalData *p_data=0; if (mdata->use_mapping_q1_on_current_cell) diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index b397bed5e4..82c028de59 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -83,7 +83,7 @@ MappingQ1::compute_shapes (const std::vector > &unit_points, template<> const unsigned int MappingQ1::vertex_mapping[2] = { 0, 1 -}; + }; template<> @@ -121,7 +121,7 @@ template<> const unsigned int MappingQ1<2>::vertex_mapping[4] = { 0, 1, 3, 2 -}; + }; template<> @@ -168,7 +168,7 @@ template<> const unsigned int MappingQ1::vertex_mapping[8] = { 0, 1, 4, 5, 3, 2, 7, 6 -}; + }; template<> @@ -276,13 +276,13 @@ MappingQ1::update_each (const UpdateFlags in) const // ignored for the interior of a // cell. if (out & (update_JxW_values - |update_normal_vectors)) + |update_normal_vectors)) out |= update_boundary_forms; if (out & (update_covariant_transformation - | update_JxW_values - | update_boundary_forms - | update_normal_vectors)) + | update_JxW_values + | update_boundary_forms + | update_normal_vectors)) out |= update_contravariant_transformation; return out; @@ -361,14 +361,14 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, // dim so we can subtract 1 // without getting negative // values. - unsigned int nindex = normal_directions[i]/2 + dim; + unsigned int nindex = this->normal_directions[i]/2 + dim; // First tangential has a // non-zero in component // (i+1)%dim, if normal is // non-zero in i. Tensor<1,dim> tangential; - tangential[(nindex+1)%dim] = (normal_directions[i]%2) ? -1 : 1; + tangential[(nindex+1)%dim] = (this->normal_directions[i]%2) ? -1 : 1; data.unit_tangentials[i].resize(n_original_q_points); std::fill (data.unit_tangentials[i].begin(), data.unit_tangentials[i].end(), @@ -561,8 +561,8 @@ MappingQ1::fill_fe_values (const typename DoFHandler::cell_iterator &c Assert (JxW_values.size() == npts, ExcDimensionMismatch(JxW_values.size(), npts)); for (unsigned int point=0; point::compute_fill_face (const typename DoFHandler::cell_iterator switch (dim) { case 2: - { - for (; result != end; ++result, ++tang1) - cross_product (*result, *tang1); - break; - }; + { + for (; result != end; ++result, ++tang1) + cross_product (*result, *tang1); + break; + }; case 3: - { - transform_contravariant(data.aux[1].begin(), - data.aux[1].end(), - data.unit_tangentials[ - face_no+GeometryInfo::faces_per_cell].begin(), - data); - typename std::vector >::const_iterator - tang2 = data.aux[1].begin(); - for (;result != end; ++result, ++tang1, ++tang2) - cross_product (*result, *tang1, *tang2); - break; - }; + { + transform_contravariant(data.aux[1].begin(), + data.aux[1].end(), + data.unit_tangentials[ + face_no+GeometryInfo::faces_per_cell].begin(), + data); + typename std::vector >::const_iterator + tang2 = data.aux[1].begin(); + for (;result != end; ++result, ++tang1, ++tang2) + cross_product (*result, *tang1, *tang2); + break; + }; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } if (update_flags & (update_normal_vectors diff --git a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc index 54e687b92a..0046e34a34 100644 --- a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc @@ -95,12 +95,12 @@ MappingQ1Eulerian::compute_mapping_support_points( // are always numbered first, // we can access them easily for (unsigned int j=0; jvertex_mapping[i]*dim+j); // compute new support point by // old (reference) value and // added shift - a[i] = cell->vertex(vertex_mapping[i]) + shift_vector; + a[i] = cell->vertex(this->vertex_mapping[i]) + shift_vector; } }