]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Change the meaning of the restriction_is_additive flag: previously it was indexed...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 16 Sep 2002 13:29:37 +0000 (13:29 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 16 Sep 2002 13:29:37 +0000 (13:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@6428 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_base.h
deal.II/deal.II/include/fe/fe_system.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_dgp.cc
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/fe_system.cc

index 6d2312c7f58ca9ffb70784d329356f706f02c3ff..e30daa0172d8b73ee37115b691ec1833900a95f8 100644 (file)
@@ -409,11 +409,17 @@ class FiniteElementBase : public Subscriptor,
   
                                     /**
                                      * Construct an object of this
-                                     * type.  You have to set some
+                                     * type. You have to set some
                                      * member variables, for example
                                      * some matrices, explicitly
                                      * after calling this base class'
-                                     * constructor.
+                                     * constructor. For this see the
+                                     * existing finite element
+                                     * classes. For the second and
+                                     * third parameter of this
+                                     * constructor, see the document
+                                     * of the respective member
+                                     * variables.
                                      */
     FiniteElementBase (const FiniteElementData<dim> &fe_data,
                       const std::vector<bool> &restriction_is_additive_flags,
@@ -816,11 +822,16 @@ class FiniteElementBase : public Subscriptor,
     component_to_base (unsigned int component) const;
 
                                     /**
-                                     * Access the @p{restriction_is_additive_flag}
-                                     * field. See there for more information on 
-                                     * its contents.
+                                     * Access the
+                                     * @p{restriction_is_additive_flag}
+                                     * field. See there for more
+                                     * information on its contents.
+                                     *
+                                     * The index must be between zero
+                                     * and the number of shape
+                                     * functions of this element.
                                      */
-    bool restriction_is_additive (const unsigned int component) const;
+    bool restriction_is_additive (const unsigned int index) const;
 
                                     /**
                                      * Return the support points of
@@ -1187,18 +1198,20 @@ class FiniteElementBase : public Subscriptor,
                                      *
                                      * This array has valid values
                                      * also in the case of
-                                     * vector-value
+                                     * vector-valued
                                      * (i.e. non-primitive) shape
                                      * functions, in contrast to the
                                      * @p{system_to_component_table}.
                                      */
-    std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> > system_to_base_table;
+    std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
+    system_to_base_table;
 
                                     /**
                                      * Likewise for the indices on
                                      * faces.
                                      */
-    std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> > face_system_to_base_table;
+    std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
+    face_system_to_base_table;
     
                                     /**
                                      * Map between component and
@@ -1278,9 +1291,48 @@ class FiniteElementBase : public Subscriptor,
                                      * build the complete matrix. The
                                      * flag should be @p{true}.
                                      *
-                                     * There is one flag per
-                                     * component in vector valued
-                                     * elements.
+                                     * For examples of use of these
+                                     * flags, see the places in the
+                                     * library where it is queried.
+                                     * 
+                                     * There is one flag per shape
+                                     * function, indicating whether
+                                     * it belongs to the class of
+                                     * shape functions that are
+                                     * additive in the restriction or
+                                     * not.
+                                     *
+                                     * Note that in previous versions
+                                     * of the library, there was one
+                                     * flag per vector component of
+                                     * the element. This is based on
+                                     * the fact that all the shape
+                                     * functions that belong to the
+                                     * same vector component must
+                                     * necessarily behave in the same
+                                     * way, to make things
+                                     * reasonable. However, the
+                                     * problem is that it is
+                                     * sometimes impossible to query
+                                     * this flag in the vector-valued
+                                     * case: this used to be done
+                                     * with the
+                                     * @p{system_to_component_index}
+                                     * function that returns which
+                                     * vector component a shape
+                                     * function is associated
+                                     * with. The point is that since
+                                     * we now support shape functions
+                                     * that are associated with more
+                                     * than one vector component (for
+                                     * example the shape functions of
+                                     * Raviart-Thomas, or Nedelec
+                                     * elements), that function can
+                                     * no more be used, so it can be
+                                     * difficult to find out which
+                                     * for vector component we would
+                                     * like to query the
+                                     * restriction-is-additive flags.
                                      */
     const std::vector<bool> restriction_is_additive_flags;
 
@@ -1558,11 +1610,11 @@ FiniteElementBase<dim>::component_to_base (unsigned int index) const
 template <int dim>
 inline
 bool
-FiniteElementBase<dim>::restriction_is_additive (const unsigned int component) const
+FiniteElementBase<dim>::restriction_is_additive (const unsigned int index) const
 {
-  Assert(component<this->n_components(),
-        ExcIndexRange(component, 0, this->n_components()));
-  return restriction_is_additive_flags[component];
+  Assert(index < this->dofs_per_cell,
+        ExcIndexRange(index, 0, this->dofs_per_cell));
+  return restriction_is_additive_flags[index];
 }
 
 
index 04f164a91ffce324b90f22132056d13ba25199e3..1247420c46151b9257a04497eb89d6e4599877a7 100644 (file)
@@ -525,12 +525,15 @@ class FESystem : public FiniteElement<dim>
 
 
                                     /**
-                                     * Helper function used in the constructor:
-                                     * takes a @p{FiniteElement} object
-                                     * and returns an boolean vector including
-                                     * the @p{restriction_is_additive_flags} of
-                                     * the mixed element consisting of @p{N}
-                                     * elements of the sub-element @p{fe}.
+                                     * Helper function used in the
+                                     * constructor: takes a
+                                     * @p{FiniteElement} object and
+                                     * returns an boolean vector
+                                     * including the
+                                     * @p{restriction_is_additive_flags}
+                                     * of the mixed element
+                                     * consisting of @p{N} elements
+                                     * of the sub-element @p{fe}.
                                      */
     static std::vector<bool>
     compute_restriction_is_additive_flags (const FiniteElement<dim> &fe,
@@ -558,6 +561,19 @@ class FESystem : public FiniteElement<dim>
                                           const FiniteElement<dim> &fe3,
                                           const unsigned int        N3);
 
+                                    /**
+                                     * Compute the named flags for a
+                                     * list of finite elements with
+                                     * multiplicities given in the
+                                     * second argument. This function
+                                     * is called from all the above
+                                     * functions.
+                                    */
+    static std::vector<bool>
+    compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim>*> &fes,
+                                           const std::vector<unsigned int>              &multiplicities);
+    
+    
                                     /**
                                      * Compute the non-zero vector
                                      * components of a composed
index 6002b08e73bdf70dbaad46b0e13641ebcb564783..d0b1c0d7b9650dc9a61e4b4f55e61df989fce5cb 100644 (file)
@@ -108,34 +108,35 @@ FiniteElementBase<dim>::InternalDataBase::~InternalDataBase ()
 
 
 template <int dim>
-FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data,
-                                          const std::vector<bool> &restriction_is_additive_flags,
-                                          const std::vector<std::vector<bool> > &nonzero_components)
+FiniteElementBase<dim>::
+FiniteElementBase (const FiniteElementData<dim> &fe_data,
+                   const std::vector<bool> &restriction_is_additive_flags,
+                   const std::vector<std::vector<bool> > &nonzero_components)
                :
                FiniteElementData<dim> (fe_data),
-  system_to_component_table(this->dofs_per_cell),
-  face_system_to_component_table(this->dofs_per_face),
-  system_to_base_table(this->dofs_per_cell),
-  face_system_to_base_table(this->dofs_per_face),              
-  component_to_system_table(this->components,
-                           std::vector<unsigned>(this->dofs_per_cell)),
-                             face_component_to_system_table(this->components,
-                                                            std::vector<unsigned>(this->dofs_per_face)),
-                                                              component_to_base_table (this->components,
-                                                                                       std::make_pair(0U, 0U)),
-                                                              restriction_is_additive_flags(restriction_is_additive_flags),
-                                                              nonzero_components (nonzero_components),
-                                                              n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
-                                                              cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
-                                                                                                n_nonzero_components_table.end(),
-                                                                                                std::bind2nd(std::not_equal_to<unsigned int>(),
-                                                                                                             1U))
-                                                                                  ==
-                                                                                  n_nonzero_components_table.end())
+                system_to_component_table (this->dofs_per_cell),
+                face_system_to_component_table(this->dofs_per_face),
+                system_to_base_table(this->dofs_per_cell),
+                face_system_to_base_table(this->dofs_per_face),                
+                component_to_system_table(this->components,
+                                          std::vector<unsigned>(this->dofs_per_cell)),
+               face_component_to_system_table(this->components,
+                                               std::vector<unsigned>(this->dofs_per_face)),
+                component_to_base_table (this->components,
+                                         std::make_pair(0U, 0U)),
+                restriction_is_additive_flags(restriction_is_additive_flags),
+               nonzero_components (nonzero_components),
+               n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
+               cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
+                                                  n_nonzero_components_table.end(),
+                                                  std::bind2nd(std::not_equal_to<unsigned int>(),
+                                                               1U))
+                                    ==
+                                    n_nonzero_components_table.end())
 {
-  Assert (restriction_is_additive_flags.size()==fe_data.components,
+  Assert (restriction_is_additive_flags.size() == this->dofs_per_cell,
          ExcDimensionMismatch(restriction_is_additive_flags.size(),
-                              fe_data.components));
+                              this->dofs_per_cell));
   Assert (nonzero_components.size() == this->dofs_per_cell,
          ExcInternalError());
   for (unsigned int i=0; i<nonzero_components.size(); ++i)
index 94f9b9de693a4150220b31be1f3347de16df8968..b8331e09931cc75ee2d26cdf839c77d68b80f4e5 100644 (file)
@@ -26,11 +26,11 @@ template <int dim>
 FE_DGP<dim>::FE_DGP (const unsigned int degree)
                :
                FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
-                                   std::vector<bool>(1,true),
+                                   std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,true),
                                    std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
                                                                    std::vector<bool>(1,true))),
-                                                                     degree(degree),
-                                                                     polynomial_space (Legendre<double>::generate_complete_basis(degree))
+                degree(degree),
+                polynomial_space (Legendre<double>::generate_complete_basis(degree))
 {
                                   // if defined, copy over matrices
                                   // from precomputed arrays
@@ -41,7 +41,15 @@ FE_DGP<dim>::FE_DGP (const unsigned int degree)
   else
     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
       this->prolongation[i].reinit(0,0);
-  
+
+                                   // restriction can be defined
+                                   // through projection for
+                                   // discontinuous elements, but is
+                                   // presently not implemented for DGP
+                                   // elements.
+                                   //
+                                   // if it were, then the following
+                                   // snippet would be the right code
 //    if ((degree < Matrices::n_projection_matrices) &&
 //        (Matrices::projection_matrices[degree] != 0))
 //      {
@@ -50,7 +58,16 @@ FE_DGP<dim>::FE_DGP (const unsigned int degree)
 //    else
 //                                  // matrix undefined, set size to zero
 //      for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-//        restriction[i].reinit(0);  
+//        restriction[i].reinit(0, 0);
+                                   // since not implemented, set to
+                                   // "empty"
+  for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
+    restriction[i].reinit(0, 0);
+
+                                   // note further, that these
+                                   // elements have neither support
+                                   // nor face-support points, so
+                                   // leave these fields empty
 };
 
 
index a638f981b0e5dc408602fa604d9aa4d9b85d3b11..6ddd18b9301929cc09d85749d228a8df20f84d57 100644 (file)
@@ -28,7 +28,8 @@ template <int dim>
 FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
                :
                FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
-                                   std::vector<bool>(1,true),
+                                   std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+                                                      true),
                                    std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
                                                                    std::vector<bool>(1,true))),
                                                                      degree(degree),
@@ -47,6 +48,12 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
                                   // from precomputed arrays and
                                   // generate all other matrices by
                                   // permutations
+                                   //
+                                   // (note: the matrix is defined if
+                                   // something was entered into the
+                                   // respective table, and what was
+                                   // entered is not a NULL pointer --
+                                   // this would allow for "holes")
   if ((degree < Matrices::n_embedding_matrices) &&
       (Matrices::embedding[degree] != 0))
     {
index 192fc8f0ed3930962544902e4bddfc9e9e43c9ae..e20e49d2e3d30f3991a21c0adac0ec590d5be52d 100644 (file)
@@ -27,7 +27,8 @@ template <int dim>
 FE_Q<dim>::FE_Q (const unsigned int degree)
                :
                FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
-                                   std::vector<bool> (1,false),
+                                   std::vector<bool> (FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+                                                       false),
                                    std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
                                                                    std::vector<bool>(1,true))),
                                                                      degree(degree),
index 69b3ba4735231da1bfa1124f270d6243340b3730..9b8e526df9c4d97e52aa8ede9505b6fff88377ab 100644 (file)
@@ -1680,11 +1680,13 @@ std::vector<bool>
 FESystem<dim>::compute_restriction_is_additive_flags (const FiniteElement<dim> &fe,
                                                      const unsigned int n_elements) 
 {
-  std::vector<bool> tmp;
-  for (unsigned int i=0; i<n_elements; ++i)
-    for (unsigned int component=0; component<fe.n_components(); ++component)
-      tmp.push_back (fe.restriction_is_additive (component));
-  return tmp;
+  std::vector<const FiniteElement<dim>*> fe_list;
+  std::vector<unsigned int>              multiplicities;
+
+  fe_list.push_back (&fe);
+  multiplicities.push_back (n_elements);
+  
+  return compute_restriction_is_additive_flags (fe_list, multiplicities);
 };
 
 
@@ -1696,14 +1698,16 @@ FESystem<dim>::compute_restriction_is_additive_flags (const FiniteElement<dim> &
                                                      const FiniteElement<dim> &fe2,
                                                      const unsigned int        N2) 
 {
-  std::vector<bool> tmp;
-  for (unsigned int i=0; i<N1; ++i)
-    for (unsigned int component=0; component<fe1.n_components(); ++component)
-      tmp.push_back (fe1.restriction_is_additive (component));
-  for (unsigned int i=0; i<N2; ++i)
-    for (unsigned int component=0; component<fe2.n_components(); ++component)
-      tmp.push_back (fe2.restriction_is_additive (component));
-  return tmp;
+  std::vector<const FiniteElement<dim>*> fe_list;
+  std::vector<unsigned int>              multiplicities;
+
+  fe_list.push_back (&fe1);
+  multiplicities.push_back (N1);
+
+  fe_list.push_back (&fe2);
+  multiplicities.push_back (N2);
+  
+  return compute_restriction_is_additive_flags (fe_list, multiplicities);
 };
 
 
@@ -1717,17 +1721,151 @@ FESystem<dim>::compute_restriction_is_additive_flags (const FiniteElement<dim> &
                                                      const FiniteElement<dim> &fe3,
                                                      const unsigned int        N3) 
 {
-  std::vector<bool> tmp;
-  for (unsigned int i=0; i<N1; ++i)
-    for (unsigned int component=0; component<fe1.n_components(); ++component)
-      tmp.push_back (fe1.restriction_is_additive (component));
-  for (unsigned int i=0; i<N2; ++i)
-    for (unsigned int component=0; component<fe2.n_components(); ++component)
-      tmp.push_back (fe2.restriction_is_additive (component));
-  for (unsigned int i=0; i<N3; ++i)
-    for (unsigned int component=0; component<fe3.n_components(); ++component)
-      tmp.push_back (fe3.restriction_is_additive (component));
-  return tmp;
+  std::vector<const FiniteElement<dim>*> fe_list;
+  std::vector<unsigned int>              multiplicities;
+
+  fe_list.push_back (&fe1);
+  multiplicities.push_back (N1);
+
+  fe_list.push_back (&fe2);
+  multiplicities.push_back (N2);
+
+  fe_list.push_back (&fe3);
+  multiplicities.push_back (N3);
+  
+  return compute_restriction_is_additive_flags (fe_list, multiplicities);
+};
+
+
+
+template <int dim>
+std::vector<bool>
+FESystem<dim>::
+compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim>*> &fes,
+                                       const std::vector<unsigned int>              &multiplicities)
+{
+  Assert (fes.size() == multiplicities.size(), ExcInternalError());
+
+                                  // first count the number of dofs
+                                  // and components that will emerge
+                                  // from the given FEs
+  unsigned int n_shape_functions = 0;
+  for (unsigned int i=0; i<fes.size(); ++i)
+    n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
+
+                                  // generate the array that will
+                                  // hold the output
+  std::vector<bool> retval (n_shape_functions, false);
+
+                                  // finally go through all the shape
+                                  // functions of the base elements,
+                                  // and copy their flags. this
+                                  // somehow copies the code in
+                                  // build_cell_table, which is not
+                                  // nice as it uses too much
+                                  // implicit knowledge about the
+                                  // layout of the individual bases
+                                  // in the composed FE, but there
+                                  // seems no way around...
+                                  //
+                                  // for each shape function, copy
+                                  // the flags from the base element
+                                  // to this one, taking into account
+                                  // multiplicities, and other
+                                  // complications
+  unsigned int total_index = 0;
+  for (unsigned int vertex_number=0;
+       vertex_number<GeometryInfo<dim>::vertices_per_cell;
+       ++vertex_number)
+    {
+      for(unsigned int base=0; base<fes.size(); ++base)
+       for (unsigned int m=0; m<multiplicities[base]; ++m)
+         for (unsigned int local_index = 0;
+              local_index < fes[base]->dofs_per_vertex;
+              ++local_index, ++total_index)
+           {
+             const unsigned int index_in_base
+               = (fes[base]->dofs_per_vertex*vertex_number + 
+                  local_index);
+
+              Assert (index_in_base < fes[base]->dofs_per_cell,
+                      ExcInternalError());
+              retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+           }
+    }
+  
+                                  // 2. Lines
+  if (GeometryInfo<dim>::lines_per_cell > 0)
+    for (unsigned int line_number= 0;
+        line_number != GeometryInfo<dim>::lines_per_cell;
+        ++line_number)
+      {
+       for (unsigned int base=0; base<fes.size(); ++base)
+         for (unsigned int m=0; m<multiplicities[base]; ++m)
+           for (unsigned int local_index = 0;
+                local_index < fes[base]->dofs_per_line;
+                ++local_index, ++total_index)
+             {
+               const unsigned int index_in_base
+                 = (fes[base]->dofs_per_line*line_number + 
+                    local_index +
+                    fes[base]->first_line_index);
+
+                Assert (index_in_base < fes[base]->dofs_per_cell,
+                        ExcInternalError());
+                retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+             }
+      }
+  
+                                  // 3. Quads
+  if (GeometryInfo<dim>::quads_per_cell > 0)
+    for (unsigned int quad_number= 0;
+        quad_number != GeometryInfo<dim>::quads_per_cell;
+        ++quad_number)
+      {
+       for (unsigned int base=0; base<fes.size(); ++base)
+         for (unsigned int m=0; m<multiplicities[base]; ++m)
+           for (unsigned int local_index = 0;
+                local_index < fes[base]->dofs_per_quad;
+                ++local_index, ++total_index)
+             {
+               const unsigned int index_in_base
+                 = (fes[base]->dofs_per_quad*quad_number + 
+                    local_index +
+                    fes[base]->first_quad_index);
+
+                Assert (index_in_base < fes[base]->dofs_per_cell,
+                        ExcInternalError());
+                retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+             }
+      }
+  
+                                  // 4. Hexes
+  if (GeometryInfo<dim>::hexes_per_cell > 0)
+    for (unsigned int hex_number= 0;
+        hex_number != GeometryInfo<dim>::hexes_per_cell;
+        ++hex_number)
+      {
+       for(unsigned int base=0; base<fes.size(); ++base)
+         for (unsigned int m=0; m<multiplicities[base]; ++m)
+           for (unsigned int local_index = 0;
+                local_index < fes[base]->dofs_per_hex;
+                ++local_index, ++total_index)
+             {
+               const unsigned int index_in_base
+                 = (fes[base]->dofs_per_hex*hex_number + 
+                    local_index +
+                    fes[base]->first_hex_index);
+
+                Assert (index_in_base < fes[base]->dofs_per_cell,
+                        ExcInternalError());
+                retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+             }
+      }
+
+  Assert (total_index == n_shape_functions, ExcInternalError());
+  
+  return retval;
 };
 
 
@@ -1970,13 +2108,14 @@ compute_nonzero_components (const std::vector<const FiniteElement<dim>*> &fes,
              }
       }
 
-  Assert (total_index == retval.size(), ExcInternalError());
+  Assert (total_index == n_shape_functions, ExcInternalError());
   
   return retval;
 };
 
 
 
+
 template <int dim>
 const FiniteElement<dim> &
 FESystem<dim>::base_element (const unsigned int index) const

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.