]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support for error integration and initial support for dirichlet boundary conditions.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 Mar 1998 17:48:41 +0000 (17:48 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 Mar 1998 17:48:41 +0000 (17:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@95 0785d39b-7218-0410-832d-ea1e28bc413d

15 files changed:
deal.II/base/include/base/function.h
deal.II/deal.II/Makefile
deal.II/deal.II/Todo
deal.II/deal.II/include/dofs/dof_accessor.h
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/grid/tria_accessor.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/include/numerics/data_io.h
deal.II/deal.II/source/dofs/dof_accessor.cc
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc

index a00e86d2a634cbba93f316758b4aabfbc3646d7d..aac5e2233612b5c5467d45544588f37f2a8a585a 100644 (file)
@@ -7,7 +7,7 @@
 
 //forward declaration
 template <int dim> class Point;
-template <class T> class vector;
+template <class T, class Alloc=alloc> class vector;
 
 
 
index ebf512092e7a7926c7a2615ee0c250f10fcf6fba..1de1718ad3deeceee24b38dacf6504b32ba08622 100644 (file)
@@ -3,7 +3,7 @@
 
 all: deal.II doc examples
 
-clean: lib-clean doc-clean examples-clean
+clean: lib-clean examples-clean doc-clean
 
 
 deal.II:
index b0d40ed8f4befdd2989b137232f9ef70170aede8..35ee247952998730e312a1002fb1ceab1eca6f78 100644 (file)
@@ -36,15 +36,21 @@ Let ParameterHandler and DataIn/Out throw exceptions. Make more tests
 Replace function objects by mem_fun, fun_ptr, ... when member function
   templates are supported.
 
+Invent a way to tell FEValues::reinit which fields to reinit. If you
+  only want to integrate a FE solution over a cell, transforming the
+  gradients may be a rather expensive operation!
 
+  Check in FEValues::shape_grad whether the gardients have been
+  reinit'd last time to avoid people taking gradients which have not
+  been computed.
+
+implement ProblemBase::distribute_cell_to_dof_vector
 
 
-DEAL:
-make dSMatrixStruct warn if max # of entries per row is overrun!
 
-let dSMatrixStruct have an empty constructor to allow for member
-  objects being initialized later
 
+
+DEAL:
 why do loops in dFMatrix run backwards?
 
 Let an empty dfmatrix assume it has dimensions zero rather than one!
index aa74ccb7a408fbd2ba534d1726666ae1b11a17b6..63df424b3991e30599c92054244124a60bb17713 100644 (file)
@@ -7,6 +7,8 @@
 
 template <int dim> class Triangulation;
 template <int dim> class DoFHandler;
+class dVector;
+
 
 #include <grid/tria_accessor.h>
 #include <vector.h>
@@ -239,7 +241,8 @@ class DoFQuadAccessor :  public BaseClass, public DoFAccessor<dim> {
                                      *  Return a pointer to the #i#th line
                                      *  bounding this #Quad#.
                                      */
-    TriaIterator<dim,DoFLineAccessor<dim,LineAccessor<dim> > > line (const unsigned int i) const;
+    TriaIterator<dim,DoFLineAccessor<dim,LineAccessor<dim> > >
+    line (const unsigned int i) const;
 
                                     /**
                                      * Return the #i#th child as a DoF quad
@@ -299,10 +302,17 @@ class DoFSubstructAccessor<1> :  public DoFLineAccessor<1,CellAccessor<1> > {
                                      * Constructor
                                      */
     DoFSubstructAccessor (Triangulation<1> *tria,
-                                const int         level,
-                                const int         index,
-                                const void       *local_data) :
+                         const int         level,
+                         const int         index,
+                         const void       *local_data) :
                    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
+                                    // it would have to ask the DoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef void * substruct_iterator;
 };
 
 
@@ -317,10 +327,17 @@ class DoFSubstructAccessor<2> : public DoFQuadAccessor<2,CellAccessor<2> > {
                                      * Constructor
                                      */
     DoFSubstructAccessor (Triangulation<2> *tria,
-                                const int         level,
-                                const int         index,
-                                const void       *local_data) :
+                         const int         level,
+                         const int         index,
+                         const void       *local_data) :
                    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
+                                    // 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;
 };
 
 
@@ -369,6 +386,16 @@ class DoFCellAccessor :  public DoFSubstructAccessor<dim> {
                                      */
     TriaIterator<dim,DoFCellAccessor<dim> > child (const unsigned int) const;
 
+                                    /**
+                                     * Return an iterator to the #i#th face
+                                     * of this cell.
+                                     *
+                                     * This function is not implemented in 1D,
+                                     * and maps to DoFQuadAccessor::line in 2D.
+                                     */
+    typename DoFSubstructAccessor<dim>::substruct_iterator
+    face (const unsigned int i) const;
+
                                     /**
                                      * Return the indices of the dofs of this
                                      * cell in the standard ordering: dofs
@@ -377,10 +404,38 @@ class DoFCellAccessor :  public DoFSubstructAccessor<dim> {
                                      * dofs on quad 0, etc.
                                      *
                                      * The dof indices are appended to the end
-                                     * of the vector. It is not assumed that
+                                     * of the vector. It is assumed that
                                      * the vector be empty beforehand.
                                      */
-    void dof_indices (vector<int> &dof_indices) const;
+    void get_dof_indices (vector<int> &dof_indices) const;
+
+                                    /**
+                                     * Return the value of the given vector
+                                     * restricted to 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 values are appended to the end
+                                     * of the vector. It is assumed that
+                                     * the vector be empty beforehand.
+                                     */
+    void get_dof_values (const dVector  &values,
+                        vector<double> &dof_values) const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcVectorNotEmpty);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcVectorDoesNotMatch);
+                                    /**
+                                     *  Exception
+                                     */
+    DeclException0 (ExcNotUsefulForThisDimension);
 };
 
 
index f660cb4ad3788c5ab277c6d3938a063e06ef29f0..360b6f434c83e0abc9262b719b3fe7cbb5b49913 100644 (file)
@@ -152,6 +152,10 @@ class DoFDimensionInfo<1> {
     typedef raw_line_iterator    raw_cell_iterator;
     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;
 };
 
 
@@ -177,6 +181,10 @@ class DoFDimensionInfo<2> {
     typedef raw_quad_iterator    raw_cell_iterator;
     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;    
 };
 
 
@@ -428,7 +436,11 @@ 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;
 
+    
                                     /**
                                      * Constructor. Take #tria# as the
                                      * triangulation to work on.
index 7757b1148f27e97be582e47b42ecba8974951eb1..f74d3e841400d8de9b404f6b142266f9dcfffefc 100644 (file)
@@ -72,6 +72,15 @@ class FEValues {
     double shape_value (const unsigned int i,
                        const unsigned int j) const;
 
+                                    /**
+                                     * Return a pointer to the matrix holding
+                                     * all values of shape functions at all
+                                     * integration points, on the present cell.
+                                     * For the format of this matrix, see the
+                                     * documentation for the matrix itself.
+                                     */
+    const dFMatrix & get_shape_values () const;
+    
                                     /**
                                      * Return the gradient of the #i#th shape
                                      * function at the #j# quadrature point.
@@ -88,18 +97,40 @@ class FEValues {
     const Point<dim> & shape_grad (const unsigned int i,
                                   const unsigned int j) const;
 
+                                    /** 
+                                     * Return a pointer to the matrix holding
+                                     * all gradients of shape functions at all
+                                     * integration points, on the present cell.
+                                     * For the format of this matrix, see the
+                                     * documentation for the matrix itself.
+                                     */
+    const vector<vector<Point<dim> > > & get_shape_grads () const;
+    
                                     /**
                                      * Return the position of the #i#th
                                      * quadrature point in real space.
                                      */
     const Point<dim> & quadrature_point (const unsigned int i) const;
 
+                                    /**
+                                     * Return a pointer to the vector of
+                                     * quadrature points.
+                                     */
+    const vector<Point<dim> > & get_quadrature_points () const;
+    
                                     /**
                                      * Return the Jacobi determinant times
                                      * the weight of the #i#th quadrature
                                      * point.
                                      */
     double JxW (const unsigned int i) const;
+
+                                    /**
+                                     * Return a pointer to the array holding
+                                     * the JxW values at the different
+                                     * quadrature points.
+                                     */
+    const vector<double> & get_JxW_values () const;
     
                                     /**
                                      * Reinitialize the gradients, Jacobi
@@ -641,6 +672,48 @@ class FiniteElement<2> : public FiniteElementBase<2> {
 
 
 
+
+
+/*------------------------ Inline functions -----------------------------------*/
+
+
+
+template <int dim>
+inline
+const dFMatrix & FEValues<dim>::get_shape_values () const {
+  return shape_values;
+};
+
+
+
+
+template <int dim>
+inline
+const vector<vector<Point<dim> > > &
+FEValues<dim>::get_shape_grads () const {
+  return shape_gradients;
+};
+
+
+
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEValues<dim>::get_quadrature_points () const {
+  return quadrature_points;
+};
+
+
+
+template <int dim>
+inline
+const vector<double> &
+FEValues<dim>::get_JxW_values () const {
+  return JxW_values;
+};
+
+
+
   
 /*----------------------------   fe.h     ---------------------------*/
 /* end of #ifndef __fe_H */
index a7fd1731199cc94ab965d777bfde57aef1491d63..2685ec001e8240206c1df28f4aee1f7ff359f605 100644 (file)
@@ -462,6 +462,10 @@ class TriaDimensionInfo;
     #Triangulation<1>#, the #quad_iterator#s are declared as
     #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
+    have no substructures apart from vertices, which are handled in
+    a different way, however.
     */
 class TriaDimensionInfo<1> {
   public:
@@ -476,6 +480,10 @@ class TriaDimensionInfo<1> {
     typedef raw_line_iterator    raw_cell_iterator;
     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;
 };
 
 
@@ -499,6 +507,10 @@ class TriaDimensionInfo<1> {
       typedef raw_quad_iterator    raw_cell_iterator;
       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;    
     \end{verbatim}
     */
 class TriaDimensionInfo<2> {
@@ -514,6 +526,10 @@ class TriaDimensionInfo<2> {
     typedef raw_quad_iterator    raw_cell_iterator;
     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;    
 };
 
 
@@ -584,10 +600,17 @@ class TriaDimensionInfo<2> {
       typedef quad_line_iterator   raw_cell_iterator;    
       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;    
     \end{verbatim}
 
     By using the cell iterators, you can write code nearly independent of
-    the spatial dimension.
+    the spatial dimension. The same applies for substructure iterators,
+    where a substructure is defined as a face of a cell. The face of a
+    cell is be a vertex in 1D and a line in 2D; however, vertices are
+    handled in a different way and therefore lines have no faces.
 
     The #Triangulation# class offers functions like #begin_active# which gives
     you an iterator to the first active cell. There are quite a lot of functions
@@ -1000,6 +1023,10 @@ 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;
+
                                     /**
                                      *  Create a triangulation and create
                                      *  the first level of the hierarchy.
index bff8f3d57f8e3896f3e4a3046a6f37a4c94bd2fb..14f4040a1c5dc69a140ab007969e6354cc595997 100644 (file)
@@ -731,10 +731,18 @@ class TriaSubstructAccessor<1> :  public LineAccessor<1> {
                                      * Constructor
                                      */
     TriaSubstructAccessor (Triangulation<1> *tria,
-                                const int         level,
-                                const int         index,
-                                const void       *local_data) :
+                          const int         level,
+                          const int         index,
+                          const void       *local_data) :
                    LineAccessor<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
+                                    // it would have to ask the DoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef void * substruct_iterator;
 };
 
 
@@ -749,10 +757,18 @@ class TriaSubstructAccessor<2> : public QuadAccessor<2> {
                                      * Constructor
                                      */
     TriaSubstructAccessor (Triangulation<2> *tria,
-                                const int         level,
-                                const int         index,
-                                const void       *local_data) :
+                          const int         level,
+                          const int         index,
+                          const void       *local_data) :
                    QuadAccessor<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
+                                    // it would have to ask the DoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef TriaIterator<2,LineAccessor<2> > substruct_iterator;
 };
 
 
@@ -853,6 +869,16 @@ class CellAccessor :  public TriaSubstructAccessor<dim> {
                                      */
     TriaIterator<dim,CellAccessor<dim> > child (const unsigned int i) const;
 
+                                    /**
+                                     * Return an iterator to the #i#th face
+                                     * of this cell.
+                                     *
+                                     * This function is not implemented in 1D,
+                                     * and maps to QuadAccessor::line in 2D.
+                                     */
+    typename TriaSubstructAccessor<dim>::substruct_iterator
+    face (const unsigned int i) const;
+    
                                     /**
                                      * Return the material id of this cell.
                                      */
index 7b2772ef41eb25513c0617ea58152bc38fa5a130..17a5873dc22c9fbe29318dd07a3dcd786d8e9e96 100644 (file)
@@ -10,6 +10,7 @@
 #include <grid/dof_constraints.h>
 
 
+
 // forward declaration
 template <int dim> class Triangulation;
 template <int dim> class DoFHandler;
@@ -18,6 +19,30 @@ template <int dim> class Quadrature;
 template <int dim> class DataOut;
 template <int dim> class Function;
 
+template <class Key, class T, class Compare = less<Key>, class Alloc = alloc> class map;
+
+
+
+
+/**
+   Denote which norm/integral is to be computed. The following possibilities
+   are implemented:
+   \begin{itemize}
+   \item #mean#: the function or difference of functions is integrated
+     on each cell.
+   \item #L1_norm#: the absolute value of the function is integrated.
+   \item #L2_norm#: the square of the function is integrated on each
+     cell; afterwards the root is taken of this value.
+   \end{itemize}
+*/
+enum NormType {
+      mean,
+      L1_norm,
+      L2_norm,
+      Linfty_norm
+};
+
+
 
 
 /**
@@ -37,12 +62,15 @@ template <int dim> class Function;
   The #assemble# member function does the assemblage of the system matrix and
   the given number of right hand sides. It does the following steps:
   \begin{itemize}
+    \item Initialize solution vector with zero entries.
     \item Create sparsity pattern of the system matrix and condense it with
       the constraints induced by hanging nodes.
     \item Initialize an assembler object.
     \item Loop over all cells and assemble matrix and vectors using the given
       quadrature formula and the equation object which contains the weak
       formulation of the equation.
+    \item Apply Dirichlet boundary conditions. See the section on boundary
+      conditions for more details.
     \item Condense the system matrix and right hand side with the constraints
       induced by hanging nodes.
   \end{itemize}
@@ -54,10 +82,120 @@ template <int dim> class Function;
   which results after having called the #assemble# function is solved. After
   this, the solution vector is distributed again, i.e. the constrained nodes
   are given their correct values.
+
+
+  {\bf Boundary conditions}
+
+  During assemblage of matrices and right hand side, use is made of dirichlet
+  boundary conditions (in short: bc) specified to the #assemble# function. You
+  can specify a list of pairs of boundary indicators (of type #unsigned char#;
+  see the section in the documentation of the \Ref{Triangulation} class for more
+  details) and the according functions denoting the dirichlet boundary values
+  of the nodes on boundary faces with this boundary indicator.
+
+  Usually, all other boundary conditions, such as inhomogeneous Neumann values
+  or mixed boundary conditions are handled in the weak formulation. No attempt
+  is made to include this into the process of assemblage therefore.
+
+  The inclusion into the assemblage process is as follows: when the matrix and
+  vectors are set up, a list of nodes subject to dirichlet bc is made and
+  matrix and vectors are changed accordingly. This is done by deleting all
+  entries in the matrix in the line of this degree of freedom, setting the
+  main diagonal entry to one and the right hand side element to the
+  boundary value at this node. This forces this node's value to be as specified.
+  To decouple the remaining linear system of equations and to make the system
+  symmetric again (at least if it was before), one Gauss elimination
+  step is performed with this line, by adding this (now almost empty) line to
+  all other lines which couple with the given degree of freedom and thus
+  eliminating all coupling between this degree of freedom and others. Now
+  also the column consists only of zeroes, apart from the main diagonal entry.
+
+  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
+  vector if the condensation process is done in-place. If done by copying
+  matrix and vectors to smaller ones, it would also be necessary to condense
+  the solution vector to preserve the preset boundary values.
+  
+  It it not clear whether the deletion of coupling between the boundary degree
+  of freedom and other dofs really forces the corresponding entry in the
+  solution vector to have the right value when using iterative solvers,
+  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
+  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
+  more stable. A refined algorithm would set the entry to the mean of the
+  other diagonal entries, but this seems to be too expensive.
+
+  Because of the mentioned question, whether or not a preset solution value
+  which does not couple with other degrees of freedom remains its value or
+  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.
+  
+
+  {\bf Computing errors}
+
+  The function #integrate_difference# performs the calculation of the error
+  between the finite element solution and a given (continuous) reference
+  function in different norms. The integration is performed using a given
+  quadrature formulae and assumes that the given finite element objects equals
+  that used for the computation of the solution.
+  
+  The result ist stored in a vector (named #difference#), where each entry
+  equals the given norm of the difference on one cell. The order of entries
+  is the same as a #cell_iterator# takes when started with #begin_active# and
+  promoted with the #++# operator.
+  
+  You can use the #distribute_cell_to_dof_vector# function to convert cell
+  based data to a data vector with values on the degrees of freedom, which
+  can then be attached to a #DataOut# object to be printed.
+  
+  Presently, there is the possibility to compute the following values from the
+  difference, on each cell: #mean#, #L1_norm#, #L2_norm#, #Linfty_norm#.
+  For the mean difference value, the reference function minus the numerical
+  solution is computed, not the other way round.
+
+  The infinity norm of the difference on a given cell returns the maximum
+  absolute value of the difference at the quadrature points given by the
+  quadrature formula parameter. This will in some cases not be too good
+  an approximation, since for example the Gauss quadrature formulae do
+  not evaluate the difference at the end or corner points of the cells.
+  You may want to chose a quadrature formula with more quadrature points
+  or one with another distribution of the quadrature points in this case.
+  
+  To get the {\it global} L_1 error, you have to sum up the entries in
+  #difference#, e.g. using the STL function
+  #accumulate (d.begin(), d.end(), 0)#, if #d# is the difference vector.
+  For the global L_2 difference, you have to sum up the squares of the
+  entries and take the root of the sum. These two operations represent the
+  l_1 and l_2 norms of the vectors, but you need not take the absolute
+  value of each entry, since the cellwise norms are already positive.
+  
+  To get the global mean difference, simply sum up the elements as above.
+  To get the L_\infty norm, take the maximum of the vector elements, e.g.
+  using the STL function #max_element (d.begin(), d.end())#.
   */
 template <int dim>
 class ProblemBase {
   public:
+                                    /**
+                                     * Declare a data type which denotes a
+                                     * mapping between a boundary indicator
+                                     * and the function denoting the boundary
+                                     * values on this part of the boundary.
+                                     * Only one boundary function may be given
+                                     * for each boundary indicator, which is
+                                     * guaranteed by the #map# data type.
+                                     *
+                                     * See the general documentation of this
+                                     * class for more detail.
+                                     */
+    typedef map<unsigned char,Function<dim> > DirichletBC;
+    
                                     /**
                                      * Constructor. Use this triangulation and
                                      * degree of freedom object during the
@@ -80,16 +218,20 @@ class ProblemBase {
                                      * Initiate the process of assemblage of
                                      * vectors and system matrix. Use the
                                      * given equation object and the given
-                                     * quadrature formula.
+                                     * quadrature formula. Also use the list
+                                     * of dirichlet boundary value functions
+                                     * (by default, no dirichlet bc are assumed
+                                     * which means that all bc are included
+                                     * into the weak formulation).
                                      *
                                      * For what exactly happens here, refer to
                                      * the general doc of this class.
                                      */
     virtual void assemble (const Equation<dim>      &equation,
                           const Quadrature<dim>    &q,
-                          const FiniteElement<dim> &fe);
-
-
+                          const FiniteElement<dim> &fe,
+                          const DirichletBC        &dirichlet_bc = DirichletBC());
+    
                                     /**
                                      * Solve the system of equations.
                                      */
@@ -100,41 +242,16 @@ class ProblemBase {
                                      * the solution computed before and
                                      * the reference solution, which
                                      * is given as a continuous function
-                                     * object. The integration is
-                                     * performed using the given quadrature
-                                     * rule and assumes that the given
-                                     * finite element objects equals
-                                     * that used for the computation of
-                                     * the solution.
-                                     *
-                                     * The result ist stored in the
-                                     * #difference# object, where each
-                                     * entry equals the L_1 norm of
-                                     * the difference on one cell. The
-                                     * order of entries is the same as
-                                     * a #cell_iterator# takes when
-                                     * started with #begin_active# and
-                                     * promoted with the #++# operator.
-                                     *
-                                     * You can use the
-                                     * #distribute_cell_to_dof_vector#
-                                     * function to convert cell based
-                                     * data to a data vector with values
-                                     * on the degrees of freedom, which
-                                     * can then be attached to a
-                                     * #DataOut# object.
+                                     * object.
                                      *
-                                     * To get the global L_1 error,
-                                     * you have to sum up the entries
-                                     * in #difference#, e.g. using the STL
-                                     * function
-                                     * #accumulate (d.begin(), d.end(), 0)#,
-                                     * if #d# is the difference vector.
-                                     */
-    void integrate_L1_difference (const Function<dim>      &exact_solution,
-                                 vector<double>           &difference,
-                                 const Quadrature<dim>    &q,
-                                 const FiniteElement<dim> &fe) const;
+                                     * See the general documentation of this
+                                     * class for more information.
+                                     */
+    void integrate_difference (const Function<dim>      &exact_solution,
+                              vector<double>           &difference,
+                              const Quadrature<dim>    &q,
+                              const FiniteElement<dim> &fe,
+                              const NormType           &norm) const;
     
                                     /**
                                      * Initialize the #DataOut# object with
@@ -178,6 +295,10 @@ class ProblemBase {
                                      * Exception
                                      */
     DeclException0 (ExcInvalidFE);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
     
   protected:
                                     /**
@@ -216,8 +337,19 @@ class ProblemBase {
                                      * hanging nodes.
                                      */
     ConstraintMatrix    constraints;
+
+                                    /**
+                                     * Apply dirichlet boundary conditions
+                                     * to the system matrix and vectors
+                                     * as described in the general
+                                     * documentation.
+                                     */
+    void apply_dirichlet_bc (dSMatrix &matrix,
+                            dVector  &solution,
+                            dVector  &right_hand_side,
+                            const DirichletBC &dirichlet_bc);
     
-  friend class Assembler<dim>;
+    friend class Assembler<dim>;
 };
 
     
index 6d35cb4956114736b7c3bb45c4ae7d9d5215580e..e77800bab8fefc13824f9a2238b9f7eb5b6354ad 100644 (file)
@@ -21,8 +21,6 @@ template <int dim> class DoFHandler;
 
 class dVector;
 
-template <class Key, class T, class Compare, class Alloc = alloc> class map;
-template <class T> struct less;
 
 
 
index cc905ae674cb8e2770a31ce35d75b0e6b3658ebd..17e7817ad898306867f59a617fa04bb9fdb519e1 100644 (file)
@@ -256,13 +256,29 @@ DoFCellAccessor<dim>::child (const unsigned int i) const {
 
 
 
+
+DoFSubstructAccessor<1>::substruct_iterator
+DoFCellAccessor<1>::face (const unsigned int) const {
+  Assert (false, ExcNotUsefulForThisDimension());
+  return 0;
+};
+
+
+
+DoFSubstructAccessor<2>::substruct_iterator
+DoFCellAccessor<2>::face (const unsigned int i) const {
+  return line(i);
+};
+
+
+
 void
-DoFCellAccessor<1>::dof_indices (vector<int> &dof_indices) const {
+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_indices.size() +
-                      dof_handler->get_selected_fe().total_dofs);
+  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));
@@ -274,23 +290,72 @@ DoFCellAccessor<1>::dof_indices (vector<int> &dof_indices) const {
 
 
 void
-DoFCellAccessor<2>::dof_indices (vector<int> &dof_indices) const {
+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_indices.size() +
-                      dof_handler->get_selected_fe().total_dofs);
+  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 {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == 0, ExcVectorNotEmpty());
+  Assert ((unsigned int)values.n() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+  
+  dof_values.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_values.push_back (values(vertex_dof_index(vertex,d)));
+  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
+    dof_values.push_back (values(dof_index(d)));
 };
 
 
 
 
+void
+DoFCellAccessor<2>::get_dof_values (const dVector  &values,
+                                   vector<double> &dof_values) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == 0, ExcVectorNotEmpty());
+  Assert ((unsigned int)values.n() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+    
+  dof_values.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_values.push_back (values(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_values.push_back (values(this->line(line)->dof_index(d)));
+  for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_quad; ++d)
+    dof_values.push_back (values(dof_index(d)));
+};
+
+
+
+
+
+
+
+
 
 // explicit instantiations
 template class DoFLineAccessor<1,CellAccessor<1> >;
index 8dec8ca80507a46bf8182dd86ae98c93199a590b..643ceb295c3ebfc4cfc819f47347928b647aeab9 100644 (file)
@@ -1037,8 +1037,8 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
                                       // both cells have no children
       {
        vector<int> old_dofs, new_dofs;
-       old_cell->dof_indices (old_dofs);
-       new_cell->dof_indices (new_dofs);
+       old_cell->get_dof_indices (old_dofs);
+       new_cell->get_dof_indices (new_dofs);
        Assert (old_dofs.size() == selected_fe->total_dofs,
                ExcInternalError ());
        Assert (new_dofs.size() == selected_fe->total_dofs,
@@ -1062,7 +1062,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
          
                                           // numbers of old dofs
          vector<int> old_dof_indices;
-         old_cell->dof_indices (old_dof_indices);
+         old_cell->get_dof_indices (old_dof_indices);
 
          Assert (old_dof_indices.size() == selected_fe->total_dofs,
                  ExcInternalError ());
@@ -1071,7 +1071,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
            {
                                               // numbers of child dofs
              vector<int> child_dof_indices;
-             child[c]->dof_indices (child_dof_indices);
+             child[c]->get_dof_indices (child_dof_indices);
 
              Assert (child_dof_indices.size() == selected_fe->total_dofs,
                      ExcInternalError ());
@@ -1094,7 +1094,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
          
                                                   // numbers of new dofs
          vector<int> new_dof_indices;
-         new_cell->dof_indices(new_dof_indices);
+         new_cell->get_dof_indices(new_dof_indices);
 
          Assert (new_dof_indices.size() == selected_fe->total_dofs,
                  ExcInternalError ());
@@ -1103,7 +1103,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
            {
                                               // numbers of child dofs
              vector<int> child_dof_indices;
-             child[c]->dof_indices (child_dof_indices);
+             child[c]->get_dof_indices (child_dof_indices);
 
              Assert (child_dof_indices.size() == selected_fe->total_dofs,
                      ExcInternalError ());
@@ -1133,8 +1133,8 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
                                       // both cells have no children
       {
        vector<int> old_dofs, new_dofs;
-       old_cell->dof_indices (old_dofs);
-       new_cell->dof_indices (new_dofs);
+       old_cell->get_dof_indices (old_dofs);
+       new_cell->get_dof_indices (new_dofs);
        Assert (old_dofs.size() == selected_fe->total_dofs,
                ExcInternalError ());
        Assert (new_dofs.size() == selected_fe->total_dofs,
@@ -1158,7 +1158,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
          
                                           // numbers of old dofs
          vector<int> old_dof_indices;
-         old_cell->dof_indices (old_dof_indices);
+         old_cell->get_dof_indices (old_dof_indices);
 
          Assert (old_dof_indices.size() == selected_fe->total_dofs,
                  ExcInternalError ());
@@ -1167,7 +1167,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
            {
                                               // numbers of child dofs
              vector<int> child_dof_indices;
-             child[c]->dof_indices (child_dof_indices);
+             child[c]->get_dof_indices (child_dof_indices);
 
              Assert (child_dof_indices.size() == selected_fe->total_dofs,
                      ExcInternalError ());
@@ -1191,7 +1191,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
          
                                                   // numbers of new dofs
          vector<int> new_dof_indices;
-         new_cell->dof_indices(new_dof_indices);
+         new_cell->get_dof_indices(new_dof_indices);
 
          Assert (new_dof_indices.size() == selected_fe->total_dofs,
                  ExcInternalError ());
@@ -1200,7 +1200,7 @@ void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterat
            {
                                               // numbers of child dofs
              vector<int> child_dof_indices;
-             child[c]->dof_indices (child_dof_indices);
+             child[c]->get_dof_indices (child_dof_indices);
 
              Assert (child_dof_indices.size() == selected_fe->total_dofs,
                      ExcInternalError ());
index c9df0fa51959909138c3a00f79c0973c5dcbe7ac..a6228c505c157857c04b77a8459b3ef38793b4dc 100644 (file)
@@ -529,6 +529,14 @@ void CellAccessor<1>::set_material_id (const unsigned char mat_id) const {
 
 
 
+Triangulation<1>::substruct_iterator
+CellAccessor<1>::face (const unsigned int) const {
+  Assert (false, ExcNotUsefulForThisDimension());
+  return 0;
+};
+
+
+
 /*------------------------ Functions: CellAccessor<2> -----------------------*/
 
 
@@ -554,6 +562,13 @@ void CellAccessor<2>::set_material_id (const unsigned char mat_id) const {
 };
 
 
+
+Triangulation<2>::substruct_iterator
+CellAccessor<2>::face (const unsigned int i) const {
+  return line(i);
+};
+
+
 /*------------------------ Functions: CellAccessor<dim> -----------------------*/
 
 
index 44f76cedada28d45c78e2a6d5ce4ba1eb31e0c09..218a66cddb5a8e5786d0c05d5c6598baf3d30ce8 100644 (file)
@@ -96,11 +96,12 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
                                                       present_level,
                                                       present_index),
                    fe);
+  const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs;
 
                                   // clear cell matrix
   if (assemble_matrix)
-    for (unsigned int i=0; i<dof_handler->get_selected_fe().total_dofs; ++i)
-      for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
+    for (unsigned int i=0; i<n_dofs; ++i)
+      for (unsigned int j=0; j<n_dofs; ++j)
        cell_matrix(i,j) = 0;
   
 
@@ -133,17 +134,17 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
 
                                   // get indices of dofs
   vector<int> dofs;
-  dof_indices (dofs);
+  get_dof_indices (dofs);
   
                                   // distribute cell matrix
   if (assemble_matrix)
-    for (unsigned int i=0; i<dof_handler->get_selected_fe().total_dofs; ++i)
-      for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
+    for (unsigned int i=0; i<n_dofs; ++i)
+      for (unsigned int j=0; j<n_dofs; ++j)
        problem.system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
 
                                   // distribute cell vector
   if (assemble_rhs)
-    for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
+    for (unsigned int j=0; j<n_dofs; ++j)
       problem.right_hand_side(dofs[j]) += cell_vector(j);
 };
 
index 0874fd4a506a0145da0bb1935ea3a51ef7b83c2e..08e75e5239f4cd9ef9e432c45a124e2f3e75f1db 100644 (file)
 #include "../../../mia/vectormemory.h"
 #include "../../../mia/cg.h"
 
+#include <algo.h>
+
+
+
+
+inline double sqr (double x) {
+  return x*x;
+};
 
 
 
@@ -22,8 +30,8 @@ ProblemBase<dim>::ProblemBase (Triangulation<dim> *tria,
                               DoFHandler<dim>    *dof) :
                tria(tria),
                dof_handler(dof),
-               system_sparsity(1,1,1),        // dummy initialisation, is later reinit'd
-               system_matrix()               // dummy initialisation, is later reinit'd
+               system_sparsity(1,1,1),   // dummy initialisation, is later reinit'd
+               system_matrix()           // dummy initialisation, is later reinit'd
 {
   Assert (tria == &dof->get_tria(), ExcDofAndTriaDontMatch());
 };
@@ -36,9 +44,10 @@ ProblemBase<dim>::~ProblemBase () {};
 
 
 template <int dim>
-void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
-                                const Quadrature<dim> &quadrature,
-                                const FiniteElement<dim> &fe) {
+void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
+                                const Quadrature<dim>    &quadrature,
+                                const FiniteElement<dim> &fe,
+                                const DirichletBC        &dirichlet_bc) {
   system_sparsity.reinit (dof_handler->n_dofs(),
                          dof_handler->n_dofs(),
                          dof_handler->max_couplings_between_dofs());
@@ -53,6 +62,8 @@ void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
 
                                   // reinite system matrix
   system_matrix.reinit (system_sparsity);
+                                  // reinit solution vector, preset
+                                  // with zeroes.
   solution.reinit (dof_handler->n_dofs());
 
                                   // create assembler
@@ -73,6 +84,11 @@ void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
     }
   while ((++assembler).state() == valid);
 
+                                  // apply Dirichlet bc as described
+                                  // in the docs
+  apply_dirichlet_bc (system_matrix, solution,
+                     right_hand_side, dirichlet_bc);
+  
                                   // condense system matrix in-place
   constraints.condense (system_matrix);
 
@@ -100,28 +116,135 @@ void ProblemBase<dim>::solve () {
 
 
 
-/*
 template <int dim>
-void ProblemBase<dim>::integrate_L1_difference (const Function<dim>      &exact_solution,
-                                               vector<double>           &difference,
-                                               const Quadrature<dim>    &q,
-                                               const FiniteElement<dim> &fe) const {
+void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_solution,
+                                            vector<double>           &difference,
+                                            const Quadrature<dim>    &q,
+                                            const FiniteElement<dim> &fe,
+                                            const NormType           &norm) const {
   Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
 
   difference.erase (difference.begin(), difference.end());
   difference.reserve (tria->n_cells());
 
+  FEValues<dim> fe_values(fe, q);
+  
                                   // loop over all cells
-  DoFHandler<dim>::active_dof_iterator cell = dof_handler->begin_active(),
-                                      endc = dof_handler->end();
-  for (; cell != end; ++cell) 
+                                  // (we need an iterator on the triangulation for
+                                  // fe_values.reinit, but we also need an iterator
+                                  // on the dofhandler; we could generate the first
+                                  // out of the second or vice versa each time
+                                  // we need them, but conversion seems to be slow.
+                                  // therefore, we keep two iterators, hoping that
+                                  // incrementing is faster than conversion...)
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active(),
+                                       endc = dof_handler->end();
+  Triangulation<dim>::active_cell_iterator tria_cell = dof_handler->get_tria().begin();
+  
+  for (; cell != endc; ++cell, ++tria_cell) 
     {
       double diff=0;
+                                      // initialize for this cell
+      fe_values.reinit (tria_cell, fe);
+
+      switch (norm) 
+       {
+         case mean:
+         case L1_norm:
+         case L2_norm:
+         case Linfty_norm:
+         {
+                                            // we need the finite element
+                                            // function \psi at the different
+                                            // integration points. Compute
+                                            // it like this:
+                                            // \psi(x_j)=\sum_i v_i \phi_i(x_j)
+                                            // with v_i the nodal values of the
+                                            // solution and \phi_i(x_j) the
+                                            // matrix of the ansatz function
+                                            // values at the integration point
+                                            // x_j. Then the vector
+                                            // of the \psi(x_j) is v*Phi with
+                                            // v being the vector of nodal
+                                            // values on this cell and Phi
+                                            // the matrix.
+                                            //
+                                            // we then need the difference:
+                                            // reference_function(x_j)-\psi_j
+                                            // and assign that to the vector
+                                            // \psi.
+           const unsigned int n_dofs = fe.total_dofs;
+           const dFMatrix & shape_values = fe_values.get_shape_values();
+           vector<double>   dof_values;
+           cell->get_dof_values (solution, dof_values);
+
+           vector<double>   psi;
+
+                                            // in praxi: first compute
+                                            // exact solution vector
+           exact_solution.value_list (fe_values.get_quadrature_points(),
+                                      psi);
+                                            // then subtract finite element
+                                            // solution
+           for (unsigned int j=0; j<n_dofs; ++j)
+             for (unsigned int i=0; i<n_dofs; ++i)
+               psi[j] -= dof_values[i]*shape_values(i,j);
+
+                                            // for L1_norm and Linfty_norm:
+                                            // take absolute
+                                            // value, for the L2_norm take
+                                            // square of psi
+           switch (norm) 
+             {
+               case mean:
+                     break;
+               case L1_norm:
+               case Linfty_norm:
+                     transform (psi.begin(), psi.end(),
+                                psi.begin(), ptr_fun(abs));
+                     break;
+               case L2_norm:
+                     transform (psi.begin(), psi.end(),
+                                psi.begin(), ptr_fun(sqr));
+                     break;
+             };
+
+
+                                            // ok, now we have the integrand,
+                                            // let's compute the integral,
+                                            // which is
+                                            // sum_j psi_j JxW_j
+                                            // (or |psi_j| or |psi_j|^2
+           switch (norm) 
+             {
+               case mean:
+               case L1_norm:
+               case L2_norm:
+                     diff = inner_product (psi.begin(), psi.end(),
+                                           fe_values.get_JxW_values().begin(),
+                                           0.0);
+                     break;
+               case Linfty_norm:
+                     diff = *max_element (psi.begin(), psi.end());
+                     break;
+             };
+
+           if (norm==L2_norm)
+             diff = sqrt(diff);
+           
+           break;
+         };
+
+         default:
+               Assert (false, ExcNotImplemented());
+       };
+
       
+                                      // append result of this cell
+                                      // to the end of the vector
       difference.push_back (diff);
     };
 };
-*/
 
 
 
@@ -136,14 +259,23 @@ void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
 };
 
 
-/*
+
 template <int dim>
 pair<char*,char*> ProblemBase<dim>::get_solution_name () const {
-  return make_pair("solution", "");
+  return pair<char*,char*>("solution", "");
+};
+
+
+
+
+template <int dim>
+void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
+                                          dVector  &solution,
+                                          dVector  &right_hand_side,
+                                          const DirichletBC &dirichlet_bc) {
+//  l;
 };
-*/
 
-  
   
 
 

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.