From 52bef21ac22c6af07dec9b6f6a0fcda0e53fb982 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 9 Jul 2017 14:04:17 -0400 Subject: [PATCH] Put internal FE functions in specific namespaces. This is the finite element equivalent to 6310b2dee3e and is also necessary for the (possible) unity build patch. --- source/fe/fe_bdm.cc | 79 ++++++----- source/fe/fe_dgp_monomial.cc | 149 ++++++++++---------- source/fe/fe_dgq.cc | 28 ++-- source/fe/fe_enriched.cc | 162 +++++++++++----------- source/fe/fe_face.cc | 30 ++-- source/fe/fe_nedelec.cc | 27 ++-- source/fe/fe_nothing.cc | 17 ++- source/fe/fe_poly_tensor.cc | 181 ++++++++++++------------ source/fe/fe_q.cc | 28 ++-- source/fe/fe_q_base.cc | 121 ++++++++-------- source/fe/fe_q_bubbles.cc | 244 +++++++++++++++++---------------- source/fe/fe_q_hierarchical.cc | 36 +++-- source/fe/fe_system.cc | 25 ++-- 13 files changed, 597 insertions(+), 530 deletions(-) diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index f0e0e22dc8..48255014e6 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -249,48 +249,54 @@ FE_BDM::get_ria_vector (const unsigned int deg) } -namespace +namespace internal { - // This function sets up the values of the polynomials we want to - // take moments with in the quadrature points. In fact, we multiply - // thos by the weights, such that the sum of function values and - // test_values over quadrature points yields the interpolated degree - // of freedom. - template - void - initialize_test_values (std::vector > &test_values, - const Quadrature &quadrature, - const unsigned int deg) + namespace FE_BDM { - PolynomialsP poly(deg); - std::vector > dummy1; - std::vector > dummy2; - std::vector > dummy3; - std::vector > dummy4; + namespace + { + // This function sets up the values of the polynomials we want to + // take moments with in the quadrature points. In fact, we multiply + // thos by the weights, such that the sum of function values and + // test_values over quadrature points yields the interpolated degree + // of freedom. + template + void + initialize_test_values (std::vector > &test_values, + const Quadrature &quadrature, + const unsigned int deg) + { + PolynomialsP poly(deg); + std::vector > dummy1; + std::vector > dummy2; + std::vector > dummy3; + std::vector > dummy4; - test_values.resize(quadrature.size()); + test_values.resize(quadrature.size()); - for (unsigned int k=0; k - void - initialize_test_values (std::vector > &, - const Quadrature<0> &, - const unsigned int) - {} + // This specialization only serves to avoid error messages. Nothing + // useful can be computed in dimension zero and thus the vector + // length stays zero. + template <> + void + initialize_test_values (std::vector > &, + const Quadrature<0> &, + const unsigned int) + {} + } + } } @@ -333,7 +339,7 @@ FE_BDM::initialize_support_points (const unsigned int deg) // point values on faces in 2D. In 3D, this is impossible, since the // moments are only taken with respect to PolynomialsP. if (dim>2) - initialize_test_values(test_values_face, face_points, deg); + internal::FE_BDM::initialize_test_values(test_values_face, face_points, deg); if (deg<=1) return; @@ -347,7 +353,7 @@ FE_BDM::initialize_support_points (const unsigned int deg) // the test functions in the // interior quadrature points - initialize_test_values(test_values_cell, cell_points, deg-2); + internal::FE_BDM::initialize_test_values(test_values_cell, cell_points, deg-2); } @@ -356,4 +362,3 @@ FE_BDM::initialize_support_points (const unsigned int deg) #include "fe_bdm.inst" DEAL_II_NAMESPACE_CLOSE - diff --git a/source/fe/fe_dgp_monomial.cc b/source/fe/fe_dgp_monomial.cc index 27316601e1..9a5f0686aa 100644 --- a/source/fe/fe_dgp_monomial.cc +++ b/source/fe/fe_dgp_monomial.cc @@ -23,88 +23,93 @@ DEAL_II_NAMESPACE_OPEN -// namespace for some functions that are used in this file. -namespace +namespace internal { - // storage of hand-chosen support - // points - // - // For dim=2, dofs_per_cell of - // FE_DGPMonomial(k) is given by - // 0.5(k+1)(k+2), i.e. - // - // k 0 1 2 3 4 5 6 7 - // dofs 1 3 6 10 15 21 28 36 - // - // indirect access of unit points: - // the points for degree k are - // located at - // - // points[start_index[k]..start_index[k+1]-1] - const unsigned int start_index2d[6]= {0,1,4,10,20,35}; - const double points2d[35][2]= + namespace FE_DGPMonomial { - {0,0}, - {0,0},{1,0},{0,1}, - {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5}, - {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5}, - {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5} - }; - - // For dim=3, dofs_per_cell of - // FE_DGPMonomial(k) is given by - // 1./6.(k+1)(k+2)(k+3), i.e. - // - // k 0 1 2 3 4 5 6 7 - // dofs 1 4 10 20 35 56 84 120 - const unsigned int start_index3d[6]= {0,1,5,15/*,35*/}; - const double points3d[35][3]= - { - {0,0,0}, - {0,0,0},{1,0,0},{0,1,0},{0,0,1}, - {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1} - }; + namespace + { + // storage of hand-chosen support + // points + // + // For dim=2, dofs_per_cell of + // FE_DGPMonomial(k) is given by + // 0.5(k+1)(k+2), i.e. + // + // k 0 1 2 3 4 5 6 7 + // dofs 1 3 6 10 15 21 28 36 + // + // indirect access of unit points: + // the points for degree k are + // located at + // + // points[start_index[k]..start_index[k+1]-1] + const unsigned int start_index2d[6]= {0,1,4,10,20,35}; + const double points2d[35][2]= + { + {0,0}, + {0,0},{1,0},{0,1}, + {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5}, + {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5}, + {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5} + }; + + // For dim=3, dofs_per_cell of + // FE_DGPMonomial(k) is given by + // 1./6.(k+1)(k+2)(k+3), i.e. + // + // k 0 1 2 3 4 5 6 7 + // dofs 1 4 10 20 35 56 84 120 + const unsigned int start_index3d[6]= {0,1,5,15/*,35*/}; + const double points3d[35][3]= + { + {0,0,0}, + {0,0,0},{1,0,0},{0,1,0},{0,0,1}, + {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1} + }; - template - void generate_unit_points (const unsigned int, - std::vector > &); + template + void generate_unit_points (const unsigned int, + std::vector > &); - template <> - void generate_unit_points (const unsigned int k, - std::vector > &p) - { - Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1)); - const double h = 1./k; - for (unsigned int i=0; i + void generate_unit_points (const unsigned int k, + std::vector > &p) + { + Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1)); + const double h = 1./k; + for (unsigned int i=0; i - void generate_unit_points (const unsigned int k, - std::vector > &p) - { - Assert(k<=4, ExcNotImplemented()); - Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError()); - for (unsigned int i=0; i + void generate_unit_points (const unsigned int k, + std::vector > &p) { - p[i](0)=points2d[start_index2d[k]+i][0]; - p[i](1)=points2d[start_index2d[k]+i][1]; + Assert(k<=4, ExcNotImplemented()); + Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError()); + for (unsigned int i=0; i - void generate_unit_points (const unsigned int k, - std::vector > &p) - { - Assert(k<=2, ExcNotImplemented()); - Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError()); - for (unsigned int i=0; i + void generate_unit_points (const unsigned int k, + std::vector > &p) { - p[i](0)=points3d[start_index3d[k]+i][0]; - p[i](1)=points3d[start_index3d[k]+i][1]; - p[i](2)=points3d[start_index3d[k]+i][2]; + Assert(k<=2, ExcNotImplemented()); + Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError()); + for (unsigned int i=0; i &source_fe, else { std::vector > unit_points(this->dofs_per_cell); - generate_unit_points(this->degree, unit_points); + internal::FE_DGPMonomial::generate_unit_points(this->degree, unit_points); FullMatrix source_fe_matrix(unit_points.size(), source_fe.dofs_per_cell); for (unsigned int j=0; j > - get_QGaussLobatto_points (const unsigned int degree) + namespace FE_DGQ { - if (degree > 0) - return QGaussLobatto<1>(degree+1).get_points(); - else - return std::vector >(1, Point<1>(0.5)); + namespace + { + std::vector > + get_QGaussLobatto_points (const unsigned int degree) + { + if (degree > 0) + return QGaussLobatto<1>(degree+1).get_points(); + else + return std::vector >(1, Point<1>(0.5)); + } + } } } @@ -48,14 +54,15 @@ template FE_DGQ::FE_DGQ (const unsigned int degree) : FE_Poly, dim, spacedim> - (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))), + (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis + (internal::FE_DGQ::get_QGaussLobatto_points(degree))), FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, true), std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))) { // Compute support points, which are the tensor product of the Lagrange // interpolation points in the constructor. - Quadrature support_quadrature(get_QGaussLobatto_points(degree)); + Quadrature support_quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree)); Assert (support_quadrature.get_points().size() > 0, (typename FiniteElement::ExcFEHasNoSupportPoints ())); this->unit_support_points = support_quadrature.get_points(); @@ -872,7 +879,8 @@ FE_DGQLegendre::clone() const template FE_DGQHermite::FE_DGQHermite (const unsigned int degree) : FE_DGQ(degree < 3 ? - Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree)) + Polynomials::generate_complete_Lagrange_basis + (internal::FE_DGQ::get_QGaussLobatto_points(degree)) : Polynomials::HermiteInterpolation::generate_complete_basis(degree)) {} diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc index b98f2e6e86..4463d3c4b2 100644 --- a/source/fe/fe_enriched.cc +++ b/source/fe/fe_enriched.cc @@ -22,93 +22,99 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - /** - * Auxiliary function to create multiplicity vector from input enrichment functions. - */ - template - std::vector - build_multiplicities(const std::vector > &functions ) + namespace FE_Enriched { - std::vector multiplicities; - multiplicities.push_back(1); // the first one is non-enriched FE - for (unsigned int i = 0; i < functions.size(); i++) - multiplicities.push_back(functions[i].size()); - - return multiplicities; - } - - - /** - * Auxiliary function to build FiniteElement's vector - */ - template - std::vector< const FiniteElement< dim, spacedim > * > - build_fes(const FiniteElement *fe_base, - const std::vector * > &fe_enriched) - { - std::vector< const FiniteElement< dim, spacedim > * > fes; - fes.push_back(fe_base); - for (unsigned int i = 0; i < fe_enriched.size(); i++) - fes.push_back(fe_enriched[i]); - - return fes; - } + namespace + { + /** + * Auxiliary function to create multiplicity vector from input enrichment functions. + */ + template + std::vector + build_multiplicities(const std::vector > &functions ) + { + std::vector multiplicities; + multiplicities.push_back(1); // the first one is non-enriched FE + for (unsigned int i = 0; i < functions.size(); i++) + multiplicities.push_back(functions[i].size()); + return multiplicities; + } - /** - * Auxiliary function which check consistency of the input parameters. - * Returns true if everything is ok. - */ - template - bool - consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes, - const std::vector< unsigned int > &multiplicities, - const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) - { - AssertThrow(fes.size() > 0, - ExcMessage("FEs size should be >=1")); - AssertThrow(fes.size() == multiplicities.size(), - ExcMessage("FEs and multiplicities should have the same size")); - AssertThrow (functions.size() == fes.size() - 1, - ExcDimensionMismatch(functions.size(), fes.size()-1)); + /** + * Auxiliary function to build FiniteElement's vector + */ + template + std::vector< const FiniteElement< dim, spacedim > * > + build_fes(const FiniteElement *fe_base, + const std::vector * > &fe_enriched) + { + std::vector< const FiniteElement< dim, spacedim > * > fes; + fes.push_back(fe_base); + for (unsigned int i = 0; i < fe_enriched.size(); i++) + fes.push_back(fe_enriched[i]); - AssertThrow(multiplicities[0] == 1, - ExcMessage("First multiplicity should be 1")); + return fes; + } - const unsigned int n_comp_base = fes[0]->n_components(); - // start from fe=1 as 0th is always non-enriched FE. - for (unsigned int fe=1; fe < fes.size(); fe++) + /** + * Auxiliary function which check consistency of the input parameters. + * Returns true if everything is ok. + */ + template + bool + consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename dealii::Triangulation::cell_iterator &) > > > &functions) { - const FE_Nothing *fe_nothing = dynamic_cast*>(fes[fe]); - if (fe_nothing) - AssertThrow (fe_nothing->is_dominating(), - ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched")); - - AssertThrow (fes[fe]->n_components() == n_comp_base, - ExcMessage("All elements must have the same number of components")); + AssertThrow(fes.size() > 0, + ExcMessage("FEs size should be >=1")); + AssertThrow(fes.size() == multiplicities.size(), + ExcMessage("FEs and multiplicities should have the same size")); + + AssertThrow (functions.size() == fes.size() - 1, + ExcDimensionMismatch(functions.size(), fes.size()-1)); + + AssertThrow(multiplicities[0] == 1, + ExcMessage("First multiplicity should be 1")); + + const unsigned int n_comp_base = fes[0]->n_components(); + + // start from fe=1 as 0th is always non-enriched FE. + for (unsigned int fe=1; fe < fes.size(); fe++) + { + const FE_Nothing *fe_nothing = dynamic_cast*>(fes[fe]); + if (fe_nothing) + AssertThrow (fe_nothing->is_dominating(), + ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched")); + + AssertThrow (fes[fe]->n_components() == n_comp_base, + ExcMessage("All elements must have the same number of components")); + } + return true; } - return true; - } - /** - * Auxiliary function which determines whether the FiniteElement will be enriched. - */ - template - bool - check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes) - { - // start from fe=1 as 0th is always non-enriched FE. - for (unsigned int fe=1; fe < fes.size(); fe++) - if (dynamic_cast*>(fes[fe]) == nullptr) - // this is not FE_Nothing => there will be enrichment - return true; + /** + * Auxiliary function which determines whether the FiniteElement will be enriched. + */ + template + bool + check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes) + { + // start from fe=1 as 0th is always non-enriched FE. + for (unsigned int fe=1; fe < fes.size(); fe++) + if (dynamic_cast*>(fes[fe]) == nullptr) + // this is not FE_Nothing => there will be enrichment + return true; - return false; + return false; + } + } } } @@ -149,8 +155,8 @@ FE_Enriched::FE_Enriched (const FiniteElement *fe_ba const std::vector * > &fe_enriched, const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) : - FE_Enriched (build_fes(fe_base,fe_enriched), - build_multiplicities(functions), + FE_Enriched (internal::FE_Enriched::build_fes(fe_base,fe_enriched), + internal::FE_Enriched::build_multiplicities(functions), functions) {} @@ -164,11 +170,11 @@ FE_Enriched::FE_Enriched (const std::vector< const FiniteElement< FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities), FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)), enrichments(functions), - is_enriched(check_if_enriched(fes)), + is_enriched(internal::FE_Enriched::check_if_enriched(fes)), fe_system(fes,multiplicities) { // descriptive error are thrown within the function. - Assert(consistency_check(fes,multiplicities,functions), + Assert(internal::FE_Enriched::consistency_check(fes,multiplicities,functions), ExcInternalError()); initialize(fes, multiplicities); diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc index dad2e319f2..448d2b3b6c 100644 --- a/source/fe/fe_face.cc +++ b/source/fe/fe_face.cc @@ -27,25 +27,31 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - std::vector > - get_QGaussLobatto_points (const unsigned int degree) + namespace FE_FaceQ { - if (degree > 0) - return QGaussLobatto<1>(degree+1).get_points(); - else - return std::vector >(1, Point<1>(0.5)); + namespace + { + std::vector > + get_QGaussLobatto_points (const unsigned int degree) + { + if (degree > 0) + return QGaussLobatto<1>(degree+1).get_points(); + else + return std::vector >(1, Point<1>(0.5)); + } + } } } template FE_FaceQ::FE_FaceQ (const unsigned int degree) : - FE_PolyFace, dim, spacedim> ( - TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))), - FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), - std::vector(1,true)) + FE_PolyFace, dim, spacedim> + (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(internal::FE_FaceQ::get_QGaussLobatto_points(degree))), + FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), + std::vector(1,true)) { // initialize unit face support points const unsigned int codim = dim-1; @@ -56,7 +62,7 @@ FE_FaceQ::FE_FaceQ (const unsigned int degree) this->unit_face_support_points[0][d] = 0.5; else { - std::vector > points = get_QGaussLobatto_points(degree); + std::vector > points = internal::FE_FaceQ::get_QGaussLobatto_points(degree); unsigned int k=0; for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz) diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index d9c557968a..d3c5d7eea6 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -44,17 +44,20 @@ DEAL_II_NAMESPACE_OPEN namespace internal { - namespace + namespace FE_Nedelec { - double - get_embedding_computation_tolerance(const unsigned int p) + namespace { - // This heuristic was computed by monitoring the worst residual - // resulting from the least squares computation when computing - // the face embedding matrices in the FE_Nedelec constructor. - // The residual growth is exponential, but is bounded by this - // function up to degree 12. - return 1.e-15*std::exp(std::pow(p,1.075)); + double + get_embedding_computation_tolerance(const unsigned int p) + { + // This heuristic was computed by monitoring the worst residual + // resulting from the least squares computation when computing + // the face embedding matrices in the FE_Nedelec constructor. + // The residual growth is exponential, but is bounded by this + // function up to degree 12. + return 1.e-15*std::exp(std::pow(p,1.075)); + } } } } @@ -107,7 +110,7 @@ FE_Nedelec::FE_Nedelec (const unsigned int order) FETools::compute_face_embedding_matrices (*this, face_embeddings, 0, 0, - internal::get_embedding_computation_tolerance(order)); + internal::FE_Nedelec::get_embedding_computation_tolerance(order)); switch (dim) { @@ -3022,7 +3025,7 @@ FE_Nedelec this_nonconst.reinit_restriction_and_prolongation_matrices (); // Fill prolongation matrices with embedding operators FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true, - internal::get_embedding_computation_tolerance(this->degree)); + internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree)); #ifdef DEBUG_NEDELEC deallog << "Restriction" << std::endl; #endif @@ -3073,7 +3076,7 @@ FE_Nedelec this_nonconst.reinit_restriction_and_prolongation_matrices (); // Fill prolongation matrices with embedding operators FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true, - internal::get_embedding_computation_tolerance(this->degree)); + internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree)); #ifdef DEBUG_NEDELEC deallog << "Restriction" << std::endl; #endif diff --git a/source/fe/fe_nothing.cc b/source/fe/fe_nothing.cc index 71216f1a05..ba059c1a62 100644 --- a/source/fe/fe_nothing.cc +++ b/source/fe/fe_nothing.cc @@ -19,10 +19,16 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - const char * - zero_dof_message = "This element has no shape functions."; + namespace FE_Nothing + { + namespace + { + const char * + zero_dof_message = "This element has no shape functions."; + } + } } @@ -85,8 +91,8 @@ double FE_Nothing::shape_value (const unsigned int /*i*/, const Point & /*p*/) const { - (void)zero_dof_message; - Assert(false,ExcMessage(zero_dof_message)); + (void)internal::FE_Nothing::zero_dof_message; + Assert(false,ExcMessage(internal::FE_Nothing::zero_dof_message)); return 0; } @@ -306,4 +312,3 @@ get_subface_interpolation_matrix (const FiniteElement & /*source_f DEAL_II_NAMESPACE_CLOSE - diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index a88e6155d7..d3d9626b44 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -23,109 +23,116 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { -//--------------------------------------------------------------------------- -// Utility method, which is used to determine the change of sign for -// the DoFs on the faces of the given cell. -//--------------------------------------------------------------------------- - - /** - * On noncartesian grids, the sign of the DoFs associated with the faces of - * the elements has to be changed in some cases. This procedure implements an - * algorithm, which determines the DoFs, which need this sign change for a - * given cell. - */ - void - get_face_sign_change_rt (const Triangulation<1>::cell_iterator &, - const unsigned int, - std::vector &face_sign) + namespace FE_PolyTensor { - // nothing to do in 1d - std::fill (face_sign.begin (), face_sign.end (), 1.0); - } - - + namespace + { + //--------------------------------------------------------------------------- + // Utility method, which is used to determine the change of sign for + // the DoFs on the faces of the given cell. + //--------------------------------------------------------------------------- + + /** + * On noncartesian grids, the sign of the DoFs associated with the faces of + * the elements has to be changed in some cases. This procedure implements an + * algorithm, which determines the DoFs, which need this sign change for a + * given cell. + */ + void + get_face_sign_change_rt (const dealii::Triangulation<1>::cell_iterator &, + const unsigned int, + std::vector &face_sign) + { + // nothing to do in 1d + std::fill (face_sign.begin (), face_sign.end (), 1.0); + } - void - get_face_sign_change_rt (const Triangulation<2>::cell_iterator &cell, - const unsigned int dofs_per_face, - std::vector &face_sign) - { - const unsigned int dim = 2; - const unsigned int spacedim = 2; - // Default is no sign - // change. I.e. multiply by one. - std::fill (face_sign.begin (), face_sign.end (), 1.0); - for (unsigned int f = GeometryInfo::faces_per_cell / 2; - f < GeometryInfo::faces_per_cell; ++f) + void + get_face_sign_change_rt (const dealii::Triangulation<2>::cell_iterator &cell, + const unsigned int dofs_per_face, + std::vector &face_sign) { - Triangulation::face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - const unsigned int nn = cell->neighbor_face_no(f); + const unsigned int dim = 2; + const unsigned int spacedim = 2; - if (nn < GeometryInfo::faces_per_cell / 2) - for (unsigned int j = 0; j < dofs_per_face; ++j) - { - Assert (f * dofs_per_face + j < face_sign.size(), - ExcInternalError()); + // Default is no sign + // change. I.e. multiply by one. + std::fill (face_sign.begin (), face_sign.end (), 1.0); -//TODO: This is probably only going to work for those elements for which all dofs are face dofs - face_sign[f * dofs_per_face + j] = -1.0; - } + for (unsigned int f = GeometryInfo::faces_per_cell / 2; + f < GeometryInfo::faces_per_cell; ++f) + { + dealii::Triangulation::face_iterator face = cell->face (f); + if (!face->at_boundary ()) + { + const unsigned int nn = cell->neighbor_face_no(f); + + if (nn < GeometryInfo::faces_per_cell / 2) + for (unsigned int j = 0; j < dofs_per_face; ++j) + { + Assert (f * dofs_per_face + j < face_sign.size(), + ExcInternalError()); + + //TODO: This is probably only going to work for those elements for which all dofs are face dofs + face_sign[f * dofs_per_face + j] = -1.0; + } + } } } - } - void - get_face_sign_change_rt (const Triangulation<3>::cell_iterator &/*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) - { - std::fill (face_sign.begin (), face_sign.end (), 1.0); -//TODO: think about what it would take here - } + void + get_face_sign_change_rt (const dealii::Triangulation<3>::cell_iterator &/*cell*/, + const unsigned int /*dofs_per_face*/, + std::vector &face_sign) + { + std::fill (face_sign.begin (), face_sign.end (), 1.0); + //TODO: think about what it would take here + } - void - get_face_sign_change_nedelec (const Triangulation<1>::cell_iterator &/*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) - { - // nothing to do in 1d - std::fill (face_sign.begin (), face_sign.end (), 1.0); - } + void + get_face_sign_change_nedelec (const dealii::Triangulation<1>::cell_iterator &/*cell*/, + const unsigned int /*dofs_per_face*/, + std::vector &face_sign) + { + // nothing to do in 1d + std::fill (face_sign.begin (), face_sign.end (), 1.0); + } - void - get_face_sign_change_nedelec (const Triangulation<2>::cell_iterator &/*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) - { - std::fill (face_sign.begin (), face_sign.end (), 1.0); -//TODO: think about what it would take here - } + void + get_face_sign_change_nedelec (const dealii::Triangulation<2>::cell_iterator &/*cell*/, + const unsigned int /*dofs_per_face*/, + std::vector &face_sign) + { + std::fill (face_sign.begin (), face_sign.end (), 1.0); + //TODO: think about what it would take here + } - void - get_face_sign_change_nedelec (const Triangulation<3>::cell_iterator &cell, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) - { - const unsigned int dim = 3; - std::fill (face_sign.begin (), face_sign.end (), 1.0); -//TODO: This is probably only going to work for those elements for which all dofs are face dofs - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - if (!(cell->line_orientation (l))) - face_sign[l] = -1.0; + void + get_face_sign_change_nedelec (const dealii::Triangulation<3>::cell_iterator &cell, + const unsigned int /*dofs_per_face*/, + std::vector &face_sign) + { + const unsigned int dim = 3; + std::fill (face_sign.begin (), face_sign.end (), 1.0); + //TODO: This is probably only going to work for those elements for which all dofs are face dofs + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + if (!(cell->line_orientation (l))) + face_sign[l] = -1.0; + } + } } } @@ -335,9 +342,9 @@ fill_fe_values std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) @@ -816,10 +823,10 @@ fill_fe_face_values std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { @@ -1295,10 +1302,10 @@ fill_fe_subface_values std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index 339205589a..a6941c5f76 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -25,19 +25,25 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - std::vector > - get_QGaussLobatto_points (const unsigned int degree) + namespace FE_Q { - if (degree > 0) - return QGaussLobatto<1>(degree+1).get_points(); - else + namespace + { + std::vector > + get_QGaussLobatto_points (const unsigned int degree) { - typedef FE_Q_Base, 1, 1> FEQ; - AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0()); + if (degree > 0) + return QGaussLobatto<1>(degree+1).get_points(); + else + { + typedef FE_Q_Base, 1, 1> FEQ; + AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0()); + } + return std::vector >(); } - return std::vector >(); + } } } @@ -47,13 +53,13 @@ template FE_Q::FE_Q (const unsigned int degree) : FE_Q_Base, dim, spacedim> - (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))), + (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(internal::FE_Q::get_QGaussLobatto_points(degree))), FiniteElementData(this->get_dpo_vector(degree), 1, degree, FiniteElementData::H1), std::vector (1, false)) { - this->initialize(get_QGaussLobatto_points(degree)); + this->initialize(internal::FE_Q::get_QGaussLobatto_points(degree)); } diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index b958a14177..ad3cc46af1 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -35,66 +35,69 @@ DEAL_II_NAMESPACE_OPEN -namespace FE_Q_Helper +namespace internal { - namespace + namespace FE_Q_Base { - // get the renumbering for faces - template - inline - std::vector - face_lexicographic_to_hierarchic_numbering (const unsigned int degree) + namespace { - std::vector dpo(dim, 1U); - for (unsigned int i=1; i face_data(dpo,1,degree); - std::vector face_renumber (face_data.dofs_per_cell); - FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber); - return face_renumber; - } + // get the renumbering for faces + template + inline + std::vector + face_lexicographic_to_hierarchic_numbering (const unsigned int degree) + { + std::vector dpo(dim, 1U); + for (unsigned int i=1; i face_data(dpo,1,degree); + std::vector face_renumber (face_data.dofs_per_cell); + FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber); + return face_renumber; + } - // dummy specialization for dim == 1 to avoid linker errors - template <> - inline - std::vector - face_lexicographic_to_hierarchic_numbering<1> (const unsigned int) - { - return std::vector(); - } + // dummy specialization for dim == 1 to avoid linker errors + template <> + inline + std::vector + face_lexicographic_to_hierarchic_numbering<1> (const unsigned int) + { + return std::vector(); + } - // in get_restriction_matrix() and get_prolongation_matrix(), want to undo - // tensorization on inner loops for performance reasons. this clears a - // dim-array - template - inline - void - zero_indices (unsigned int (&indices)[dim]) - { - for (unsigned int d=0; d + inline + void + zero_indices (unsigned int (&indices)[dim]) + { + for (unsigned int d=0; d - inline - void - increment_indices (unsigned int (&indices)[dim], - const unsigned int dofs1d) - { - ++indices[0]; - for (int d=0; d + inline + void + increment_indices (unsigned int (&indices)[dim], + const unsigned int dofs1d) + { + ++indices[0]; + for (int d=0; d::Implementation const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); + internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); @@ -329,7 +332,7 @@ struct FE_Q_Base::Implementation const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); + internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); @@ -928,7 +931,7 @@ void FE_Q_Base::initialize_unit_face_support_points // find renumbering of faces and assign from values of quadrature std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_degree); + internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering(q_degree); Quadrature<1> support_1d(points); Quadrature support_quadrature(support_1d); this->unit_face_support_points.resize(support_quadrature.size()); @@ -1284,11 +1287,11 @@ FE_Q_Base // now expand from 1D info. block innermost dimension (x_0) in order to // avoid difficult checks at innermost loop unsigned int j_indices[dim]; - FE_Q_Helper::zero_indices (j_indices); + internal::FE_Q_Base::zero_indices (j_indices); for (unsigned int j=0; j (i_indices); + internal::FE_Q_Base::zero_indices (i_indices); for (unsigned int i=0; i // update indices that denote the tensor product position. a bit // fuzzy and therefore not done for innermost x_0 direction - FE_Q_Helper::increment_indices (i_indices, dofs1d); + internal::FE_Q_Base::increment_indices (i_indices, dofs1d); } Assert (i_indices[dim-1] == 1, ExcInternalError()); - FE_Q_Helper::increment_indices (j_indices, dofs1d); + internal::FE_Q_Base::increment_indices (j_indices, dofs1d); } // the discontinuous node is simply mapped on the discontinuous node on @@ -1428,7 +1431,7 @@ FE_Q_Base this->poly_space.compute_value(index_map_inverse[j], point); } unsigned int j_indices[dim]; - FE_Q_Helper::zero_indices (j_indices); + internal::FE_Q_Base::zero_indices (j_indices); double sum_check = 0; for (unsigned int j = 0; j my_restriction(mother_dof,child_dof)=val; sum_check += val; } - FE_Q_Helper::increment_indices (j_indices, dofs1d); + internal::FE_Q_Base::increment_indices (j_indices, dofs1d); } Assert (std::fabs(sum_check-1) < std::max(eps, 5e-16*std::sqrt(this->dofs_per_cell)), diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index c82b9cd83d..4312fb0542 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -38,131 +38,135 @@ DEAL_II_NAMESPACE_OPEN -namespace FE_Q_Bubbles_Helper +namespace internal { - namespace + namespace FE_Q_Bubbles { - template - inline - void - compute_embedding_matrices(const FE_Q_Bubbles &fe, - std::vector > > &matrices, - const bool isotropic_only) + namespace { - const unsigned int dpc = fe.dofs_per_cell; - const unsigned int degree = fe.degree; + template + inline + void + compute_embedding_matrices(const dealii::FE_Q_Bubbles &fe, + std::vector > > &matrices, + const bool isotropic_only) + { + const unsigned int dpc = fe.dofs_per_cell; + const unsigned int degree = fe.degree; - // Initialize quadrature formula on fine cells - std::unique_ptr > q_fine; - Quadrature<1> q_dummy(std::vector >(1), std::vector (1,1.)); - switch (dim) - { - case 1: - if (spacedim==1) - q_fine.reset(new QGauss (degree+1)); - else if (spacedim==2) - q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), q_dummy)); - else - q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), q_dummy, q_dummy)); - break; - case 2: - if (spacedim==2) + // Initialize quadrature formula on fine cells + std::unique_ptr > q_fine; + Quadrature<1> q_dummy(std::vector >(1), std::vector (1,1.)); + switch (dim) + { + case 1: + if (spacedim==1) + q_fine.reset(new QGauss (degree+1)); + else if (spacedim==2) + q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), q_dummy)); + else + q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), q_dummy, q_dummy)); + break; + case 2: + if (spacedim==2) + q_fine.reset(new QGauss (degree+1)); + else + q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy)); + break; + case 3: q_fine.reset(new QGauss (degree+1)); - else - q_fine.reset(new QAnisotropic(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy)); - break; - case 3: - q_fine.reset(new QGauss (degree+1)); - break; - default: - Assert(false, ExcInternalError()); - } + break; + default: + Assert(false, ExcInternalError()); + } - Assert(q_fine.get() != nullptr, ExcInternalError()); - const unsigned int nq = q_fine->size(); + Assert(q_fine.get() != nullptr, ExcInternalError()); + const unsigned int nq = q_fine->size(); - // loop over all possible refinement cases - unsigned int ref_case = (isotropic_only) - ? RefinementCase::isotropic_refinement - : RefinementCase::cut_x; - for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - { - const unsigned int nc - = GeometryInfo::n_children(RefinementCase(ref_case)); - - for (unsigned int i=0; i tr; - GridGenerator::hyper_cube (tr, 0, 1); - tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); - tr.execute_coarsening_and_refinement(); - - DoFHandler dh(tr); - dh.distribute_dofs(fe); - - FEValues fine (StaticMappingQ1::mapping, fe, *q_fine, - update_quadrature_points - | update_JxW_values | update_values); - - const unsigned int n_dofs = dh.n_dofs(); - - FullMatrix fine_mass(n_dofs); - FullMatrix coarse_rhs_matrix(n_dofs, dpc); - - std::vector > child_ldi - (nc, std::vector(fe.dofs_per_cell)); - - //now create the mass matrix and all the right_hand sides - unsigned int child_no = 0; - typename dealii::DoFHandler::active_cell_iterator cell - = dh.begin_active(); - for (; cell!=dh.end(); ++cell, ++child_no) - { - fine.reinit(cell); - cell->get_dof_indices(child_ldi[child_no]); - - for (unsigned int q=0; q quad_tmp; - for (unsigned int k=0; k solution (n_dofs, dpc); - fine_mass.gauss_jordan(); - fine_mass.mmult(solution, coarse_rhs_matrix); - - //and distribute to the fine cell matrices - for (unsigned int child_no=0; child_no 1.e-12) - matrices[ref_case-1][child_no](i,j)=solution(gdi, j); - } - } + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + { + const unsigned int nc + = GeometryInfo::n_children(RefinementCase(ref_case)); + + for (unsigned int i=0; i tr; + GridGenerator::hyper_cube (tr, 0, 1); + tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); + tr.execute_coarsening_and_refinement(); + + dealii::DoFHandler dh(tr); + dh.distribute_dofs(fe); + + dealii::FEValues fine + (StaticMappingQ1::mapping, fe, *q_fine, + update_quadrature_points + | update_JxW_values | update_values); + + const unsigned int n_dofs = dh.n_dofs(); + + FullMatrix fine_mass(n_dofs); + FullMatrix coarse_rhs_matrix(n_dofs, dpc); + + std::vector > child_ldi + (nc, std::vector(fe.dofs_per_cell)); + + //now create the mass matrix and all the right_hand sides + unsigned int child_no = 0; + typename dealii::DoFHandler::active_cell_iterator cell + = dh.begin_active(); + for (; cell!=dh.end(); ++cell, ++child_no) + { + fine.reinit(cell); + cell->get_dof_indices(child_ldi[child_no]); + + for (unsigned int q=0; q quad_tmp; + for (unsigned int k=0; k solution (n_dofs, dpc); + fine_mass.gauss_jordan(); + fine_mass.mmult(solution, coarse_rhs_matrix); + + //and distribute to the fine cell matrices + for (unsigned int child_no=0; child_no 1.e-12) + matrices[ref_case-1][child_no](i,j)=solution(gdi, j); + } + } + } } } } @@ -196,7 +200,7 @@ FE_Q_Bubbles::FE_Q_Bubbles (const unsigned int q_degree) this->reinit_restriction_and_prolongation_matrices(); if (dim == spacedim) { - FE_Q_Bubbles_Helper::compute_embedding_matrices + internal::FE_Q_Bubbles::compute_embedding_matrices (*this, this->prolongation, false); // Fill restriction matrices with L2-projection FETools::compute_projection_matrices (*this, this->restriction); @@ -234,7 +238,7 @@ FE_Q_Bubbles::FE_Q_Bubbles (const Quadrature<1> &points) this->reinit_restriction_and_prolongation_matrices(); if (dim == spacedim) { - FE_Q_Bubbles_Helper::compute_embedding_matrices + internal::FE_Q_Bubbles::compute_embedding_matrices (*this, this->prolongation, false); // Fill restriction matrices with L2-projection FETools::compute_projection_matrices (*this, this->restriction); diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 2e55a25841..184e0a5c71 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -28,23 +28,29 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - /** - * A function which maps in[i] to i,i.e. output[in[i]] = i; - */ - inline - std::vector - invert_numbering (const std::vector &in) + namespace FE_Q_Hierarchical { - std::vector out (in.size()); - for (unsigned int i=0; i + invert_numbering (const std::vector &in) { - Assert (in[i] < out.size(), - dealii::ExcIndexRange(in[i],0,out.size())); - out[in[i]]=i; + std::vector out (in.size()); + for (unsigned int i=0; i:: face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree) { FiniteElementData fe_data(FE_Q_Hierarchical::get_dpo_vector(degree),1,degree); - return invert_numbering(FE_Q_Hierarchical:: - hierarchic_to_fe_q_hierarchical_numbering (fe_data)); + return internal::FE_Q_Hierarchical::invert_numbering + (FE_Q_Hierarchical::hierarchic_to_fe_q_hierarchical_numbering (fe_data)); } diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 0f93679872..40ee23cd78 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -30,18 +30,21 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - bool IsNonZero (unsigned int i) + namespace FESystem { - return i>0; - } - - unsigned int count_nonzeros(const std::vector &vec) - { - return std::count_if(vec.begin(), vec.end(), IsNonZero); + namespace + { + unsigned int count_nonzeros(const std::vector &vec) + { + return std::count_if(vec.begin(), vec.end(), [](const unsigned int i) + { + return i > 0; + }); + } + } } - } /* ----------------------- FESystem::InternalData ------------------- */ @@ -275,7 +278,7 @@ FESystem::FESystem ( FiniteElement (FETools::Compositing::multiply_dof_numbers(fes, multiplicities), FETools::Compositing::compute_restriction_is_additive_flags (fes, multiplicities), FETools::Compositing::compute_nonzero_components(fes, multiplicities)), - base_elements(count_nonzeros(multiplicities)) + base_elements(internal::FESystem::count_nonzeros(multiplicities)) { initialize(fes, multiplicities); } @@ -1458,7 +1461,7 @@ void FESystem::initialize (const std::vector 0, ExcMessage ("Need to pass at least one finite element.")); - Assert (count_nonzeros(multiplicities) > 0, + Assert (internal::FESystem::count_nonzeros(multiplicities) > 0, ExcMessage("You only passed FiniteElements with multiplicity 0.")); // Note that we need to skip every fe with multiplicity 0 in the following block of code -- 2.39.5