]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Most functions in FESystem implemented, the rest may be shouldnt
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Feb 1999 13:40:14 +0000 (13:40 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Feb 1999 13:40:14 +0000 (13:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@813 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.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_system.cc
deal.II/lac/include/lac/dsmatrix.h

index a7ac0c29303e65b06e363a0020140660f8a8d4cd..22ae8255644ca888ead029418793c645d7fdaac6 100644 (file)
@@ -435,6 +435,11 @@ class FiniteElementBase :
                                      * Map between linear dofs and component dofs.
                                      */
     vector< pair<unsigned, unsigned> > system_to_component_table;
+
+                                    /**
+                                     * Map between component and linear dofs.
+                                     */
+    vector< vector<unsigned> > component_to_system_table;
 };
 
 
@@ -1130,11 +1135,13 @@ class FiniteElement : public FiniteElementBase<dim> {
 
 template <int dim>
 inline unsigned
-FiniteElementBase<dim>::component_to_system_index (unsigned /*component*/,
+FiniteElementBase<dim>::component_to_system_index (unsigned component,
                                                   unsigned component_index) const
 {
-  Assert(false, ExcNotImplemented());
-  return component_index;
+  Assert(component<n_components, ExcInvalidIndex(component));
+  Assert(component_index<component_to_system_table[component].size(),
+        ExcInvalidIndex(component_index));
+  return component_to_system_table[component][component_index];
 }
 
 template <int dim>  
index 5a5beaabf80bbbd48a681e164bf9f32d657f8e6d..d7026356834e3b7d0cdb83717d646abb2e8f1188 100644 (file)
@@ -101,7 +101,7 @@ class FESystem : public FiniteElement<dim>
                                     /** 
                                      * Constructor for mixed
                                      * discretizations with two
-                                     * components.
+                                     * base elements.
                                      *
                                      * See the other constructor.
                                      */
@@ -317,7 +317,7 @@ class FESystem : public FiniteElement<dim>
                                     vector<Point<dim> >         &normal_vectors) const;
 
                                     /** 
-                                     *Number of different composing
+                                     * Number of different base
                                      * elements of this object.
                                      *
                                      * Since these objects can have
@@ -327,7 +327,7 @@ class FESystem : public FiniteElement<dim>
                                      * of finite elements composed
                                      * into this structure.
                                      */
-    unsigned n_component_elements() const;
+    unsigned n_base_elements() const;
 
                                     /**
                                      * How often is a composing element used.
@@ -372,6 +372,21 @@ class FESystem : public FiniteElement<dim>
                                      */
     vector< ElementPair > base_elements;
 
+                                    /**
+                                     * The base element establishing a
+                                     * component.
+                                     *
+                                     * This table converts a
+                                     * component number to the
+                                     * #base_element# number. While
+                                     * component information contains
+                                     * multiplicity of base elements,
+                                     * the result allows access to
+                                     * shape functions of the base
+                                     * element.
+                                     *
+                                     */
+    vector<unsigned> component_to_base_table;
     
                                     /**
                                      * Helper function used in the constructor:
@@ -382,8 +397,7 @@ class FESystem : public FiniteElement<dim>
                                      * multiplied by #n#. Don't touch the
                                      * number of functions for the
                                      * transformation from unit to real
-                                     * cell.
-                                     */
+                                     * cell.  */
     static FiniteElementData<dim>
     multiply_dof_numbers (const FiniteElementData<dim> &fe_data,
                          const unsigned int            N);
@@ -421,7 +435,7 @@ class FESystem : public FiniteElement<dim>
 
 template<int dim>
 inline unsigned
-FESystem<dim>::n_component_elements() const
+FESystem<dim>::n_base_elements() const
 {
   return base_elements.size();
 }
@@ -447,7 +461,8 @@ template <int dim>
 template <class FE>
 FESystem<dim>::FESystem (const FE &fe, const unsigned int n_elements) :
                FiniteElement (multiply_dof_numbers(fe, n_elements)),
-               base_elements(1)
+               base_elements(1),
+               component_to_base_table(n_components)
 {
   base_elements[0] = ElementPair(new FE, n_elements);
   base_elements[0].first -> subscribe ();
@@ -460,7 +475,8 @@ FESystem<dim>::FESystem (const FE1 &fe1, const unsigned int n1,
                         const FE2 &fe2, const unsigned int n2)
                :
                FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2)),
-               base_elements(2)
+               base_elements(2),
+               component_to_base_table(n_components)
 {
   base_elements[0] = ElementPair(new FE1, n1);
   base_elements[0].first -> subscribe ();
index f8e488b8fd54b3dc70f6d45b7138ac9c7ff2894d..ce0a9290e51ac36413b77dd9dc940ace6d3c2605 100644 (file)
@@ -125,7 +125,8 @@ bool FiniteElementData<dim>::operator== (const FiniteElementData<2> &f) const
 template <>
 FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) :
                FiniteElementData<1> (fe_data),
-               system_to_component_table(total_dofs)
+               system_to_component_table(total_dofs),
+               component_to_system_table(n_components, vector<unsigned>(total_dofs))
 {
   const unsigned int dim=1;
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i) 
@@ -137,7 +138,10 @@ FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) :
   interface_constraints(0,0)=1.;
 
   for (unsigned int j=0 ; j<total_dofs ; ++j)
-    system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+    {
+      system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+      component_to_system_table[0][j] = j;
+    }
 };
 
 #endif
@@ -148,7 +152,8 @@ FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) :
 template <>
 FiniteElementBase<2>::FiniteElementBase (const FiniteElementData<2> &fe_data) :
                FiniteElementData<2> (fe_data),
-               system_to_component_table(total_dofs)
+               system_to_component_table(total_dofs),
+               component_to_system_table(n_components, vector<unsigned>(total_dofs))
 {
   const unsigned int dim=2;
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i) 
@@ -159,7 +164,10 @@ FiniteElementBase<2>::FiniteElementBase (const FiniteElementData<2> &fe_data) :
   interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line,
                                2*dofs_per_vertex+dofs_per_line);
   for (unsigned int j=0 ; j<total_dofs ; ++j)
-    system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+    {
+      system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+      component_to_system_table[0][j] = j;
+    }
 };
 
 #endif
index fc4d87be23638f2e0a932dffbbee19f8f953e6f2..94991e57d4a6e1aae54c73576c7febc24625f68b 100644 (file)
@@ -23,30 +23,39 @@ FESystem<dim>::~FESystem ()
 template <int dim>
 void FESystem<dim>::initialize ()
 {
+                                  // Inintialize mapping from component
+                                  // to base element
+  unsigned total_index = 0;
+  for (unsigned base=0 ; base < n_base_elements() ; ++base)
+    for (unsigned m = 0; m < element_multiplicity(base); ++m)
+      component_to_base_table[total_index++] = base;
+
                                   // Initialize index table
-                                  // Multi-component base elements have to be thought of.
+                                  // Multi-component base elements have
+                                  // to be thought of.
+  
                                   // 1. Vertices
-  unsigned total_index = 0;
+      total_index = 0;
   for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<2>::vertices_per_cell ;
        ++vertex_number)
     {
       unsigned comp_start = 0;
-      for(unsigned comp = 0; comp < n_component_elements() ;
-         ++comp)
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
        {
-         for (unsigned m = 0; m < element_multiplicity(comp); ++m)
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
            {
              for (unsigned local_index = 0 ;
-                  local_index < base_element(comp).dofs_per_vertex ;
+                  local_index < base_element(base).dofs_per_vertex ;
                   ++local_index)
                {
                  system_to_component_table[total_index++]
                    = pair<unsigned,unsigned>
                    (comp_start+m,
-                    vertex_number*base_element(comp).dofs_per_vertex+local_index);
+                    vertex_number*base_element(base).dofs_per_vertex+local_index);
                }
            }
-         comp_start += element_multiplicity(comp);
+         comp_start += element_multiplicity(base);
        }
     }
   
@@ -55,23 +64,23 @@ void FESystem<dim>::initialize ()
        ++line_number)
     {
       unsigned comp_start = 0;
-      for(unsigned comp = 0; comp < n_component_elements() ;
-         ++comp)
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
        {
-         for (unsigned m = 0; m < element_multiplicity(comp); ++m)
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
            {
              for (unsigned local_index = 0 ;
-                  local_index < base_element(comp).dofs_per_line ;
+                  local_index < base_element(base).dofs_per_line ;
                   ++local_index)
                {
                  system_to_component_table[total_index++]
                    = pair<unsigned,unsigned>
                    (comp_start+m,
-                    line_number*base_element(comp).dofs_per_line
-                    +local_index+base_element(comp).first_line_index);
+                    line_number*base_element(base).dofs_per_line
+                    +local_index+base_element(base).first_line_index);
                }
            }
-         comp_start += element_multiplicity(comp);
+         comp_start += element_multiplicity(base);
        }
     }
   
@@ -80,23 +89,23 @@ void FESystem<dim>::initialize ()
        ++quad_number)
     {
       unsigned comp_start = 0;
-      for(unsigned comp = 0; comp < n_component_elements() ;
-         ++comp)
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
        {
-         for (unsigned m = 0; m < element_multiplicity(comp); ++m)
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
            {
              for (unsigned local_index = 0 ;
-                  local_index < base_element(comp).dofs_per_quad ;
+                  local_index < base_element(base).dofs_per_quad ;
                   ++local_index)
                {
                  system_to_component_table[total_index++]
                    = pair<unsigned,unsigned>
                    (comp_start+m,
-                    quad_number*base_element(comp).dofs_per_quad
-                    +local_index+base_element(comp).first_quad_index);
+                    quad_number*base_element(base).dofs_per_quad
+                    +local_index+base_element(base).first_quad_index);
                }
            }
-         comp_start += element_multiplicity(comp);
+         comp_start += element_multiplicity(base);
        }
     }
   
@@ -105,25 +114,37 @@ void FESystem<dim>::initialize ()
        ++hex_number)
     {
       unsigned comp_start = 0;
-      for(unsigned comp = 0; comp < n_component_elements() ;
-         ++comp, comp_start += element_multiplicity(comp))
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
        {
-         for (unsigned m = 0; m < element_multiplicity(comp); ++m)
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
            {
              for (unsigned local_index = 0 ;
-                  local_index < base_element(comp).dofs_per_hex ;
+                  local_index < base_element(base).dofs_per_hex ;
                   ++local_index)
                {
                  system_to_component_table[total_index++]
                    = pair<unsigned,unsigned>
                    (comp_start+m,
-                    hex_number*base_element(comp).dofs_per_hex
-                    +local_index+base_element(comp).first_hex_index);
+                    hex_number*base_element(base).dofs_per_hex
+                    +local_index+base_element(base).first_hex_index);
                }
            }
+         comp_start += element_multiplicity(base);
+         
        }
     }
   
+                                  // Initialize mapping from components to
+                                  // linear index. Fortunately, this is
+                                  // the inverse of what we just did.
+  for (unsigned comp=0 ; comp<n_components ; ++comp)
+    component_to_system_table[comp]
+      .resize(base_element(component_to_base_table[comp]).total_dofs);
+
+  for (unsigned sys=0 ; sys < total_dofs ; ++sys)
+    component_to_system_table[system_to_component_table[sys].first]
+      [system_to_component_table[sys].second] = sys;
   
 
                                   // distribute the matrices of the base
@@ -214,16 +235,16 @@ FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe1,
 
 
 template <int dim>
-double FESystem<dim>::shape_value (const unsigned int i,
-                                  const Point<dim>  &p) const
+double
+FESystem<dim>::shape_value (const unsigned int i,
+                           const Point<dim>  &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return 0.;
-  
   Assert((i<total_dofs), ExcInvalidIndex(i));
 
-
-//  return base_element->shape_value (i / n_sub_elements, p);
+  pair<unsigned,unsigned> comp = system_to_component_index(i);
+  
+  return base_element(component_to_base_table[comp.first])
+    .shape_value(comp.second, p);
 };
 
 
@@ -233,12 +254,12 @@ Tensor<1,dim>
 FESystem<dim>::shape_grad (const unsigned int  i,
                           const Point<dim>   &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return Tensor<1,dim>();
-
   Assert((i<total_dofs), ExcInvalidIndex(i));
 
-//  return base_element->shape_grad (i / n_sub_elements, p);
+  pair<unsigned,unsigned> comp = system_to_component_index(i);
+  
+  return base_element(component_to_base_table[comp.first])
+    .shape_grad(comp.second, p);
 };
 
 
@@ -248,20 +269,23 @@ Tensor<2,dim>
 FESystem<dim>::shape_grad_grad (const unsigned int  i,
                                const Point<dim>   &p) const
  {
-  Assert(false, ExcNotImplemented());
-  return Tensor<2,dim>();
-
   Assert((i<total_dofs), ExcInvalidIndex(i));
 
-//  return base_element->shape_grad_grad (i / n_sub_elements, p);
+
+  pair<unsigned,unsigned> comp = system_to_component_index(i);
+  
+  return base_element(component_to_base_table[comp.first])
+    .shape_grad_grad(comp.second, p);
 };
 
 
 
 template <int dim>
-void FESystem<dim>::get_unit_support_points (vector<Point<dim> > &support_points) const
+void FESystem<dim>::get_unit_support_points (vector<Point<dim> > &/*support_points*/) const
 {
-/*  Assert (support_points.size() == total_dofs,
+  Assert(false, ExcNotImplemented());
+/*
+  Assert (support_points.size() == total_dofs,
          ExcWrongFieldDimension (support_points.size(), total_dofs));
 
   vector<Point<dim> > base_support_points (base_element->total_dofs);
@@ -276,9 +300,9 @@ void FESystem<dim>::get_unit_support_points (vector<Point<dim> > &support_points
 
 
 template <int dim>
-void FESystem<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &cell,
-                                       const Boundary<dim> &boundary,
-                                       vector<Point<dim> > &support_points) const
+void FESystem<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &/*cell*/,
+                                       const Boundary<dim> &/*boundary*/,
+                                       vector<Point<dim> > &/*support_points*/) const
 {
   Assert(false, ExcNotImplemented());
 /*
@@ -297,9 +321,9 @@ void FESystem<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &ce
 
 
 template <int dim>
-void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterator &face,
-                                            const Boundary<dim> &boundary,
-                                            vector<Point<dim> > &support_points) const
+void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterator &/*face*/,
+                                            const Boundary<dim> &/*boundary*/,
+                                            vector<Point<dim> > &/*support_points*/) const
 {
   Assert(false, ExcNotImplemented());
 /*
@@ -318,9 +342,9 @@ void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterato
 
 
 template <int dim>
-void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
-                                          const Boundary<dim> &boundary,
-                                          dFMatrix            &local_mass_matrix) const
+void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &/*cell*/,
+                                          const Boundary<dim> &/*boundary*/,
+                                          dFMatrix            &/*local_mass_matrix*/) const
 {
   Assert(false, ExcNotImplemented());
 /*
@@ -350,19 +374,17 @@ void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator
 
 
 template <int dim>
-double FESystem<dim>::shape_value_transform (const unsigned int i,
-                                            const Point<dim>  &p) const
+double FESystem<dim>::shape_value_transform (const unsigned int /*i*/,
+                                            const Point<dim>  &/*p*/) const
 {
   Assert(false, ExcNotImplemented());
   return 0.;
-  
-//  return base_element->shape_value_transform (i, p);
 };
 
 
 template <int dim>
-Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int i,
-                                                  const Point<dim>  &p) const
+Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int /*i*/,
+                                                  const Point<dim>  &/*p*/) const
 {
   Assert(false, ExcNotImplemented());
   return Tensor<1,dim>();
@@ -373,10 +395,10 @@ Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int i,
 
 
 template <int dim>
-void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                       const Boundary<dim>         &boundary,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<double>      &face_jacobi_determinants) const
+void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
+                                       const Boundary<dim>         &/*boundary*/,
+                                       const vector<Point<dim-1> > &/*unit_points*/,
+                                       vector<double>      &/*face_jacobi_determinants*/) const
 {
   Assert(false, ExcNotImplemented());
 
@@ -386,10 +408,10 @@ void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &fa
 
 
 template <int dim>
-void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                          const unsigned int           subface_no,
-                                          const vector<Point<dim-1> > &unit_points,
-                                          vector<double>      &face_jacobi_determinants) const
+void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
+                                          const unsigned int           /*subface_no*/,
+                                          const vector<Point<dim-1> > &/*unit_points*/,
+                                          vector<double>      &/*face_jacobi_determinants*/) const
 {
   Assert(false, ExcNotImplemented());
 
@@ -400,11 +422,11 @@ void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator
 
 
 template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                       const unsigned int          face_no,
-                                       const Boundary<dim>         &boundary,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<Point<dim> >         &normal_vectors) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
+                                       const unsigned int          /*face_no*/,
+                                       const Boundary<dim>         &/*boundary*/,
+                                       const vector<Point<dim-1> > &/*unit_points*/,
+                                       vector<Point<dim> >         &/*normal_vectors*/) const
 {
   Assert(false, ExcNotImplemented());
 
@@ -414,11 +436,11 @@ void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &ce
 
 
 template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                       const unsigned int          face_no,
-                                       const unsigned int          subface_no,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<Point<dim> >         &normal_vectors) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
+                                       const unsigned int          /*face_no*/,
+                                       const unsigned int          /*subface_no*/,
+                                       const vector<Point<dim-1> > &/*unit_points*/,
+                                       vector<Point<dim> >         &/*normal_vectors*/) const
 {
   Assert(false, ExcNotImplemented());
 
index e08e21ccd4c8b4aab9968d215399ad6ef05647cd..30484ce44a2d49b65a7e4cbbed775c5096e2fcbf 100644 (file)
@@ -475,7 +475,8 @@ class dSMatrix
                                     /**
                                      * Set the element #(i,j)# to #value#.
                                      * Throws an error if the entry does
-                                     * not exist.
+                                     * not exist. Still, it is allowed to store
+                                     * zero values in non-existent fields.
                                      */
     void set (const unsigned int i, const unsigned int j,
              const double value);
@@ -483,7 +484,8 @@ class dSMatrix
                                     /**
                                      * Add #value# to the element #(i,j)#.
                                      * Throws an error if the entry does
-                                     * not exist.
+                                     * not exist. Still, it is allowed to store
+                                     * zero values in non-existent fields.
                                      */
     void add (const unsigned int i, const unsigned int j,
              const double value);
@@ -784,9 +786,12 @@ inline
 void dSMatrix::set (const unsigned int i, const unsigned int j,
                    const double value) {
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->operator()(i,j) != -1,
+  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
          ExcInvalidIndex(i,j));
-  val[cols->operator()(i,j)] = value;
+
+  const int index = cols->operator()(i,j);
+
+  if (index >= 0) val[index] = value;
 };
 
 
@@ -795,9 +800,12 @@ inline
 void dSMatrix::add (const unsigned int i, const unsigned int j,
                    const double value) {
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->operator()(i,j) != -1,
+  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
          ExcInvalidIndex(i,j));
-  val[cols->operator()(i,j)] += value;
+
+  const int index = cols->operator()(i,j);
+  
+  if (index >= 0) val[index] += value;
 };
 
 

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.