]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
mixed discretization with Dirichlet boundary works
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 1999 09:16:20 +0000 (09:16 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 1999 09:16:20 +0000 (09:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@997 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/include/numerics/vectors.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_system.cc
deal.II/deal.II/source/numerics/vectors.cc

index bc8c05745a03291285bbcc8fef845b1ad3c89b5f..459711ffe970852f27b4d72d933689aaa2ae52e2 100644 (file)
@@ -70,6 +70,16 @@ class FiniteElementData
                                      */
     const unsigned int first_hex_index;
     
+                                    /**
+                                     * First index of dof on a line for face data.
+                                     */
+    const unsigned int first_face_line_index;
+    
+                                    /**
+                                     * First index of dof on a quad for face data.
+                                     */
+    const unsigned int first_face_quad_index;
+
                                     /**
                                      * Total number of degrees of freedom
                                      * on a cell. This information is
@@ -239,8 +249,8 @@ class FiniteElementBase :
                                     /**
                                      * Compute system index from components.
                                      */
-    unsigned component_to_system_index (unsigned component,
-                                       unsigned component_index) const;
+    unsigned int component_to_system_index (unsigned int component,
+                                       unsigned int component_index) const;
   
                                     /**
                                      * Compute component and index from system index.
@@ -249,7 +259,23 @@ class FiniteElementBase :
                                      * component and second index in
                                      * component.
                                      */
-    pair<unsigned,unsigned> system_to_component_index (unsigned index) const; 
+    pair<unsigned int,unsigned int> system_to_component_index (unsigned int index) const; 
+    
+                                    /**
+                                     * Compute system index from components on a face.
+                                     */
+    unsigned int face_component_to_system_index (unsigned int component,
+                                                unsigned int component_index) const;
+  
+                                    /**
+                                     * Compute component and index from system
+                                     * index for a face.
+                                     *
+                                     * Return value contains first
+                                     * component and second index in
+                                     * component.
+                                     */
+    pair<unsigned int,unsigned int> face_system_to_component_index (unsigned int index) const; 
     
 
                                     /**
@@ -430,12 +456,22 @@ class FiniteElementBase :
                                     /**
                                      * Map between linear dofs and component dofs.
                                      */
-    vector< pair<unsigned, unsigned> > system_to_component_table;
+    vector< pair<unsigned int, unsigned int> > system_to_component_table;
+
+                                    /**
+                                     * Map between linear dofs and component dofs on face.
+                                     */
+    vector< pair<unsigned int, unsigned int> > face_system_to_component_table;
 
                                     /**
                                      * Map between component and linear dofs.
                                      */
-    vector< vector<unsigned> > component_to_system_table;
+    vector< vector<unsigned int> > component_to_system_table;
+
+                                    /**
+                                     * Map between component and linear dofs on a face.
+                                     */
+    vector< vector<unsigned int> > face_component_to_system_table;
 };
 
 
@@ -1333,7 +1369,7 @@ class FiniteElement : public FiniteElementBase<dim>
                                      * components of mixed
                                      * discretizations.
                                      */
-    virtual const FiniteElement<dim>& base_element(unsigned index) const;
+    virtual const FiniteElement<dim>& base_element(unsigned int index) const;
     
                                     /**
                                      * Exception
@@ -1362,9 +1398,9 @@ class FiniteElement : public FiniteElementBase<dim>
 };
 
 template <int dim>
-inline unsigned
-FiniteElementBase<dim>::component_to_system_index (unsigned component,
-                                                  unsigned component_index) const
+inline unsigned int
+FiniteElementBase<dim>::component_to_system_index (unsigned int component,
+                                                  unsigned int component_index) const
 {
   Assert(component<n_components, ExcInvalidIndex(component));
   Assert(component_index<component_to_system_table[component].size(),
@@ -1373,13 +1409,32 @@ FiniteElementBase<dim>::component_to_system_index (unsigned component,
 }
 
 template <int dim>  
-inline pair<unsigned,unsigned>
-FiniteElementBase<dim>::system_to_component_index (unsigned index) const
+inline pair<unsigned int,unsigned int>
+FiniteElementBase<dim>::system_to_component_index (unsigned int index) const
 {
   Assert(index < system_to_component_table.size(), ExcInvalidIndex(index));
   return system_to_component_table[index];
 }
 
+template <int dim>
+inline unsigned int
+FiniteElementBase<dim>::face_component_to_system_index (unsigned int component,
+                                                       unsigned int component_index) const
+{
+  Assert(component<n_components, ExcInvalidIndex(component));
+  Assert(component_index<face_component_to_system_table[component].size(),
+        ExcInvalidIndex(component_index));
+  return face_component_to_system_table[component][component_index];
+}
+
+template <int dim>  
+inline pair<unsigned int,unsigned int>
+FiniteElementBase<dim>::face_system_to_component_index (unsigned int index) const
+{
+  Assert(index < system_to_component_table.size(), ExcInvalidIndex(index));
+  return face_system_to_component_table[index];
+}
+
 /*----------------------------   fe.h     ---------------------------*/
 /* end of #ifndef __fe_H */
 #endif
index 19e51b3aaad5d38e12a96de5d15cb3d44157021a..8d6d5d385926c2c3419ad5e2df0f0586127520f2 100644 (file)
@@ -445,6 +445,17 @@ class FESystem //<dim>
                                      * the #.h# file.
                                      */
     void initialize();
+
+                                    /**
+                                     * Used by #initialize#.
+                                     */
+    void build_cell_table();
+    
+                                    /**
+                                     * Used by #initialize#.
+                                     */
+    void build_face_table();
+    
                                     /**
                                      *Exception.
                                      */
index 21e09386ba237764a5523929867bc99e51943265..2e1d5e432aa53e5635f29f938a1a09a47f9ce93d 100644 (file)
@@ -243,7 +243,8 @@ enum NormType {
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim>
-class VectorTools {
+class VectorTools //<dim>
+{
   public:
                                     /**
                                      * Declare a data type which denotes a
@@ -259,6 +260,14 @@ class VectorTools {
                                      */
     typedef map<unsigned char,const Function<dim>*> FunctionMap;
 
+                                    /**
+                                     * Data type for vector valued boundary function map.
+                                     */
+    typedef map<unsigned char,const VectorFunction<dim>*> VectorFunctionMap;
+
+//    typedef map<pair<unsigned char, unsigned int> ,const Function<dim>*>
+//    VectorFunctionMap;
+    
                                     /**
                                      * Compute the interpolation of
                                      * #function# at the support points to
@@ -329,10 +338,10 @@ class VectorTools {
                                        Vector<double>           &rhs_vector);
     
                                     /**
-                                     * Make up the list of node subject
+                                     * Make up the list of nodes subject
                                      * to Dirichlet boundary conditions
-                                     * and the values they are to be
-                                     * assigned, by interpolation around
+                                     * and the values to be
+                                     * assigned to them, by interpolation around
                                      * the boundary. If the
                                      * #boundary_values# contained values
                                      * before, the new ones are added, or
@@ -347,7 +356,16 @@ class VectorTools {
                                             const FunctionMap     &dirichlet_bc,
                                             const Boundary<dim>   &boundary,
                                             map<int,double>       &boundary_values);
-    
+
+                                    /**
+                                     * Create boundary value information for vector
+                                     * valued functions.
+                                     * See the other #interpolate_boundary_values#.
+                                     */
+    static void interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                            const VectorFunctionMap     &dirichlet_bc,
+                                            const Boundary<dim>   &boundary,
+                                            map<int,double>       &boundary_values);
                                     /**
                                      * Project #function# to the boundary
                                      * of the domain, using the given quadrature
index f28668aeb96c4643cf58b79dc9330482b4cfe5fd..f4dbb5fddec31de23bc7f264a553168d71a92c50 100644 (file)
@@ -36,11 +36,19 @@ FiniteElementData<dim>::FiniteElementData (const unsigned int dofs_per_vertex,
                dofs_per_face(GeometryInfo<dim>::vertices_per_face * dofs_per_vertex+
                              GeometryInfo<dim>::lines_per_face * dofs_per_line +
                              dofs_per_quad),
-               first_line_index(GeometryInfo<dim>::vertices_per_cell * dofs_per_vertex),
+               first_line_index(GeometryInfo<dim>::vertices_per_cell
+                                * dofs_per_vertex),
                first_quad_index(first_line_index+
-                                12*dofs_per_line),
+                                GeometryInfo<dim>::lines_per_cell
+                                * dofs_per_line),
                first_hex_index(first_quad_index+
-                               GeometryInfo<dim>::faces_per_cell * dofs_per_quad),
+                               GeometryInfo<dim>::faces_per_cell
+                               * dofs_per_quad),
+               first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+                                     * dofs_per_vertex),
+               first_face_quad_index(first_line_index+
+                                     GeometryInfo<dim-1>::lines_per_cell
+                                     * dofs_per_line),
                total_dofs (first_hex_index+dofs_per_hex),
                n_transform_functions (n_transform_functions),
                n_components(n_components)
@@ -67,6 +75,11 @@ FiniteElementData<dim>::FiniteElementData (const unsigned int dofs_per_vertex,
                                 GeometryInfo<dim>::lines_per_cell * dofs_per_line),
                first_hex_index(first_quad_index+
                                GeometryInfo<dim>::quads_per_cell*dofs_per_quad),
+               first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+                                     * dofs_per_vertex),
+               first_face_quad_index(first_line_index+
+                                     GeometryInfo<dim-1>::lines_per_cell
+                                     * dofs_per_line),
                total_dofs (first_quad_index+dofs_per_quad),
                n_transform_functions (n_transform_functions),
                n_components(n_components)
@@ -91,6 +104,11 @@ FiniteElementData<dim>::FiniteElementData (const unsigned int dofs_per_vertex,
                                 GeometryInfo<dim>::lines_per_cell * dofs_per_line),
                first_hex_index(first_quad_index+
                                GeometryInfo<dim>::quads_per_cell*dofs_per_quad),
+               first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+                                     * dofs_per_vertex),
+               first_face_quad_index(first_line_index+
+                                     GeometryInfo<dim-1>::lines_per_cell
+                                     * dofs_per_line),
                total_dofs (first_line_index+dofs_per_line),
                n_transform_functions (n_transform_functions),
                n_components(n_components)
@@ -120,7 +138,9 @@ template <int dim>
 FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data) :
                FiniteElementData<dim> (fe_data),
                system_to_component_table(total_dofs),
-               component_to_system_table(n_components, vector<unsigned>(total_dofs))
+               face_system_to_component_table(dofs_per_face),
+               component_to_system_table(n_components, vector<unsigned>(total_dofs)),
+               face_component_to_system_table(n_components, vector<unsigned>(dofs_per_face))
 {
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i) 
     {
@@ -160,6 +180,11 @@ FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data
       system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
       component_to_system_table[0][j] = j;
     }
+  for (unsigned int j=0 ; j<dofs_per_face ; ++j)
+    {
+      face_system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+      face_component_to_system_table[0][j] = j;
+    }
 };
 
 
index 3bd2efb2042774288e2d6db1df71cf18b0ab11af..c7cbb58a4863f7af2f725a0c57c3b4f21b04e795 100644 (file)
@@ -19,12 +19,10 @@ FESystem<dim>::~FESystem ()
 };
 
 
-
 template <int dim>
-void FESystem<dim>::initialize ()
+void
+FESystem<dim>::build_cell_table()
 {
-                                  // 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)
@@ -35,7 +33,7 @@ void FESystem<dim>::initialize ()
                                   // to be thought of.
   
                                   // 1. Vertices
-      total_index = 0;
+  total_index = 0;
   for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_cell ;
        ++vertex_number)
     {
@@ -137,7 +135,8 @@ void FESystem<dim>::initialize ()
          
        }
     }
-  
+                                  // Inintialize mapping from component
+                                  // to base element
                                   // Initialize mapping from components to
                                   // linear index. Fortunately, this is
                                   // the inverse of what we just did.
@@ -148,8 +147,120 @@ void FESystem<dim>::initialize ()
   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;
+}
+
+
+template <int dim>
+void
+FESystem<dim>::build_face_table()
+{
+  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.
+  
+                                  // 1. Vertices
+  total_index = 0;
+  for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_face ;
+       ++vertex_number)
+    {
+      unsigned comp_start = 0;
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
+       {
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
+           {
+             for (unsigned local_index = 0 ;
+                  local_index < base_element(base).dofs_per_vertex ;
+                  ++local_index)
+               {
+                 face_system_to_component_table[total_index++]
+                   = pair<unsigned,unsigned>
+                   (comp_start+m,
+                    vertex_number*base_element(base).dofs_per_vertex+local_index);
+               }
+           }
+         comp_start += element_multiplicity(base);
+       }
+    }
+  
+                                  // 2. Lines
+  for (unsigned line_number= 0 ; ((line_number != GeometryInfo<dim>::lines_per_face) &&
+                                 (GeometryInfo<dim>::lines_per_cell > 0));
+       ++line_number)
+    {
+      unsigned comp_start = 0;
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
+       {
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
+           {
+             for (unsigned local_index = 0 ;
+                  local_index < base_element(base).dofs_per_line ;
+                  ++local_index)
+               {
+                 face_system_to_component_table[total_index++]
+                   = pair<unsigned,unsigned>
+                   (comp_start+m,
+                    line_number*base_element(base).dofs_per_line
+                    +local_index+base_element(base).first_face_line_index);
+               }
+           }
+         comp_start += element_multiplicity(base);
+       }
+    }
+  
+                                  // 3. Quads
+  for (unsigned quad_number= 0 ; ((quad_number != GeometryInfo<dim>::quads_per_face) &&
+                                 (GeometryInfo<dim>::quads_per_cell > 0));
+       ++quad_number)
+    {
+      unsigned comp_start = 0;
+      for(unsigned base = 0; base < n_base_elements() ;
+         ++base)
+       {
+         for (unsigned m = 0; m < element_multiplicity(base); ++m)
+           {
+             for (unsigned local_index = 0 ;
+                  local_index < base_element(base).dofs_per_quad ;
+                  ++local_index)
+               {
+                 face_system_to_component_table[total_index++]
+                   = pair<unsigned,unsigned>
+                   (comp_start+m,
+                    quad_number*base_element(base).dofs_per_quad
+                    +local_index+base_element(base).first_face_quad_index);
+               }
+           }
+         comp_start += element_multiplicity(base);
+       }
+    }
   
+                                  // Inintialize mapping from component
+                                  // to base element
+                                  // 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)
+    face_component_to_system_table[comp]
+      .resize(base_element(component_to_base_table[comp]).total_dofs);
 
+  for (unsigned sys=0 ; sys < dofs_per_face ; ++sys)
+    face_component_to_system_table[face_system_to_component_table[sys].first]
+      [face_system_to_component_table[sys].second] = sys;
+}
+
+
+
+template <int dim>
+void FESystem<dim>::initialize ()
+{
+  build_cell_table();
+  build_face_table();
                                   // distribute the matrices of the base
                                   // finite element to the matrices of
                                   // this object
@@ -324,23 +435,31 @@ void FESystem<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &/*
 
 
 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());
-/*
   Assert (support_points.size() == dofs_per_face,
          ExcWrongFieldDimension (support_points.size(), dofs_per_face));
 
-  vector<Point<dim> > base_support_points (base_element->dofs_per_face);
-  base_element->get_face_support_points (face, boundary, base_support_points);
-  
-  for (unsigned int i=0; i<base_element->dofs_per_face; ++i) 
-    for (unsigned int n=0; n<n_sub_elements; ++n)
-      support_points[i*n_sub_elements+n] = base_support_points[i];
-*/
-};
+  vector<Point<dim> > base_support_points (base_element(0).dofs_per_face);
+  unsigned int comp = 0;
+  for (unsigned int base=0 ; base<n_base_elements(); ++base)
+    {
+      base_support_points.resize(base_element(base).dofs_per_face);
+      base_element(base).get_face_support_points (face, boundary, base_support_points);
+      for (unsigned int inbase = 0 ; inbase < element_multiplicity(base); ++inbase)
+       {
+         for (unsigned int i=0; i<base_element(base).dofs_per_face; ++i)
+           {
+             support_points[face_component_to_system_index(comp,i)]
+               = base_support_points[i];
+           }
+         
+             ++comp;
+       }
+    }
+}
 
 
 
@@ -358,15 +477,15 @@ void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator
 
                                   // first get the local mass matrix for
                                   // the base object
-  FullMatrix<double> base_mass_matrix (base_element->total_dofs,
-                            base_element->total_dofs);
-  base_element->get_local_mass_matrix (cell, boundary, base_mass_matrix);
+  FullMatrix<double> base_mass_matrix (base_element.total_dofs,
+                            base_element.total_dofs);
+  base_element.get_local_mass_matrix (cell, boundary, base_mass_matrix);
 
 
                                   // now distribute it to the mass matrix
                                   // of this object
-  for (unsigned int i=0; i<base_element->total_dofs; ++i)
-    for (unsigned int j=0; j<base_element->total_dofs; ++j)
+  for (unsigned int i=0; i<base_element.total_dofs; ++i)
+    for (unsigned int j=0; j<base_element.total_dofs; ++j)
       for (unsigned int n=0; n<n_sub_elements; ++n)
                                         // only fill diagonals of the blocks
        local_mass_matrix (i*n_sub_elements + n,
index 68785d54dd4060b2c1c1d3f4c23c789d696bfa28..426f4af048028ef94d6b6b5e8ec83b9f936b6f6e 100644 (file)
@@ -241,7 +241,17 @@ void
 VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
                                             const FunctionMap &,
                                             const Boundary<1> &,
-                                            map<int,double>   &) {
+                                            map<int,double>   &)
+{
+  Assert (false, ExcNotImplemented());
+};
+
+template <>
+void VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
+                                                 const VectorFunctionMap&,
+                                                 const Boundary<1>&,
+                                                 map<int,double>&)
+{
   Assert (false, ExcNotImplemented());
 };
 
@@ -292,6 +302,56 @@ VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
       };
 };
 
+template <int dim>
+void
+VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                              const VectorFunctionMap     &dirichlet_bc,
+                                              const Boundary<dim>      &boundary,
+                                              map<int,double>   &boundary_values)
+{
+  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
+         ExcInvalidBoundaryIndicator());
+
+  const FiniteElement<dim> &fe = dof.get_fe();
+
+                                  // use two face iterators, since we need
+                                  // a DoF-iterator for the dof indices, but
+                                  // a Tria-iterator for the fe object
+  DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+                                       endf = dof.end_face();
+  
+  typename VectorFunctionMap::const_iterator function_ptr;
+
+                                  // field to store the indices of dofs
+  vector<int>         face_dofs (fe.dofs_per_face);
+  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+  vector< Vector<double> > dof_values (fe.dofs_per_face, Vector<double>(fe.n_components));
+       
+  for (; face!=endf; ++face)
+    if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
+       dirichlet_bc.end()) 
+                                      // face is subject to one of the
+                                      // bc listed in #dirichlet_bc#
+      {
+                                        // get indices, physical location and
+                                        // boundary values of dofs on this
+                                        // face
+       face->get_dof_indices (face_dofs);
+       fe.get_face_support_points (face, boundary, dof_locations);
+       function_ptr->second->value_list (dof_locations, dof_values);
+
+                                        // enter into list
+
+       for (unsigned int i=0; i<face_dofs.size(); ++i)
+         {
+           pair<unsigned int, unsigned int>
+             index = fe.face_system_to_component_index(i);
+                                                   
+           boundary_values[face_dofs[i]] = dof_values[i](index.first);
+         }
+      }
+}
+
 
 
 template <int dim>

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.