]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support for faces of cells and starting support for Dirichlet boundary conditions...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 27 Mar 1998 15:32:42 +0000 (15:32 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 27 Mar 1998 15:32:42 +0000 (15:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@99 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
deal.II/base/include/base/function.h
deal.II/deal.II/Attic/examples/poisson/poisson.cc
deal.II/deal.II/include/dofs/dof_accessor.h
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/dofs/dof_accessor.cc
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/numerics/base.cc
tests/big-tests/poisson/poisson.cc

index aac5e2233612b5c5467d45544588f37f2a8a585a..5ecd440b722c8d6376a6ab735c552e33743a09f4 100644 (file)
@@ -5,19 +5,75 @@
 /*----------------------------   function.h     ---------------------------*/
 
 
-//forward declaration
-template <int dim> class Point;
-template <class T, class Alloc=alloc> class vector;
+#include <base/exceptions.h>
+#include <vector.h>
+#include <grid/point.h>
 
 
 
 
+
+/**
+   This class is a model for a continuous function. It returns the value
+   of the function at a given point through the #operator ()# member function,
+   which is virtual. It also has a function to return a whole list of function
+   values at different points to reduce the overhead of the virtual function
+   calls; this function is preset to successively call the function returning
+   one value at a time.
+
+   There are other functions return the gradient of the function at one or
+   several points. You only have to overload those functions you need; the
+   functions returning several values at a time will call those returning
+   only one value, while those ones will throw an exception when called but
+   not overloaded.
+
+   Unless only called a very small number of times, you should overload
+   both those functions returning only one value as well as those returning
+   a whole array, since the cost of evaluation of a point value is often
+   less than the virtual function call itself.
+*/
 template <int dim>
 class Function {
   public:
-    virtual double operator () (const Point<dim> &p) const = 0;
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const;
+
+                                    /**
+                                     * Set #values# to the point values
+                                     * of the function at the #points#.
+                                     * It is assumed that #values# be
+                                     * empty.
+                                     */
     virtual void value_list (const vector<Point<dim> > &points,
                             vector<double>            &values) const;
+
+                                    /**
+                                     * Return the gradient of the function
+                                     * at the given point.
+                                     */
+    virtual Point<dim> gradient (const Point<dim> &p) const;
+
+                                    /**
+                                     * Set #gradients# to the gradients of
+                                     * the function at the #points#.
+                                     * It is assumed that #values# be
+                                     * empty.
+                                     */
+    virtual void gradient_list (const vector<Point<dim> > &points,
+                               vector<Point<dim> >       &gradients) const;
+
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcPureFunctionCalled);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcVectorNotEmpty);
 };
 
 
index d7761bbc5995dcff6c812c18027240cf340e05a7..a1f644440822e5915d8d0670e5f213ecb7061aef 100644 (file)
@@ -9,12 +9,15 @@
 #include <grid/tria_boundary.h>
 #include <grid/dof_constraints.h>
 #include <basic/data_io.h>
+#include <basic/function.h>
 #include <fe/fe_lib.h>
 #include <fe/quadrature_lib.h>
 #include <numerics/base.h>
 #include <numerics/assembler.h>
 #include <lac/dsmatrix.h>
 
+
+#include <map.h>
 #include <fstream.h>
 #include <cmath>
 extern "C" {
index 63df424b3991e30599c92054244124a60bb17713..925090e51b853c713f981b3449806c58336abff2 100644 (file)
@@ -55,6 +55,10 @@ class DoFAccessor {
                    << "Invalid index " << arg1
                    << ", index must be between " << arg2
                    << " and " << arg3 << ".");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcVectorNotEmpty);
 
   protected:
                                     /**
@@ -164,6 +168,18 @@ class DoFLineAccessor :  public BaseClass, public DoFAccessor<dim> {
                               const unsigned int i,
                               const int          index) const;
 
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * line in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, 
+                                     * dofs on line.
+                                     *
+                                     * The dof indices are appended to the end
+                                     * of the vector. It is assumed that
+                                     * the vector be empty beforehand.
+                                     */
+    void get_dof_indices (vector<int> &dof_indices) const;
+
                                     /**
                                      * Return the #i#th child as a DoF line
                                      * iterator. This function is needed since
@@ -237,6 +253,19 @@ class DoFQuadAccessor :  public BaseClass, public DoFAccessor<dim> {
                               const unsigned int i,
                               const int          index) const;
 
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * quad in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * The dof indices are appended to the end
+                                     * of the vector. It is assumed that
+                                     * the vector be empty beforehand.
+                                     */
+    void get_dof_indices (vector<int> &dof_indices) const;
+
                                     /**
                                      *  Return a pointer to the #i#th line
                                      *  bounding this #Quad#.
@@ -308,11 +337,11 @@ class DoFSubstructAccessor<1> :  public DoFLineAccessor<1,CellAccessor<1> > {
                    DoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
-                                    // what a substruct_iterator is. Otherwise
+                                    // what a face_iterator is. Otherwise
                                     // it would have to ask the DoFHandler<dim>
                                     // which would need another big include
                                     // file and maybe cyclic interdependence
-    typedef void * substruct_iterator;
+    typedef void * face_iterator;
 };
 
 
@@ -333,11 +362,11 @@ class DoFSubstructAccessor<2> : public DoFQuadAccessor<2,CellAccessor<2> > {
                    DoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
-                                    // what a substruct_iterator is. Otherwise
+                                    // what a face_iterator is. Otherwise
                                     // it would have to ask the DoFHandler<dim>
                                     // which would need another big include
                                     // file and maybe cyclic interdependence
-    typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > substruct_iterator;
+    typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > face_iterator;
 };
 
 
@@ -393,22 +422,9 @@ class DoFCellAccessor :  public DoFSubstructAccessor<dim> {
                                      * This function is not implemented in 1D,
                                      * and maps to DoFQuadAccessor::line in 2D.
                                      */
-    typename DoFSubstructAccessor<dim>::substruct_iterator
+    typename DoFSubstructAccessor<dim>::face_iterator
     face (const unsigned int i) const;
 
-                                    /**
-                                     * Return the indices of the dofs of this
-                                     * cell in the standard ordering: dofs
-                                     * on vertex 0, dofs on vertex 1, etc,
-                                     * dofs on line 0, dofs on line 1, etc,
-                                     * dofs on quad 0, etc.
-                                     *
-                                     * The dof indices are appended to the end
-                                     * of the vector. It is assumed that
-                                     * the vector be empty beforehand.
-                                     */
-    void get_dof_indices (vector<int> &dof_indices) const;
-
                                     /**
                                      * Return the value of the given vector
                                      * restricted to the dofs of this
@@ -424,10 +440,6 @@ class DoFCellAccessor :  public DoFSubstructAccessor<dim> {
     void get_dof_values (const dVector  &values,
                         vector<double> &dof_values) const;
 
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcVectorNotEmpty);
                                     /**
                                      * Exception
                                      */
index 360b6f434c83e0abc9262b719b3fe7cbb5b49913..20aebc1c9f93d2b1d9ecf589a6364a4584c435d4 100644 (file)
@@ -153,9 +153,9 @@ class DoFDimensionInfo<1> {
     typedef line_iterator        cell_iterator;
     typedef active_line_iterator active_cell_iterator;
 
-    typedef void * raw_substruct_iterator;
-    typedef void * substruct_iterator;
-    typedef void * active_substruct_iterator;
+    typedef void * raw_face_iterator;
+    typedef void * face_iterator;
+    typedef void * active_face_iterator;
 };
 
 
@@ -182,9 +182,9 @@ class DoFDimensionInfo<2> {
     typedef quad_iterator        cell_iterator;
     typedef active_quad_iterator active_cell_iterator;
 
-    typedef raw_line_iterator    raw_substruct_iterator;
-    typedef line_iterator        substruct_iterator;
-    typedef active_line_iterator active_substruct_iterator;    
+    typedef raw_line_iterator    raw_face_iterator;
+    typedef line_iterator        face_iterator;
+    typedef active_line_iterator active_face_iterator;    
 };
 
 
@@ -436,9 +436,9 @@ class DoFHandler : public DoFDimensionInfo<dim> {
     typedef typename DoFDimensionInfo<dim>::cell_iterator cell_iterator;
     typedef typename DoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
 
-    typedef typename DoFDimensionInfo<dim>::raw_substruct_iterator raw_substruct_iterator;
-    typedef typename DoFDimensionInfo<dim>::substruct_iterator substruct_iterator;
-    typedef typename DoFDimensionInfo<dim>::active_substruct_iterator active_substruct_iterator;
+    typedef typename DoFDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+    typedef typename DoFDimensionInfo<dim>::face_iterator face_iterator;
+    typedef typename DoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
 
     
                                     /**
@@ -703,6 +703,108 @@ class DoFHandler : public DoFDimensionInfo<dim> {
 
                                     /*---------------------------------------*/
 
+                                    /**
+                                     *  @name Face iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first face, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no faces, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 2D and #begin_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    begin_raw_face   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used face
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 2D and #begin_quad# in 3D.
+                                     */
+    face_iterator        begin_face       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  face on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 2D and #begin_active_quad# in 3D.
+                                     */
+    active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 2D and #end_quad# in 3D.
+                                     */
+    raw_face_iterator    end_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last face, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  face of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face (const unsigned int level) const;
+                                    //@}
+
+    
+                                    /*---------------------------------------*/
+
                                     /**
                                      *  @name Line iterator functions
                                      */
@@ -928,6 +1030,10 @@ class DoFHandler : public DoFDimensionInfo<dim> {
     DeclException1 (ExcMatrixHasWrongSize,
                    int,
                    << "The matrix has the wrong dimension " << arg1);
+                                    /**
+                                     *  Exception
+                                     */
+    DeclException0 (ExcFunctionNotUseful);
     
   protected:
                                     /**
index 2685ec001e8240206c1df28f4aee1f7ff359f605..97e3e131038206025e7175e43e5ff24cb01b7820 100644 (file)
@@ -463,7 +463,7 @@ class TriaDimensionInfo;
     #void *#. Thus these types exist, but are useless and will
     certainly make any involuntary use visible.
 
-    The same applies for the #substruct_iterator# types, since lines
+    The same applies for the #face_iterator# types, since lines
     have no substructures apart from vertices, which are handled in
     a different way, however.
     */
@@ -481,9 +481,9 @@ class TriaDimensionInfo<1> {
     typedef line_iterator        cell_iterator;
     typedef active_line_iterator active_cell_iterator;
 
-    typedef void * raw_substruct_iterator;
-    typedef void * substruct_iterator;
-    typedef void * active_substruct_iterator;
+    typedef void * raw_face_iterator;
+    typedef void * face_iterator;
+    typedef void * active_face_iterator;
 };
 
 
@@ -508,9 +508,9 @@ class TriaDimensionInfo<1> {
       typedef quad_iterator        cell_iterator;
       typedef active_quad_iterator active_cell_iterator;
 
-      typedef raw_line_iterator    raw_substruct_iterator;
-      typedef line_iterator        substruct_iterator;
-      typedef active_line_iterator active_substruct_iterator;    
+      typedef raw_line_iterator    raw_face_iterator;
+      typedef line_iterator        face_iterator;
+      typedef active_line_iterator active_face_iterator;    
     \end{verbatim}
     */
 class TriaDimensionInfo<2> {
@@ -527,9 +527,9 @@ class TriaDimensionInfo<2> {
     typedef quad_iterator        cell_iterator;
     typedef active_quad_iterator active_cell_iterator;
 
-    typedef raw_line_iterator    raw_substruct_iterator;
-    typedef line_iterator        substruct_iterator;
-    typedef active_line_iterator active_substruct_iterator;    
+    typedef raw_line_iterator    raw_face_iterator;
+    typedef line_iterator        face_iterator;
+    typedef active_line_iterator active_face_iterator;    
 };
 
 
@@ -601,9 +601,9 @@ class TriaDimensionInfo<2> {
       typedef quad_iterator        cell_iterator;
       typedef active_quad_iterator active_cell_iterator;
 
-      typedef raw_line_iterator    raw_substruct_iterator;
-      typedef line_iterator        substruct_iterator;
-      typedef active_line_iterator active_substruct_iterator;    
+      typedef raw_line_iterator    raw_face_iterator;
+      typedef line_iterator        face_iterator;
+      typedef active_line_iterator active_face_iterator;    
     \end{verbatim}
 
     By using the cell iterators, you can write code nearly independent of
@@ -1023,9 +1023,9 @@ class Triangulation : public TriaDimensionInfo<dim> {
     typedef typename TriaDimensionInfo<dim>::cell_iterator cell_iterator;
     typedef typename TriaDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
 
-    typedef typename TriaDimensionInfo<dim>::raw_substruct_iterator raw_substruct_iterator;
-    typedef typename TriaDimensionInfo<dim>::substruct_iterator substruct_iterator;
-    typedef typename TriaDimensionInfo<dim>::active_substruct_iterator active_substruct_iterator;
+    typedef typename TriaDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+    typedef typename TriaDimensionInfo<dim>::face_iterator face_iterator;
+    typedef typename TriaDimensionInfo<dim>::active_face_iterator active_face_iterator;
 
                                     /**
                                      *  Create a triangulation and create
@@ -1338,6 +1338,109 @@ class Triangulation : public TriaDimensionInfo<dim> {
                                     //@}
 
                                     /*---------------------------------------*/
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Face iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first face, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no faces, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 2D and #begin_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    begin_raw_face   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used face
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 2D and #begin_quad# in 3D.
+                                     */
+    face_iterator        begin_face       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  face on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 2D and #begin_active_quad# in 3D.
+                                     */
+    active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 2D and #end_quad# in 3D.
+                                     */
+    raw_face_iterator    end_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last face, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  face of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face (const unsigned int level) const;
+                                    //@}
+
+    
+                                    /*---------------------------------------*/
 
                                     /**
                                      *  @name Line iterator functions
index 17a5873dc22c9fbe29318dd07a3dcd786d8e9e96..572239610dcba04f289b3c213509cb341c3408c5 100644 (file)
@@ -18,6 +18,8 @@ template <int dim> class FiniteElement;
 template <int dim> class Quadrature;
 template <int dim> class DataOut;
 template <int dim> class Function;
+template <int dim> class Equation;
+template <int dim> class Assembler;
 
 template <class Key, class T, class Compare = less<Key>, class Alloc = alloc> class map;
 
@@ -110,6 +112,16 @@ enum NormType {
   eliminating all coupling between this degree of freedom and others. Now
   also the column consists only of zeroes, apart from the main diagonal entry.
 
+  It seems as if we had to make clear not to overwrite the lines of other
+  boundary nodes when doing the Gauss elimination step. However, since we
+  reset the right hand side when passing such a node, it is not a problem
+  to change the right hand side values of other boundary nodes not yet
+  processed. It would be a problem to change those entries of nodes already
+  processed, but since the matrix entry of the present column on the row
+  of an already processed node is zero, the Gauss step does not change
+  the right hand side. We need therefore not take special care of other
+  boundary nodes.
+  
   To make solving faster, we preset the solution vector with the right boundary
   values. Since boundary nodes can never be hanging nodes, and since all other
   entries of the solution vector are zero, we need not condense the solution
@@ -123,7 +135,9 @@ enum NormType {
   since their search directions may contains components in the direction
   of the boundary node. For this reason, we perform a very simple line
   balancing by not setting the main diagonal entry to unity, but rather
-  to the value it had before deleting this line. Of course we have to change
+  to the value it had before deleting this line, or to the first nonzero
+  main diagonal entry if it is zero from a previous Gauss elimination
+  step. Of course we have to change
   the right hand side appropriately. This is not a very good
   strategy, but it at least should give the main diagonal entry a value
   in the right order of dimension, which makes the solving process a bit
@@ -135,6 +149,13 @@ enum NormType {
   not during solving iteratively, it may or may not be necessary to set
   the correct value after solving again. This question is an open one as of
   now and may be answered by future experience.
+
+  At present, boundary values are interpolated, i.e. a node is given the
+  point value of the boundary function. In some cases, it may be necessary
+  to use the L2-projection of the boundary function or any other method.
+  This can be done by overloading the virtual function
+  #make_boundary_value_list# which must return a list of boundary dofs
+  and their corresponding values.
   
 
   {\bf Computing errors}
@@ -194,7 +215,7 @@ class ProblemBase {
                                      * See the general documentation of this
                                      * class for more detail.
                                      */
-    typedef map<unsigned char,Function<dim> > DirichletBC;
+    typedef map<unsigned char,Function<dim>* > DirichletBC;
     
                                     /**
                                      * Constructor. Use this triangulation and
@@ -282,6 +303,18 @@ class ProblemBase {
                                      * want anything else.
                                      */
     virtual pair<char*,char*> get_solution_name () const;
+
+                                    /**
+                                     * Make up the list of node subject
+                                     * to Dirichlet boundary conditions
+                                     * and the values they are to be
+                                     * assigned.
+                                     *
+                                     * See the general doc for more
+                                     * information.
+                                     */
+    virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
+                                          map<int,double>   &boundary_values) const;
     
                                     /**
                                      * Exception
@@ -299,6 +332,10 @@ class ProblemBase {
                                      * Exception
                                      */
     DeclException0 (ExcNotImplemented);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidBoundaryIndicator);
     
   protected:
                                     /**
@@ -344,9 +381,9 @@ class ProblemBase {
                                      * as described in the general
                                      * documentation.
                                      */
-    void apply_dirichlet_bc (dSMatrix &matrix,
-                            dVector  &solution,
-                            dVector  &right_hand_side,
+    void apply_dirichlet_bc (dSMatrix          &matrix,
+                            dVector           &solution,
+                            dVector           &right_hand_side,
                             const DirichletBC &dirichlet_bc);
     
     friend class Assembler<dim>;
index 17e7817ad898306867f59a617fa04bb9fdb519e1..8de61c1d6e8845f5fb37260b4caf010ab1782e1e 100644 (file)
@@ -83,6 +83,23 @@ void DoFLineAccessor<dim,BaseClass>::set_vertex_dof_index (const unsigned int ve
 
 
 
+template <int dim, class BaseClass>
+void
+DoFLineAccessor<dim,BaseClass>::get_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
+  
+  dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
+      dof_indices.push_back (vertex_dof_index(vertex,d));
+  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
+    dof_indices.push_back (dof_index(d));
+};
+
+
+
 template <int dim, class BaseClass>
 TriaIterator<dim,DoFLineAccessor<dim,BaseClass> >
 DoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
@@ -177,6 +194,27 @@ void DoFQuadAccessor<dim,BaseClass>::set_vertex_dof_index (const unsigned int ve
 
 
 
+template <int dim, class BaseClass>
+void
+DoFQuadAccessor<dim,BaseClass>::get_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
+  
+  dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
+      dof_indices.push_back (vertex_dof_index(vertex,d));
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
+      dof_indices.push_back (this->line(line)->dof_index(d));
+  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_quad; ++d)
+    dof_indices.push_back (dof_index(d));
+};
+
+
+
+
 template <int dim, class BaseClass>
 TriaIterator<dim,DoFLineAccessor<dim,LineAccessor<dim> > >
 DoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
@@ -257,7 +295,7 @@ DoFCellAccessor<dim>::child (const unsigned int i) const {
 
 
 
-DoFSubstructAccessor<1>::substruct_iterator
+DoFSubstructAccessor<1>::face_iterator
 DoFCellAccessor<1>::face (const unsigned int) const {
   Assert (false, ExcNotUsefulForThisDimension());
   return 0;
@@ -265,50 +303,13 @@ DoFCellAccessor<1>::face (const unsigned int) const {
 
 
 
-DoFSubstructAccessor<2>::substruct_iterator
+DoFSubstructAccessor<2>::face_iterator
 DoFCellAccessor<2>::face (const unsigned int i) const {
   return line(i);
 };
 
 
 
-void
-DoFCellAccessor<1>::get_dof_indices (vector<int> &dof_indices) const {
-  Assert (dof_handler != 0, ExcInvalidObject());
-  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
-  Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
-  
-  dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
-  for (unsigned int vertex=0; vertex<2; ++vertex)
-    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
-      dof_indices.push_back (vertex_dof_index(vertex,d));
-  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
-    dof_indices.push_back (dof_index(d));
-};
-
-
-
-
-void
-DoFCellAccessor<2>::get_dof_indices (vector<int> &dof_indices) const {
-  Assert (dof_handler != 0, ExcInvalidObject());
-  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
-  Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
-  
-  dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
-  for (unsigned int vertex=0; vertex<4; ++vertex)
-    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
-      dof_indices.push_back (vertex_dof_index(vertex,d));
-  for (unsigned int line=0; line<4; ++line)
-    for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
-      dof_indices.push_back (this->line(line)->dof_index(d));
-  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_quad; ++d)
-    dof_indices.push_back (dof_index(d));
-};
-
-
-
-
 void
 DoFCellAccessor<1>::get_dof_values (const dVector  &values,
                                    vector<double> &dof_values) const {
index 1ec514f779bd89a801f4f4435e958ccf147bc182..ba2d1eaf33cf501799745e5a997ec63b765591a1 100644 (file)
@@ -87,7 +87,7 @@ void ConstraintMatrix::clear () {
 void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
                                 dSMatrixStruct       &condensed) const {
   Assert (sorted == true, ExcMatrixNotClosed());
-  Assert (uncondensed.compressed == true, ExcMatrixNotClosed());
+  Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed());
   Assert (uncondensed.n_rows() == uncondensed.n_cols(),
          ExcMatrixNotSquare());
 
@@ -116,21 +116,23 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
 
  
   next_constraint = lines.begin();
-  for (int row=0; row<uncondensed.rows; ++row)
+  for (int row=0; row<uncondensed.n_rows(); ++row)
     if (new_line[row] != -1)
                                       // line not constrained
                                       // copy entries if column will not
                                       // be condensed away, distribute
                                       // otherwise
-      for (int j=uncondensed.rowstart[row]; j<uncondensed.rowstart[row+1]; ++j)
-       if (new_line[uncondensed.colnums[j]] != -1)
-         condensed.add (new_line[row], new_line[uncondensed.colnums[j]]);
+      for (int j=uncondensed.get_rowstart_indices()[row];
+          j<uncondensed.get_rowstart_indices()[row+1]; ++j)
+       if (new_line[uncondensed.get_column_numbers()[j]] != -1)
+         condensed.add (new_line[row], new_line[uncondensed.get_column_numbers()[j]]);
        else 
          {
                                             // let c point to the constraint
                                             // of this column
            vector<ConstraintLine>::const_iterator c = lines.begin();
-           while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c;
+           while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j])
+             ++c;
 
            for (unsigned int q=0; q!=(*c).entries.size(); ++q) 
              condensed.add (new_line[row], new_line[(*c).entries[q].first]);
@@ -138,13 +140,14 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
     else
                                       // line must be distributed
       {
-       for (int j=uncondensed.rowstart[row]; j<uncondensed.rowstart[row+1]; ++j)
+       for (int j=uncondensed.get_rowstart_indices()[row];
+            j<uncondensed.get_rowstart_indices()[row+1]; ++j)
                                           // for each entry: distribute
-         if (new_line[uncondensed.colnums[j]] != -1)
+         if (new_line[uncondensed.get_column_numbers()[j]] != -1)
                                             // column is not constrained
            for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q) 
                condensed.add (new_line[(*next_constraint).entries[q].first],
-                              new_line[uncondensed.colnums[j]]);
+                              new_line[uncondensed.get_column_numbers()[j]]);
        
          else
                                             // not only this line but
@@ -153,7 +156,7 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
                                               // let c point to the constraint
                                               // of this column
              vector<ConstraintLine>::const_iterator c = lines.begin();
-             while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c;
+             while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j]) ++c;
              
              for (unsigned int p=0; p!=(*c).entries.size(); ++p)
                for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
@@ -171,7 +174,7 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
 
 void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const {
   Assert (sorted == true, ExcMatrixNotClosed());
-  Assert (sparsity.compressed == false, ExcMatrixIsClosed());
+  Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
   Assert (sparsity.n_rows() == sparsity.n_cols(),
          ExcMatrixNotSquare());
   
@@ -199,46 +202,52 @@ void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const {
   for (int row=0; row<n_rows; ++row)
     if (distribute[row] == -1)
                                       // regular line. loop over cols
-      for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+      for (int j=sparsity.get_rowstart_indices()[row];
+          j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                         // end of row reached?
-       if (sparsity.colnums[j] == -1)
+       if (sparsity.get_column_numbers()[j] == -1)
          break;
        else
          {
-           if (distribute[sparsity.colnums[j]] != -1)
+           if (distribute[sparsity.get_column_numbers()[j]] != -1)
                                               // distribute entry at regular
                                               // row #row# and irregular column
                                               // sparsity.colnums[j]
              for (unsigned int q=0;
-                  q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q) 
+                  q!=lines[distribute[sparsity.get_column_numbers()[j]]].entries.size();
+                  ++q) 
                sparsity.add (row,
-                             lines[distribute[sparsity.colnums[j]]].entries[q].first);
+                             lines[distribute[sparsity.get_column_numbers()[j]]].
+                             entries[q].first);
          }
     else
                                       // row must be distributed
-      for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+      for (int j=sparsity.get_rowstart_indices()[row];
+          j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                         // end of row reached?
-       if (sparsity.colnums[j] == -1)
+       if (sparsity.get_column_numbers()[j] == -1)
          break;
        else
          {
-           if (distribute[sparsity.colnums[j]] == -1)
+           if (distribute[sparsity.get_column_numbers()[j]] == -1)
                                               // distribute entry at irregular
                                               // row #row# and regular column
                                               // sparsity.colnums[j]
              for (unsigned int q=0;
                   q!=lines[distribute[row]].entries.size(); ++q) 
                sparsity.add (lines[distribute[row]].entries[q].first,
-                             sparsity.colnums[j]);
+                             sparsity.get_column_numbers()[j]);
            else
                                               // distribute entry at irregular
                                               // row #row# and irregular column
-                                              // sparsity.colnums[j]
+                                              // sparsity.get_column_numbers()[j]
              for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
                for (unsigned int q=0;
-                    q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+                    q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+                                   .entries.size(); ++q)
                  sparsity.add (lines[distribute[row]].entries[p].first,
-                               lines[distribute[sparsity.colnums[j]]].entries[q].first);
+                               lines[distribute[sparsity.get_column_numbers()[j]]]
+                               .entries[q].first);
          };
   sparsity.compress();
 };
@@ -250,8 +259,8 @@ void ConstraintMatrix::condense (const dSMatrix &uncondensed,
   const dSMatrixStruct &uncondensed_struct = uncondensed.get_sparsity_pattern ();
   
   Assert (sorted == true, ExcMatrixNotClosed());
-  Assert (uncondensed_struct.compressed == true, ExcMatrixNotClosed());
-  Assert (condensed.get_sparsity_pattern().compressed == true, ExcMatrixNotClosed());
+  Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed());
+  Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed());
   Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(),
          ExcMatrixNotSquare());
   Assert (condensed.n() == condensed.m(),
@@ -284,40 +293,43 @@ void ConstraintMatrix::condense (const dSMatrix &uncondensed,
 
  
   next_constraint = lines.begin();
-  for (int row=0; row<uncondensed_struct.rows; ++row)
+  for (int row=0; row<uncondensed_struct.n_rows(); ++row)
     if (new_line[row] != -1)
                                       // line not constrained
                                       // copy entries if column will not
                                       // be condensed away, distribute
                                       // otherwise
-      for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
-       if (new_line[uncondensed_struct.colnums[j]] != -1)
-         condensed.add (new_line[row], new_line[uncondensed_struct.colnums[j]],
-                        uncondensed.val[j]);
+      for (int j=uncondensed_struct.get_rowstart_indices()[row];
+          j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
+       if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
+         condensed.add (new_line[row], new_line[uncondensed_struct.get_column_numbers()[j]],
+                        uncondensed.global_entry(j));
        else 
          {
                                             // let c point to the constraint
                                             // of this column
            vector<ConstraintLine>::const_iterator c = lines.begin();
-           while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+           while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c;
 
            for (unsigned int q=0; q!=(*c).entries.size(); ++q)
                                               // distribute to rows with
                                               // appropriate weight
              condensed.add (new_line[row], new_line[(*c).entries[q].first],
-                            uncondensed.val[j] * (*c).entries[q].second);
+                            uncondensed.global_entry(j) * (*c).entries[q].second);
          }
     else
                                       // line must be distributed
       {
-       for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
+       for (int j=uncondensed_struct.get_rowstart_indices()[row];
+            j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
                                           // for each column: distribute
-         if (new_line[uncondensed_struct.colnums[j]] != -1)
+         if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
                                             // column is not constrained
            for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q) 
                condensed.add (new_line[(*next_constraint).entries[q].first],
-                              new_line[uncondensed_struct.colnums[j]],
-                              uncondensed.val[j] * (*next_constraint).entries[q].second);
+                              new_line[uncondensed_struct.get_column_numbers()[j]],
+                              uncondensed.global_entry(j) *
+                              (*next_constraint).entries[q].second);
        
          else
                                             // not only this line but
@@ -326,13 +338,13 @@ void ConstraintMatrix::condense (const dSMatrix &uncondensed,
                                               // let c point to the constraint
                                               // of this column
              vector<ConstraintLine>::const_iterator c = lines.begin();
-             while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+             while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c;
              
              for (unsigned int p=0; p!=(*c).entries.size(); ++p)
                for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
                    condensed.add (new_line[(*next_constraint).entries[q].first],
                                   new_line[(*c).entries[p].first],
-                                  uncondensed.val[j] *
+                                  uncondensed.global_entry(j) *
                                   (*next_constraint).entries[q].second *
                                   (*c).entries[p].second);
            };
@@ -347,7 +359,7 @@ void ConstraintMatrix::condense (dSMatrix &uncondensed) const {
   const dSMatrixStruct &sparsity = uncondensed.get_sparsity_pattern ();
 
   Assert (sorted == true, ExcMatrixNotClosed());
-  Assert (sparsity.compressed == true, ExcMatrixNotClosed());
+  Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
   Assert (sparsity.n_rows() == sparsity.n_cols(),
          ExcMatrixNotSquare());
   
@@ -378,67 +390,75 @@ void ConstraintMatrix::condense (dSMatrix &uncondensed) const {
       
     if (distribute[row] == -1)
                                       // regular line. loop over cols
-      for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+      for (int j=sparsity.get_rowstart_indices()[row];
+          j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                         // end of row reached?
-       if (sparsity.colnums[j] == -1)
+       if (sparsity.get_column_numbers()[j] == -1)
          break;
        else
          {
-           if (distribute[sparsity.colnums[j]] != -1)
+           if (distribute[sparsity.get_column_numbers()[j]] != -1)
                                               // distribute entry at regular
                                               // row #row# and irregular column
-                                              // sparsity.colnums[j]; set old
+                                              // sparsity.get_column_numbers()[j]; set old
                                               // entry to zero
              {
                for (unsigned int q=0;
-                    q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q) 
+                    q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+                                   .entries.size(); ++q) 
                  uncondensed.add (row,
-                                  lines[distribute[sparsity.colnums[j]]].entries[q].first,
-                                  uncondensed.val[j] *
-                                  lines[distribute[sparsity.colnums[j]]].entries[q].second);
-               uncondensed.val[j] = 0.;
+                                  lines[distribute[sparsity.get_column_numbers()[j]]]
+                                  .entries[q].first,
+                                  uncondensed.global_entry(j) *
+                                  lines[distribute[sparsity.get_column_numbers()[j]]]
+                                  .entries[q].second);
+               uncondensed.global_entry(j) = 0.;
              };
          }
     else
                                       // row must be distributed
-      for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+      for (int j=sparsity.get_rowstart_indices()[row];
+          j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                         // end of row reached?
-       if (sparsity.colnums[j] == -1)
+       if (sparsity.get_column_numbers()[j] == -1)
          break;
        else
          {
-           if (distribute[sparsity.colnums[j]] == -1)
+           if (distribute[sparsity.get_column_numbers()[j]] == -1)
                                               // distribute entry at irregular
                                               // row #row# and regular column
-                                              // sparsity.colnums[j]. set old
+                                              // sparsity.get_column_numbers()[j]. set old
                                               // entry to zero
              {
                for (unsigned int q=0;
                     q!=lines[distribute[row]].entries.size(); ++q) 
                  uncondensed.add (lines[distribute[row]].entries[q].first,
-                                  sparsity.colnums[j],
-                                  uncondensed.val[j] *
+                                  sparsity.get_column_numbers()[j],
+                                  uncondensed.global_entry(j) *
                                   lines[distribute[row]].entries[q].second);
                
-               uncondensed.val[j] = 0.;
+               uncondensed.global_entry(j) = 0.;
              }
            else
                                               // distribute entry at irregular
                                               // row #row# and irregular column
-                                              // sparsity.colnums[j]
+                                              // sparsity.get_column_numbers()[j]
                                               // set old entry to one if on main
                                               // diagonal, zero otherwise
              {
                for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
                  for (unsigned int q=0;
-                      q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+                      q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+                                     .entries.size(); ++q)
                    uncondensed.add (lines[distribute[row]].entries[p].first,
-                                    lines[distribute[sparsity.colnums[j]]].entries[q].first,
-                                    uncondensed.val[j] *
+                                    lines[distribute[sparsity.get_column_numbers()[j]]]
+                                    .entries[q].first,
+                                    uncondensed.global_entry(j) *
                                     lines[distribute[row]].entries[p].second *
-                                    lines[distribute[sparsity.colnums[j]]].entries[q].second);
+                                    lines[distribute[sparsity.get_column_numbers()[j]]]
+                                    .entries[q].second);
                
-               uncondensed.val[j] = (row == sparsity.colnums[j] ?
+               uncondensed.global_entry(j) = (row == sparsity.get_column_numbers()[j] ?
                                      1. : 0. );
              };
          };
index a47282da07c2de6abbf5cd2fbb2c1c7a1365c5a4..53002bf983829af1a8a3a62cee31b32a49f555d5 100644 (file)
@@ -174,6 +174,164 @@ DoFHandler<2>::last_active (const unsigned int level) const {
 
 //------------------------------------------------------------------
 
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::begin_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::begin_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::begin_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::end_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::last_raw_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::last_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::last_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::last_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::last_active_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::last_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::begin_face (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::begin_active_face (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::end_face () const {
+  return end_line ();
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::last_raw_face () const {
+  return last_raw_line ();
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::last_raw_face (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::last_face () const {
+  return last_line ();
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::last_face (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::last_active_face () const {
+  return last_active_line ();
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::last_active_face (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+
+
+// ---------------------------------------------------------------
+
+
+
 template <int dim>
 typename DoFHandler<dim>::raw_line_iterator
 DoFHandler<dim>::begin_raw_line (const unsigned int level) const {
index fe955a0547b639ed5019f614c18ec573ad39fecc..613077495777061132f46eb97564747d2faa7c74 100644 (file)
@@ -956,6 +956,161 @@ Triangulation<2>::last_active (const unsigned int level) const {
 
 
 
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::begin_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::begin_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::begin_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::end_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::last_raw_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::last_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::last_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::last_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::last_active_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::last_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::begin_face (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::begin_active_face (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::end_face () const {
+  return end_line ();
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::last_raw_face () const {
+  return last_raw_line ();
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::last_raw_face (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::last_face () const {
+  return last_line ();
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::last_face (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::last_active_face () const {
+  return last_active_line ();
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::last_active_face (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+
+
+
+
 template <int dim>
 typename TriaDimensionInfo<dim>::raw_line_iterator
 Triangulation<dim>::begin_raw_line (unsigned int level) const {
index a6228c505c157857c04b77a8459b3ef38793b4dc..ff74f023cdb860a0000b2a874ada581a0a3e4fe6 100644 (file)
@@ -529,7 +529,7 @@ void CellAccessor<1>::set_material_id (const unsigned char mat_id) const {
 
 
 
-Triangulation<1>::substruct_iterator
+Triangulation<1>::face_iterator
 CellAccessor<1>::face (const unsigned int) const {
   Assert (false, ExcNotUsefulForThisDimension());
   return 0;
@@ -563,7 +563,7 @@ void CellAccessor<2>::set_material_id (const unsigned char mat_id) const {
 
 
 
-Triangulation<2>::substruct_iterator
+Triangulation<2>::face_iterator
 CellAccessor<2>::face (const unsigned int i) const {
   return line(i);
 };
index 08e75e5239f4cd9ef9e432c45a124e2f3e75f1db..84c946b70867d3d13f0db234da3dab48fb71461f 100644 (file)
@@ -11,6 +11,7 @@
 #include "../../../mia/vectormemory.h"
 #include "../../../mia/cg.h"
 
+#include <map.h>
 #include <algo.h>
 
 
@@ -273,7 +274,134 @@ void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
                                           dVector  &solution,
                                           dVector  &right_hand_side,
                                           const DirichletBC &dirichlet_bc) {
-//  l;
+  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
+         ExcInvalidBoundaryIndicator());
+  
+                                  // first make up a list of dofs subject
+                                  // to any boundary condition and which
+                                  // value they take; if a node occurs
+                                  // with two bc (e.g. a corner node, with
+                                  // the lines in 2D being subject to
+                                  // different bc's), the last value is taken
+  map<int,double> boundary_values;
+  make_boundary_value_list (dirichlet_bc, boundary_values);
+
+  map<int,double>::const_iterator dof, endd;
+  const unsigned int n_dofs   = (unsigned int)matrix.m();
+  const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern();
+  const int *sparsity_rowstart = sparsity.get_rowstart_indices(),
+           *sparsity_colnums  = sparsity.get_column_numbers();
+
+  for (dof=boundary_values.begin(), endd=boundary_values.end(); dof != endd; ++dof)
+    {
+                                      // for each boundary dof:
+
+                                      // set entries of this line
+                                      // to zero
+      for (int j=sparsity_rowstart[(*dof).first];
+          j<sparsity_rowstart[(*dof).first+1]; ++j)
+       if (sparsity_colnums[j] != (*dof).first)
+                                          // if not main diagonal entry
+         matrix.global_entry(j) = 0.;
+      
+                                      // set right hand side to
+                                      // wanted value: if main diagonal
+                                      // entry nonzero, don't touch it
+                                      // and scale rhs accordingly. If
+                                      // zero, take the first main
+                                      // diagonal entry we can find, or
+                                      // one if no nonzero main diagonal
+                                      // element exists.
+      if (matrix.diag_element((*dof).first) != 0.0)
+       right_hand_side((*dof).first) = (*dof).second *
+                                       matrix.diag_element((*dof).first);
+      else
+       {
+         double first_diagonal_entry = 1;
+         for (unsigned int i=0; i<n_dofs; ++i)
+           if (matrix.diag_element(i) != 0)
+             {
+               first_diagonal_entry = matrix.diag_element(i);
+               break;
+             };
+         
+         matrix.set((*dof).first, (*dof).first,
+                    first_diagonal_entry);
+         right_hand_side((*dof).first) = (*dof).second * first_diagonal_entry;
+       };
+
+                                      // store the only nonzero entry
+                                      // of this line for the Gauss
+                                      // elimination step
+      const double diagonal_entry = matrix.diag_element((*dof).first);
+      
+
+                                      // do the Gauss step
+      for (unsigned int row=0; row<n_dofs; ++row)
+       for (int j=sparsity_rowstart[row];
+            j<sparsity_rowstart[row+1]; ++row)
+         if ((sparsity_colnums[j] == (signed int)(*dof).first) &&
+             ((signed int)row != (*dof).first))
+                                            // this line has an entry
+                                            // in the regarding column
+                                            // but this is not the main
+                                            // diagonal entry
+           {
+                                              // correct right hand side
+             right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
+                                     (*dof).second;
+             
+                                              // set matrix entry to zero
+             matrix.global_entry(j) = 0.;
+           };
+      
+      
+                                      // preset solution vector
+      solution((*dof).first) = (*dof).second;
+    };
+};
+
+
+
+
+
+void
+ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
+                                         map<int,double>   &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
+template <int dim>
+void
+ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
+                                           map<int,double>   &boundary_values) const {
+  DoFHandler<dim>::active_face_iterator face = dof_handler->begin_active_face(),
+                                       endf = dof_handler->end_face();
+  DirichletBC::const_iterator function_ptr;
+  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
+       vector<int>         face_dofs;
+       vector<Point<dim> > dof_locations;
+       vector<double>      dof_values;
+       
+       face->get_dof_indices (face_dofs);
+//physical location
+       (*function_ptr).second->value_list (dof_locations, dof_values);
+
+                                        // enter into list
+       for (unsigned int i=0; i<face_dofs.size(); ++i) 
+         boundary_values[face_dofs[i]] = dof_values[i];
+      };
 };
 
   
index d7761bbc5995dcff6c812c18027240cf340e05a7..a1f644440822e5915d8d0670e5f213ecb7061aef 100644 (file)
@@ -9,12 +9,15 @@
 #include <grid/tria_boundary.h>
 #include <grid/dof_constraints.h>
 #include <basic/data_io.h>
+#include <basic/function.h>
 #include <fe/fe_lib.h>
 #include <fe/quadrature_lib.h>
 #include <numerics/base.h>
 #include <numerics/assembler.h>
 #include <lac/dsmatrix.h>
 
+
+#include <map.h>
 #include <fstream.h>
 #include <cmath>
 extern "C" {

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.