From 14e8ab23f571fbf1789834d7283653e42a98e366 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Fri, 28 Nov 2003 11:52:35 +0000 Subject: [PATCH] new function tensor_degree in FiniteElementData git-svn-id: https://svn.dealii.org/trunk@8207 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_base.h | 40 ++++++++++++++++++- deal.II/deal.II/source/fe/fe_data.cc | 6 ++- deal.II/deal.II/source/fe/fe_dgp.cc | 6 +-- deal.II/deal.II/source/fe/fe_dgq.cc | 16 ++++---- deal.II/deal.II/source/fe/fe_q.cc | 6 +-- .../deal.II/source/fe/fe_q_hierarchical.cc | 6 +-- .../deal.II/source/fe/fe_raviart_thomas.cc | 4 +- deal.II/deal.II/source/fe/fe_system.cc | 21 +++++++--- deal.II/doc/news/c-4-0.html | 11 ++++- tests/fe/fe_data_test.cc | 4 +- 10 files changed, 89 insertions(+), 31 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index b3561e90b7..e3e6c43a76 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -45,7 +45,7 @@ template class FESystem; * @ref{FiniteElementBase} should be private base classes of * @ref{FiniteElement}. * - * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2001 + * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2001, 2003 */ template class FiniteElementData @@ -137,6 +137,13 @@ class FiniteElementData */ const unsigned int components; + /** + * Maximal polynomial degree of a + * shape function in a single + * coordinate direction. + */ + const unsigned int degree; + /** * Default * constructor. Constructs an @@ -162,9 +169,17 @@ class FiniteElementData * * Hence this constructor requires * @p{dofs_per_object.size()==dim+1}. + * + * @param n_components Number of + * vector components of the + * element. + * @param degree + * Maximal polynomial degree in a + * single direction. */ FiniteElementData (const std::vector &dofs_per_object, - const unsigned int n_components); + const unsigned int n_components, + const unsigned int degree = deal_II_numbers::invalid_unsigned_int); /** * Return the @p{dofs_per_vertex}. @@ -201,6 +216,17 @@ class FiniteElementData */ unsigned int n_components () const; + /** + * Maximal polynomial degree of a + * shape function in a single + * coordinate direction. + * + * This function can be used to + * determine the optimal + * quadrature rule. + */ + unsigned int tensor_degree () const; + /** * Comparison operator. */ @@ -1778,6 +1804,16 @@ FiniteElementData::n_components () const +template +inline +unsigned int +FiniteElementData::tensor_degree () const +{ + return degree; +} + + + template inline std::pair diff --git a/deal.II/deal.II/source/fe/fe_data.cc b/deal.II/deal.II/source/fe/fe_data.cc index bd13f0ee25..6ecadf0b5a 100644 --- a/deal.II/deal.II/source/fe/fe_data.cc +++ b/deal.II/deal.II/source/fe/fe_data.cc @@ -18,7 +18,8 @@ template FiniteElementData::FiniteElementData (const std::vector &dofs_per_object, - const unsigned int n_components) : + const unsigned int n_components, + const unsigned int degree) : dofs_per_vertex(dofs_per_object[0]), dofs_per_line(dofs_per_object[1]), dofs_per_quad(dim>1? dofs_per_object[2]:0), @@ -45,7 +46,8 @@ FiniteElementData::FiniteElementData (const std::vector &dofs GeometryInfo::lines_per_cell * dofs_per_line + GeometryInfo::quads_per_cell * dofs_per_quad + GeometryInfo::hexes_per_cell * dofs_per_hex), - components(n_components) + components(n_components), + degree(degree) { Assert(dofs_per_object.size()==dim+1, ExcDimensionMismatch(dofs_per_object.size()-1,dim)); } diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 90d06302b1..9230f796cf 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -30,9 +30,9 @@ template FE_DGP::FE_DGP (const unsigned int degree) : - FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell,true), - std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + FiniteElement (FiniteElementData(get_dpo_vector(degree), 1, degree), + 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))), degree(degree), polynomial_space (Polynomials::Legendre::generate_complete_basis(degree)) diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 7e1d50b4c4..f128640556 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -45,7 +45,7 @@ namespace // given an integer N, compute its // integer square root (if it // exists, otherwise give up) - unsigned int int_sqrt (const unsigned int N) + inline unsigned int int_sqrt (const unsigned int N) { for (unsigned int i=0; i<=N; ++i) if (i*i == N) @@ -58,7 +58,7 @@ namespace // given an integer N, compute its // integer cube root (if it // exists, otherwise give up) - unsigned int int_cuberoot (const unsigned int N) + inline unsigned int int_cuberoot (const unsigned int N) { for (unsigned int i=0; i<=N; ++i) if (i*i*i == N) @@ -71,7 +71,7 @@ namespace // given N, generate i=0...N-1 // equidistant points in the // interior of the interval [0,1] - Point<1> + inline Point<1> generate_unit_point (const unsigned int i, const unsigned int N, const int2type<1> ) @@ -90,7 +90,7 @@ namespace // given N, generate i=0...N-1 // equidistant points in the domain // [0,1]^2 - Point<2> + inline Point<2> generate_unit_point (const unsigned int i, const unsigned int N, const int2type<2> ) @@ -116,7 +116,7 @@ namespace // given N, generate i=0...N-1 // equidistant points in the domain // [0,1]^3 - Point<3> + inline Point<3> generate_unit_point (const unsigned int i, const unsigned int N, const int2type<3> ) @@ -144,10 +144,10 @@ namespace template FE_DGQ::FE_DGQ (const unsigned int degree) : - FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + FiniteElement (FiniteElementData(get_dpo_vector(degree), 1, degree), + std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, true), - std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector >(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))), degree(degree), polynomial_space (Polynomials::LagrangeEquidistant::generate_complete_basis(degree)) diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index db4bfa1b1c..53e7be3370 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -431,10 +431,10 @@ namespace FE_Q_Helper template FE_Q::FE_Q (const unsigned int degree) : - FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector (FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + FiniteElement (FiniteElementData(get_dpo_vector(degree),1, degree), + std::vector (FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, false), - std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector >(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))), degree(degree), renumber(lexicographic_to_hierarchic_numbering (*this, degree)), diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 7a3413c53a..edf384f409 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -47,10 +47,10 @@ namespace template FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) : - FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector (FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + FiniteElement (FiniteElementData(get_dpo_vector(degree),1, degree), + std::vector (FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, false), - std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector >(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))), degree(degree), renumber(lexicographic_to_hierarchic_numbering (*this, degree)), diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index b63f5b854b..ddcc7c65d9 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -128,9 +128,9 @@ template FE_RaviartThomas::FE_RaviartThomas (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree), - dim), + dim, degree+1), get_ria_vector (degree), - std::vector >(FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell, + std::vector >(FiniteElementData(get_dpo_vector(degree),dim,degree+1).dofs_per_cell, std::vector(dim,true))), degree(degree), polynomials (create_polynomials(degree)), diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 6aac1b39ff..9ae5a5cbb0 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -20,7 +20,7 @@ #include #include #include - +#include #ifdef HAVE_STD_STRINGSTREAM # include #else @@ -1732,7 +1732,7 @@ FESystem::multiply_dof_numbers (const FiniteElementData &fe_data, if (dim>1) dpo.push_back(fe_data.dofs_per_quad * N); if (dim>2) dpo.push_back(fe_data.dofs_per_hex * N); - return FiniteElementData (dpo, fe_data.n_components() * N); + return FiniteElementData (dpo, fe_data.n_components() * N, fe_data.tensor_degree()); } @@ -1843,10 +1843,17 @@ FESystem::multiply_dof_numbers (const FiniteElementData &fe1, dpo.push_back(fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2); if (dim>1) dpo.push_back(fe1.dofs_per_quad * N1 + fe2.dofs_per_quad * N2); if (dim>2) dpo.push_back(fe1.dofs_per_hex * N1 + fe2.dofs_per_hex * N2); - + + // degree is the maximal degree of + // the components. max also makes + // sure that one unknown degree + // makes the degree of the system + // unknown. + unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); return FiniteElementData (dpo, fe1.n_components() * N1 + - fe2.n_components() * N2); + fe2.n_components() * N2, + degree); } @@ -1873,11 +1880,13 @@ FESystem::multiply_dof_numbers (const FiniteElementData &fe1, if (dim>2) dpo.push_back(fe1.dofs_per_hex * N1 + fe2.dofs_per_hex * N2 + fe3.dofs_per_hex * N3); - + // degree is the maximal degree of the components + unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); + degree = std::max(degree, fe3.tensor_degree()); return FiniteElementData (dpo, fe1.n_components() * N1 + fe2.n_components() * N2 + - fe3.n_components() * N3); + fe3.n_components() * N3, degree); } diff --git a/deal.II/doc/news/c-4-0.html b/deal.II/doc/news/c-4-0.html index 7de7e42175..22a5c5f8b9 100644 --- a/deal.II/doc/news/c-4-0.html +++ b/deal.II/doc/news/c-4-0.html @@ -303,10 +303,19 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

deal.II

    + +
  1. + Improved: FiniteElementData has a function + tensor_degree(), returning the degree of the + polynomial space suitable for choosing a tensor product quadrature formula. +
    + (GK 2003/11/28) +

    +
  2. New: Long requested but never implemented before in the library: there is now a function GridTool::find_active_cell_around_point + GridTool::find_active_cell_around_point that, given a point, finds the active cell in which this point lies.
    diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc index 04240d56ea..e542e044a0 100644 --- a/tests/fe/fe_data_test.cc +++ b/tests/fe/fe_data_test.cc @@ -95,7 +95,8 @@ void test_fe_datas() deallog << "first_face_quad_index=" << fe_data->first_face_quad_index << std::endl; deallog << "dofs_per_face=" << fe_data->dofs_per_face << std::endl; deallog << "dofs_per_cell=" << fe_data->dofs_per_cell << std::endl; - deallog << "components=" << fe_data->components << std::endl; + deallog << "components=" << fe_data->components << std::endl + << "degree=" << fe_data->tensor_degree() << std::endl; for (unsigned int f=0;f::faces_per_cell;++f) { @@ -105,6 +106,7 @@ void test_fe_datas() deallog << '\t' << s; deallog << std::endl; } + deallog << std::endl; } // delete all FiniteElementDatas -- 2.39.5