From: heister Date: Fri, 28 Sep 2012 19:32:18 +0000 (+0000) Subject: FESystem can now deal with n_elements==0 for a block. Remove a ton of dublicated... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=278de496c399c983b8bfcb0b277f1507fa94d960;p=dealii-svn.git FESystem can now deal with n_elements==0 for a block. Remove a ton of dublicated code in the process. git-svn-id: https://svn.dealii.org/trunk@26862 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 318bd56f1c..d98241be0e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -87,6 +87,10 @@ DoFHandler, in particular removal of specializations.

Specific improvements

    +
  1. Fixed: FESystem can now deal with n_elements==0 for a block. +
    +(Timo Heister, 2012/09/28) +
  2. Fixed: ParameterHandler::print_parameters, when using ParameterHandler::OutputStyle::LaTeX would always print a list of parameters in each section as a latex itemized environment. diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index 9a0b56ccce..d8df1cf10c 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -850,59 +850,16 @@ class FESystem : public FiniteElement * cell. */ static FiniteElementData - multiply_dof_numbers (const FiniteElementData &fe_data, - const unsigned int N); - - /** - * Same as above for mixed elements - * with two different sub-elements. - */ - static FiniteElementData - multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2); - - /** - * Same as above for mixed elements - * with three different sub-elements. - */ - static FiniteElementData - multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3); - - /** - * with 4 different sub-elements - */ - static FiniteElementData - multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3, - const FiniteElementData &fe4, - const unsigned int N4); - - - /** - * and with 5 different sub-elements - */ - static FiniteElementData - multiply_dof_numbers (const FiniteElementData &fe1, + multiply_dof_numbers (const FiniteElement *fe1, const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3, - const FiniteElementData &fe4, - const unsigned int N4, - const FiniteElementData &fe5, - const unsigned int N5); + const FiniteElement *fe2=NULL, + const unsigned int N2=0, + const FiniteElement *fe3=NULL, + const unsigned int N3=0, + const FiniteElement *fe4=NULL, + const unsigned int N4=0, + const FiniteElement *fe5=NULL, + const unsigned int N5=0); /** * Same as above but for @@ -1083,7 +1040,8 @@ class FESystem : public FiniteElement * @p prolongation * matrices. */ - void initialize(); + void initialize (const std::vector*> &fes, + const std::vector &multiplicities); /** * Used by @p initialize. diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 1c024921bd..707b22f672 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -27,6 +27,16 @@ DEAL_II_NAMESPACE_OPEN +namespace +{ + bool IsNonZero (unsigned int i) { return i>0; } + + unsigned int count_nonzeros(const std::vector & vec) + { + return std::count_if(vec.begin(), vec.end(), IsNonZero); + } + +}; /* ----------------------- FESystem::InternalData ------------------- */ @@ -149,16 +159,16 @@ const unsigned int FESystem::invalid_face_number; template FESystem::FESystem (const FiniteElement &fe, const unsigned int n_elements) : - FiniteElement (multiply_dof_numbers(fe, n_elements), + 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((n_elements>0)) { - this->base_to_block_indices.reinit(0, 0); - this->base_to_block_indices.push_back(n_elements); - base_elements[0] = ElementPair(fe.clone(), n_elements); - base_elements[0].first->subscribe (typeid(*this).name()); - initialize (); + std::vector*> fes; + fes.push_back(&fe); + std::vector multiplicities; + multiplicities.push_back(n_elements); + initialize(fes, multiplicities); } @@ -168,21 +178,20 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, const FiniteElement &fe2, const unsigned int n2) : - FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2), + FiniteElement (multiply_dof_numbers(&fe1, n1, &fe2, n2), compute_restriction_is_additive_flags (fe1, n1, fe2, n2), compute_nonzero_components(fe1, n1, fe2, n2)), - base_elements(2) + base_elements((n1>0)+(n2>0)) { - this->base_to_block_indices.reinit(0, 0); - this->base_to_block_indices.push_back(n1); - this->base_to_block_indices.push_back(n2); - base_elements[0] = ElementPair(fe1.clone(), n1); - base_elements[0].first->subscribe (typeid(*this).name()); - base_elements[1] = ElementPair(fe2.clone(), n2); - base_elements[1].first->subscribe (typeid(*this).name()); - initialize (); + std::vector*> fes; + fes.push_back(&fe1); + fes.push_back(&fe2); + std::vector multiplicities; + multiplicities.push_back(n1); + multiplicities.push_back(n2); + initialize(fes, multiplicities); } @@ -194,28 +203,26 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n2, const FiniteElement &fe3, const unsigned int n3) : - FiniteElement (multiply_dof_numbers(fe1, n1, - fe2, n2, - fe3, n3), + FiniteElement (multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3), compute_restriction_is_additive_flags (fe1, n1, fe2, n2, fe3, n3), compute_nonzero_components(fe1, n1, fe2, n2, fe3, n3)), - base_elements(3) + base_elements((n1>0)+(n2>0)+(n3>0)) { - this->base_to_block_indices.reinit(0, 0); - this->base_to_block_indices.push_back(n1); - this->base_to_block_indices.push_back(n2); - this->base_to_block_indices.push_back(n3); - base_elements[0] = ElementPair(fe1.clone(), n1); - base_elements[0].first->subscribe (typeid(*this).name()); - base_elements[1] = ElementPair(fe2.clone(), n2); - base_elements[1].first->subscribe (typeid(*this).name()); - base_elements[2] = ElementPair(fe3.clone(), n3); - base_elements[2].first->subscribe (typeid(*this).name()); - initialize (); + std::vector*> fes; + fes.push_back(&fe1); + fes.push_back(&fe2); + fes.push_back(&fe3); + std::vector multiplicities; + multiplicities.push_back(n1); + multiplicities.push_back(n2); + multiplicities.push_back(n3); + initialize(fes, multiplicities); } @@ -228,10 +235,10 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n3, const FiniteElement &fe4, const unsigned int n4) : - FiniteElement (multiply_dof_numbers(fe1, n1, - fe2, n2, - fe3, n3, - fe4, n4), + FiniteElement (multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4), compute_restriction_is_additive_flags (fe1, n1, fe2, n2, fe3, n3, @@ -240,22 +247,19 @@ FESystem::FESystem (const FiniteElement &fe1, fe2, n2, fe3, n3, fe4 ,n4)), - base_elements(4) + base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0)) { - this->base_to_block_indices.reinit(0, 0); - this->base_to_block_indices.push_back(n1); - this->base_to_block_indices.push_back(n2); - this->base_to_block_indices.push_back(n3); - this->base_to_block_indices.push_back(n4); - base_elements[0] = ElementPair(fe1.clone(), n1); - base_elements[0].first->subscribe (typeid(*this).name()); - base_elements[1] = ElementPair(fe2.clone(), n2); - base_elements[1].first->subscribe (typeid(*this).name()); - base_elements[2] = ElementPair(fe3.clone(), n3); - base_elements[2].first->subscribe (typeid(*this).name()); - base_elements[3] = ElementPair(fe4.clone(), n4); - base_elements[3].first->subscribe (typeid(*this).name()); - initialize (); + std::vector*> fes; + fes.push_back(&fe1); + fes.push_back(&fe2); + fes.push_back(&fe3); + fes.push_back(&fe4); + std::vector multiplicities; + multiplicities.push_back(n1); + multiplicities.push_back(n2); + multiplicities.push_back(n3); + multiplicities.push_back(n4); + initialize(fes, multiplicities); } template @@ -269,11 +273,11 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n4, const FiniteElement &fe5, const unsigned int n5) : - FiniteElement (multiply_dof_numbers(fe1, n1, - fe2, n2, - fe3, n3, - fe4, n4, - fe5, n5), + FiniteElement (multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4, + &fe5, n5), compute_restriction_is_additive_flags (fe1, n1, fe2, n2, fe3, n3, @@ -284,26 +288,19 @@ FESystem::FESystem (const FiniteElement &fe1, fe3, n3, fe4 ,n4, fe5, n5)), - base_elements(5) + base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0)+(n5>0)) { - this->base_to_block_indices.reinit(0, 0); - this->base_to_block_indices.push_back(n1); - this->base_to_block_indices.push_back(n2); - this->base_to_block_indices.push_back(n3); - this->base_to_block_indices.push_back(n4); - this->base_to_block_indices.push_back(n5); - - base_elements[0] = ElementPair(fe1.clone(), n1); - base_elements[0].first->subscribe (typeid(*this).name()); - base_elements[1] = ElementPair(fe2.clone(), n2); - base_elements[1].first->subscribe (typeid(*this).name()); - base_elements[2] = ElementPair(fe3.clone(), n3); - base_elements[2].first->subscribe (typeid(*this).name()); - base_elements[3] = ElementPair(fe4.clone(), n4); - base_elements[3].first->subscribe (typeid(*this).name()); - base_elements[4] = ElementPair(fe5.clone(), n5); - base_elements[4].first->subscribe (typeid(*this).name()); - initialize (); + std::vector*> fes; + fes.push_back(&fe1); + fes.push_back(&fe2); + fes.push_back(&fe3); + fes.push_back(&fe4); + std::vector multiplicities; + multiplicities.push_back(n1); + multiplicities.push_back(n2); + multiplicities.push_back(n3); + multiplicities.push_back(n4); + initialize(fes, multiplicities); } @@ -316,28 +313,12 @@ FESystem::FESystem ( FiniteElement (multiply_dof_numbers(fes, multiplicities), compute_restriction_is_additive_flags (fes, multiplicities), compute_nonzero_components(fes, multiplicities)), - base_elements(fes.size()) + base_elements(count_nonzeros(multiplicities)) { - Assert (fes.size() == multiplicities.size(), - ExcDimensionMismatch (fes.size(), multiplicities.size()) ); - Assert (fes.size() > 0, - ExcMessage ("Need to pass at least one finite element.")); - - this->base_to_block_indices.reinit(0, 0); - - for (unsigned int i=0; ibase_to_block_indices.push_back( multiplicities[i] ); - - for (unsigned int i=0; iclone(), multiplicities[i]); - base_elements[i].first->subscribe (typeid(*this).name()); - } - initialize (); + initialize(fes, multiplicities); } - template FESystem::~FESystem () { @@ -1872,8 +1853,37 @@ void FESystem::build_interface_constraints () template -void FESystem::initialize () +void FESystem::initialize (const std::vector*> &fes, + const std::vector &multiplicities) { + Assert (fes.size() == multiplicities.size(), + ExcDimensionMismatch (fes.size(), multiplicities.size()) ); + Assert (fes.size() > 0, + ExcMessage ("Need to pass at least one finite element.")); + Assert (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 + + this->base_to_block_indices.reinit(0, 0); + + for (unsigned int i=0; i0) + this->base_to_block_indices.push_back( multiplicities[i] ); + + unsigned int ind=0; + for (unsigned int i=0; i0) + { + base_elements[ind] = ElementPair(fes[i]->clone(), multiplicities[i]); + base_elements[ind].first->subscribe (typeid(*this).name()); + ++ind; + } + } + + Assert(ind>0, ExcInternalError()); + build_cell_tables(); build_face_tables(); @@ -1985,24 +1995,6 @@ void FESystem::initialize () -template -FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElementData &fe_data, - const unsigned int N) -{ - std::vector dpo; - dpo.push_back(fe_data.dofs_per_vertex * N); - dpo.push_back(fe_data.dofs_per_line * N); - 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, - fe_data.tensor_degree(), - fe_data.conforming_space, - N); -} - @@ -2640,167 +2632,31 @@ compare_for_face_domination (const FiniteElement &fe_other) const template FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2) +FESystem::multiply_dof_numbers (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2, + const unsigned int N2, + const FiniteElement *fe3, + const unsigned int N3, + const FiniteElement *fe4, + const unsigned int N4, + const FiniteElement *fe5, + const unsigned int N5) { - std::vector dpo; - dpo.push_back(fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2); - 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, - degree, - typename FiniteElementData::Conformity(fe1.conforming_space - & fe2.conforming_space), - N1 + N2); -} + std::vector*> fes; + fes.push_back(fe1); + fes.push_back(fe2); + fes.push_back(fe3); + fes.push_back(fe4); + fes.push_back(fe5); - - -template -FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3) -{ - std::vector dpo; - dpo.push_back(fe1.dofs_per_vertex * N1 + - fe2.dofs_per_vertex * N2 + - fe3.dofs_per_vertex * N3); - dpo.push_back(fe1.dofs_per_line * N1 + - fe2.dofs_per_line * N2 + - fe3.dofs_per_line * N3); - if (dim>1) dpo.push_back(fe1.dofs_per_quad * N1 + - fe2.dofs_per_quad * N2 + - fe3.dofs_per_quad * N3); - 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 - const unsigned int - degree = std::max (std::max(fe1.tensor_degree(), - fe2.tensor_degree()), - fe3.tensor_degree()); - return FiniteElementData ( - dpo, - fe1.n_components() * N1 + fe2.n_components() * N2 + fe3.n_components() * N3, - degree, - typename FiniteElementData::Conformity(fe1.conforming_space - & fe2.conforming_space - & fe3.conforming_space), - N1 + N2 + N3); -} - - - -template -FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3, - const FiniteElementData &fe4, - const unsigned int N4) -{ - std::vector dpo; - dpo.push_back(fe1.dofs_per_vertex * N1 + - fe2.dofs_per_vertex * N2 + - fe3.dofs_per_vertex * N3 + - fe4.dofs_per_vertex * N4); - dpo.push_back(fe1.dofs_per_line * N1 + - fe2.dofs_per_line * N2 + - fe3.dofs_per_line * N3 + - fe4.dofs_per_line * N4); - if (dim>1) dpo.push_back(fe1.dofs_per_quad * N1 + - fe2.dofs_per_quad * N2 + - fe3.dofs_per_quad * N3 + - fe4.dofs_per_quad * N4); - if (dim>2) dpo.push_back(fe1.dofs_per_hex * N1 + - fe2.dofs_per_hex * N2 + - fe3.dofs_per_hex * N3 + - fe4.dofs_per_hex * N4); - // degree is the maximal degree of the components - const unsigned int - degree = std::max(std::max (std::max(fe1.tensor_degree(),fe2.tensor_degree()),fe3.tensor_degree()), fe4.tensor_degree()); - - return FiniteElementData ( - dpo, - fe1.n_components() * N1 + fe2.n_components() * N2 + fe3.n_components() * N3 + - fe4.n_components() * N4, - degree, - typename FiniteElementData::Conformity(fe1.conforming_space - & fe2.conforming_space - & fe3.conforming_space - & fe4.conforming_space), - N1 + N2 + N3 + N4); -} - - -template -FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElementData &fe1, - const unsigned int N1, - const FiniteElementData &fe2, - const unsigned int N2, - const FiniteElementData &fe3, - const unsigned int N3, - const FiniteElementData &fe4, - const unsigned int N4, - const FiniteElementData &fe5, - const unsigned int N5) -{ - std::vector dpo; - dpo.push_back(fe1.dofs_per_vertex * N1 + - fe2.dofs_per_vertex * N2 + - fe3.dofs_per_vertex * N3 + - fe4.dofs_per_vertex * N4 + - fe5.dofs_per_vertex * N5); - dpo.push_back(fe1.dofs_per_line * N1 + - fe2.dofs_per_line * N2 + - fe3.dofs_per_line * N3 + - fe4.dofs_per_line * N4 + - fe5.dofs_per_line * N5); - if (dim>1) dpo.push_back(fe1.dofs_per_quad * N1 + - fe2.dofs_per_quad * N2 + - fe3.dofs_per_quad * N3 + - fe4.dofs_per_quad * N4 + - fe5.dofs_per_quad * N5); - if (dim>2) dpo.push_back(fe1.dofs_per_hex * N1 + - fe2.dofs_per_hex * N2 + - fe3.dofs_per_hex * N3 + - fe4.dofs_per_hex * N4 + - fe5.dofs_per_hex * N5); - // degree is the maximal degree of the components - const unsigned int - degree = std::max(std::max(std::max (std::max(fe1.tensor_degree(),fe2.tensor_degree()),fe3.tensor_degree()), fe4.tensor_degree()), fe5.tensor_degree()); - - return FiniteElementData ( - dpo, - fe1.n_components() * N1 + fe2.n_components() * N2 + fe3.n_components() * N3 + - fe4.n_components() * N4 + fe5.n_components() * N5, - degree, - typename FiniteElementData::Conformity(fe1.conforming_space - & fe2.conforming_space - & fe3.conforming_space - & fe4.conforming_space - & fe5.conforming_space), - N1 + N2 + N3 + N4 + N5); + std::vector mult; + mult.push_back(N1); + mult.push_back(N2); + mult.push_back(N3); + mult.push_back(N4); + mult.push_back(N5); + return multiply_dof_numbers(fes, mult); } @@ -2825,24 +2681,25 @@ FESystem::multiply_dof_numbers ( unsigned int summed_multiplicities = 0; for (unsigned int i=0; idofs_per_vertex * multiplicities[i]; - multiplied_dofs_per_line+=fes[i]->dofs_per_line * multiplicities[i]; - multiplied_dofs_per_quad+=fes[i]->dofs_per_quad * multiplicities[i]; - multiplied_dofs_per_hex+=fes[i]->dofs_per_hex * multiplicities[i]; - - multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; + if (multiplicities[i]>0) + { + multiplied_dofs_per_vertex+=fes[i]->dofs_per_vertex * multiplicities[i]; + multiplied_dofs_per_line+=fes[i]->dofs_per_line * multiplicities[i]; + multiplied_dofs_per_quad+=fes[i]->dofs_per_quad * multiplicities[i]; + multiplied_dofs_per_hex+=fes[i]->dofs_per_hex * multiplicities[i]; - degree = std::max(degree, fes[i]->tensor_degree() ); + multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; - summed_multiplicities += multiplicities[i]; - } + degree = std::max(degree, fes[i]->tensor_degree() ); - unsigned char total_conformity = fes[0]->conforming_space; + summed_multiplicities += multiplicities[i]; + } - for (unsigned int i=1; iconforming_space; + unsigned char total_conformity = 0xff;; + for (unsigned int i=0; i0) + total_conformity &= fes[i]->conforming_space; std::vector dpo; dpo.push_back(multiplied_dofs_per_vertex);