From 594ec5cff95f644a91580e4ab3ed9b228dfe62d3 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 27 Jan 2012 23:34:40 +0000 Subject: [PATCH] Patch by Jason Sheldon: Allow the construction of FESystem objects with arbitrarily many base elements. git-svn-id: https://svn.dealii.org/trunk@24945 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 8 ++- deal.II/include/deal.II/fe/fe_system.h | 29 +++++++- deal.II/source/fe/fe_system.cc | 97 ++++++++++++++++++++++++-- deal.II/source/fe/fe_tools.cc | 44 +++--------- 4 files changed, 135 insertions(+), 43 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 6ab07c62fd..d5d7e084b8 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -81,9 +81,13 @@ enabled due to a missing include file in file

Specific improvements

    +
  1. Improved: There is now a constructor for FESystem that allows to +create collections of finite elements of arbitrary length. +
    +(Jason Sheldon, 2012/01/27) -
  2. Improved: VectorTools::point_value now also works within the hp framework. -Fixed: GridTools::find_active_cell_around_point for the hp-case works now also with MappingCollections containing only one mapping, as is the standard case in other functions using hp. +
  3. Improved: VectorTools::point_value() now also works within the hp framework. +Fixed: GridTools::find_active_cell_around_point() for the hp-case works now also with MappingCollections containing only one mapping, as is the standard case in other functions using hp.
    (Christian Goll, 2012/01/26) diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index 38c7109688..65271cc487 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -220,6 +220,21 @@ class FESystem : public FiniteElement const FiniteElement &fe4, const unsigned int n4, const FiniteElement &fe5, const unsigned int n5); + /** + * Same as above but for any + * number of base + * elements. Pointers to the base + * elements and their + * multiplicities are passed as + * vectors to this + * constructor. The length of + * these vectors is assumed to be + * equal. + */ + + FESystem (const std::vector*> &fes, + const std::vector &multiplicities); + /** * Destructor. */ @@ -424,7 +439,7 @@ class FESystem : public FiniteElement virtual void get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; - + /** * Access to a composing * element. The index needs to be @@ -889,6 +904,16 @@ class FESystem : public FiniteElement const FiniteElementData &fe5, const unsigned int N5); + /** + * Same as above but for + * any number of sub-elements. + */ + static FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities); + + + /** * Helper function used in the * constructor: takes a diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 054846bf93..d2f015ee48 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -292,7 +292,7 @@ FESystem::FESystem (const FiniteElement &fe1, 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); @@ -307,6 +307,35 @@ FESystem::FESystem (const FiniteElement &fe1, } + +template +FESystem::FESystem ( + const std::vector*> &fes, + const std::vector &multiplicities) + : + FiniteElement (multiply_dof_numbers(fes, multiplicities), + compute_restriction_is_additive_flags (fes, multiplicities), + compute_nonzero_components(fes, multiplicities)), + base_elements(fes.size()) +{ + Assert (fes.size() == multiplicities.size(), + ExcDimensionMismatch (fes.size(), multiplicities.size()) ); + + 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 (); +} + + + template FESystem::~FESystem () { @@ -1129,7 +1158,7 @@ FESystem::compute_fill ( // those shape functions // that belong to a given // base element -//TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style +//TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style const UpdateFlags base_flags(dim_1==dim ? base_fe_data.current_update_flags() : base_fe_data.update_flags); @@ -1183,7 +1212,7 @@ FESystem::compute_fill ( for (unsigned int q=0; q::build_cell_tables() unsigned total_index = 0; this->block_indices_data.reinit(0,0); - + for (unsigned int base=0; base < this->n_base_elements(); ++base) for (unsigned int m = 0; m < this->element_multiplicity(base); ++m) { @@ -1256,7 +1285,7 @@ FESystem::build_cell_tables() } Assert (total_index == this->component_to_base_table.size(), ExcInternalError()); - + // Initialize index tables. // Multi-component base elements // have to be thought of. For @@ -2805,6 +2834,62 @@ FESystem::multiply_dof_numbers (const FiniteElementData &fe1, } + +template +FiniteElementData +FESystem::multiply_dof_numbers ( + const std::vector*> &fes, + const std::vector &multiplicities) +{ + Assert (fes.size() == multiplicities.size(), ExcDimensionMismatch (fes.size(), multiplicities.size())); + + unsigned int multiplied_dofs_per_vertex = 0; + unsigned int multiplied_dofs_per_line = 0; + unsigned int multiplied_dofs_per_quad = 0; + unsigned int multiplied_dofs_per_hex = 0; + + unsigned int multiplied_n_components = 0; + + unsigned int degree = 0; // degree is the maximal degree of the components + + 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]; + + degree = std::max(degree, fes[i]->tensor_degree() ); + + summed_multiplicities += multiplicities[i]; + } + + unsigned char total_conformity = fes[0]->conforming_space; + + for (unsigned int i=1; iconforming_space; + + + std::vector dpo; + dpo.push_back(multiplied_dofs_per_vertex); + dpo.push_back(multiplied_dofs_per_line); + if (dim>1) dpo.push_back(multiplied_dofs_per_quad); + if (dim>2) dpo.push_back(multiplied_dofs_per_hex); + + return FiniteElementData ( + dpo, + multiplied_n_components, + degree, + typename FiniteElementData::Conformity(total_conformity), + summed_multiplicities); +} + + + template std::vector FESystem::compute_restriction_is_additive_flags (const FiniteElement &fe, diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc index c80ee34957..efd9422845 100644 --- a/deal.II/source/fe/fe_tools.cc +++ b/deal.II/source/fe/fe_tools.cc @@ -2027,6 +2027,10 @@ namespace FETools // have to figure out what the // base elements are. this can // only be done recursively + + + + if (name_part == "FESystem") { // next we have to get at the @@ -2038,7 +2042,7 @@ namespace FETools // recursive calls if one of // these calls should throw // an exception - std::vector*> base_fes; + std::vector*> base_fes; std::vector base_multiplicities; try { @@ -2112,38 +2116,12 @@ namespace FETools // generate the composed // element FiniteElement *system_element = 0; - switch (base_fes.size()) - { - case 1: - { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0]); - break; - } - - case 2: - { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0], - *base_fes[1], - base_multiplicities[1]); - break; - } - - case 3: - { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0], - *base_fes[1], - base_multiplicities[1], - *base_fes[2], - base_multiplicities[2]); - break; - } - - default: - AssertThrow (false, ExcNotImplemented()); - } + + // uses new FESystem constructor + // which is independent of + // the number of FEs in the system + system_element = new FESystem(base_fes, + base_multiplicities); // now we don't need the // list of base elements -- 2.39.5