]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Jason Sheldon: Allow the construction of FESystem objects with arbitrarily...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 27 Jan 2012 23:34:40 +0000 (23:34 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 27 Jan 2012 23:34:40 +0000 (23:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@24945 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/fe_system.h
deal.II/source/fe/fe_system.cc
deal.II/source/fe/fe_tools.cc

index 6ab07c62fd7cf0169e90489e4a26e090b27c0667..d5d7e084b87e15be92def7690fd4d7fd88c21ae4 100644 (file)
@@ -81,9 +81,13 @@ enabled due to a missing include file in file
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Improved: There is now a constructor for FESystem that allows to
+create collections of finite elements of arbitrary length.
+<br>
+(Jason Sheldon, 2012/01/27)
 
-<li> Improved: <code>VectorTools::point_value</code> now also works within the hp framework.
-Fixed: <code>GridTools::find_active_cell_around_point</code> for the hp-case works now also with MappingCollections containing only one mapping, as is the standard case in other functions using hp.
+<li> 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.
 <br>
 (Christian Goll, 2012/01/26)
 
index 38c7109688222b31338118502f3d87ab5f1ec1ff..65271cc48777dcedf3e0a99047c7ba629e0ae73c 100644 (file)
@@ -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<dim,spacedim>
              const FiniteElement<dim,spacedim> &fe4, const unsigned int n4,
              const FiniteElement<dim,spacedim> &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<const FiniteElement<dim,spacedim>*> &fes,
+             const std::vector<unsigned int>                   &multiplicities);
+
                                     /**
                                      * Destructor.
                                      */
@@ -424,7 +439,7 @@ class FESystem : public FiniteElement<dim,spacedim>
     virtual void
     get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                              FullMatrix<double>           &matrix) const;
-    
+
                                     /**
                                      * Access to a composing
                                      * element. The index needs to be
@@ -889,6 +904,16 @@ class FESystem : public FiniteElement<dim,spacedim>
                          const FiniteElementData<dim> &fe5,
                          const unsigned int            N5);
 
+                                  /**
+                                   * Same as above but for
+                                   * any number of sub-elements.
+                                   */
+    static FiniteElementData<dim>
+    multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+                         const std::vector<unsigned int>                       &multiplicities);
+
+
+
                                     /**
                                      * Helper function used in the
                                      * constructor: takes a
index 054846bf93b395261028f5f31bf62cec3fb030dc..d2f015ee48e149f0ca6e07778e7e865e954a9f7e 100644 (file)
@@ -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<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &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<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe1,
 }
 
 
+
+template <int dim, int spacedim>
+FESystem<dim,spacedim>::FESystem (
+  const std::vector<const FiniteElement<dim,spacedim>*>  &fes,
+  const std::vector<unsigned int>                  &multiplicities)
+               :
+               FiniteElement<dim,spacedim> (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; i<fes.size(); i++)
+    this->base_to_block_indices.push_back( multiplicities[i] );
+
+  for (unsigned int i=0; i<fes.size(); i++)
+    {
+      base_elements[i] = ElementPair(fes[i]->clone(), multiplicities[i]);
+      base_elements[i].first->subscribe (typeid(*this).name());
+    }
+  initialize ();
+}
+
+
+
 template <int dim, int spacedim>
 FESystem<dim,spacedim>::~FESystem ()
 {
@@ -1129,7 +1158,7 @@ FESystem<dim,spacedim>::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<dim,spacedim>::compute_fill (
                        for (unsigned int q=0; q<n_q_points; ++q)
                          data.shape_values[out_index+s][q] =
                            base_data.shape_values(in_index+s,q);
-                     
+
                      if (base_flags & update_gradients)
                        for (unsigned int q=0; q<n_q_points; ++q)
                          data.shape_gradients[out_index+s][q]=
@@ -1245,7 +1274,7 @@ FESystem<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::multiply_dof_numbers (const FiniteElementData<dim> &fe1,
 }
 
 
+
+template <int dim, int spacedim>
+FiniteElementData<dim>
+FESystem<dim,spacedim>::multiply_dof_numbers (
+  const std::vector<const FiniteElement<dim,spacedim>*>   &fes,
+  const std::vector<unsigned int>                   &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; i<fes.size(); i++)
+    {
+      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];
+
+      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; i<fes.size(); i++)
+    total_conformity &= fes[i]->conforming_space;
+
+
+  std::vector<unsigned int> 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<dim> (
+    dpo,
+    multiplied_n_components,
+    degree,
+    typename FiniteElementData<dim>::Conformity(total_conformity),
+    summed_multiplicities);
+}
+
+
+
 template <int dim, int spacedim>
 std::vector<bool>
 FESystem<dim,spacedim>::compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> &fe,
index c80ee3495747740ca84f2d4d867d3131b265972b..efd94228452f11cd2cf127d704d86a58d03ab4f5 100644 (file)
@@ -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<FiniteElement<dim,spacedim>*> base_fes;
+           std::vector<const FiniteElement<dim,spacedim>*> base_fes;
            std::vector<unsigned int>        base_multiplicities;
            try
              {
@@ -2112,38 +2116,12 @@ namespace FETools
                                                 // generate the composed
                                                 // element
                FiniteElement<dim,spacedim> *system_element = 0;
-               switch (base_fes.size())
-                 {
-                   case 1:
-                   {
-                     system_element = new FESystem<dim>(*base_fes[0],
-                                                        base_multiplicities[0]);
-                     break;
-                   }
-
-                   case 2:
-                   {
-                     system_element = new FESystem<dim>(*base_fes[0],
-                                                        base_multiplicities[0],
-                                                        *base_fes[1],
-                                                        base_multiplicities[1]);
-                     break;
-                   }
-
-                   case 3:
-                   {
-                     system_element = new FESystem<dim>(*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<dim>(base_fes,
+                                       base_multiplicities);
 
                                                 // now we don't need the
                                                 // list of base elements

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.