]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some functions from ProblemBase to more useful locations.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 19 May 1998 13:51:14 +0000 (13:51 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 19 May 1998 13:51:14 +0000 (13:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@315 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc

index f914af98658060fb5b5839e3eca09cfa37a033aa..82e76641d08a7a38e19f70bcc0867bc2d571ba52 100644 (file)
@@ -25,27 +25,6 @@ template <int dim> class StraightBoundary;
 template <int dim> class Function;
 
 
-/**
- *  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,
-      H1_seminorm,
-      H1_norm
-};
-
-
 
 
 /**
@@ -102,136 +81,11 @@ enum NormType {
  * 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 these 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 actually apply the boundary conditions, use is made of the
+ * #MatrixTools::apply_boundary_values# function and by interpolation of
+ * the boundary_values using the #MatrixTool::interpolate_boundary_values#
+ * function. See there for more information.
  *
- * 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
- * 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, 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
- * 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.
- *
- * 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.
- *
- * You should be aware that the boundary function may be evaluated at nodes
- * on the interior of faces. These, however, need not be on the true
- * boundary, but rather are on the approximation of the boundary represented
- * by teh mapping of the unit cell to the real cell. Since this mapping will
- * in most cases not be the exact one at the face, the boundary function is
- * evaluated at points which are not on the boundary and you should make
- * sure that the returned values are reasonable in some sense anyway.
- * 
- *
- * \subsection{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 of the #DoFHandler#
- * class 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#,
- * #H1_seminorm#.
- * 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.
- * You should also take into account the superconvergence properties of finite
- * elements in some points: for example in 1D, the standard finite element
- * method is a collocation method and should return the exact value at nodal
- * points. Therefore, the trapezoidal rule should always return a vanishing
- * L-infinity error. Conversely, in 2D the maximum L-infinity error should
- * be located at the vertices or at the center of the cell, which would make
- * it plausible to use the Simpson quadrature rule. On the other hand, there
- * may be superconvergence at Gauss integration points. These examples are not
- * intended as a rule of thumb, rather they are though to illustrate that the
- * use of the wrong quadrature formula may show a significantly wrong result
- * and care should be taken to chose the right formula.
- *
- * The $H_1$ seminorm is the $L_2$ norm of the gradient of the difference. The
- * full $H_1$ norm is the sum of the seminorm and the $L_2$ norm.
- * 
- * To get the {\it global} L_1 error, you have to sum up the entries in
- * #difference#, e.g. using #dVector::l1_norm# function.
- * For the global L_2 difference, you have to sum up the squares of the
- * entries and take the root of the sum, e.g. using #dVector::l2_norm.
- * 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 #dVector::linfty_norm# function.
- *
- * For the global $H_1$ norm and seminorm, the same rule applies as for the
- * $L_2$ norm: compute the $l_2$ norm of the cell error vector.
  */
 template <int dim>
 class ProblemBase {
@@ -311,27 +165,11 @@ class ProblemBase {
                           const Boundary<dim>      &boundary = StraightBoundary<dim>());
     
                                     /**
-                                     * Solve the system of equations.
+                                     * Solve the system of equations. This uses
+                                     * a simple CG method.
                                      */
     virtual void solve ();
 
-                                    /**
-                                     * Integrate the difference between
-                                     * the solution computed before and
-                                     * the reference solution, which
-                                     * is given as a continuous function
-                                     * object.
-                                     *
-                                     * See the general documentation of this
-                                     * class for more information.
-                                     */
-    void integrate_difference (const Function<dim>      &exact_solution,
-                              dVector                  &difference,
-                              const Quadrature<dim>    &q,
-                              const FiniteElement<dim> &fe,
-                              const NormType           &norm,
-                              const Boundary<dim> &boundary=StraightBoundary<dim>()) const;
-
                                     /**
                                      * Initialize the #DataOut# object with
                                      * the grid and DoF handler used for this
@@ -362,20 +200,6 @@ class ProblemBase {
                                      */
     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 FunctionMap        &dirichlet_bc,
-                                          const FiniteElement<dim> &fe,
-                                          const Boundary<dim>      &boundary,
-                                          map<int,double>          &boundary_values) const;
-    
                                     /**
                                      * Exception
                                      */
@@ -383,22 +207,6 @@ class ProblemBase {
                                     /**
                                      * Exception
                                      */
-    DeclException0 (ExcNoMemory);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInvalidFE);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotImplemented);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInvalidBoundaryIndicator);
-                                    /**
-                                     * Exception
-                                     */
     DeclException0 (ExcNoTriaSelected);
     
   protected:
@@ -439,19 +247,6 @@ class ProblemBase {
                                      */
     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 FunctionMap &dirichlet_bc,
-                            const FiniteElement<dim> &fe,
-                            const Boundary<dim>      &boundary);
-    
     friend class Assembler<dim>;
 };
 
index de9b94d9db613e076773853e9798e6f03fddf44b..501f29cd63211bdeebeef815bb854a2ab6c1c53b 100644 (file)
@@ -7,6 +7,8 @@
 
 
 #include <base/exceptions.h>
+#include <map>
+
 
 template <int dim> class Triangulation;
 template <int dim> class DoFHandler;
@@ -171,6 +173,161 @@ class MatrixCreator {
 
 
 
+/**
+ * Provide a collection of functions operating on matrices. These include
+ * the application of boundary conditions to a linear system of equations
+ * and others.
+ *
+ *
+ * \subsection{Boundary conditions}
+ *
+ * The #apply_boundar_values# function inserts boundary conditions of
+ * into a system of equations.  To actually do this you have to specify
+ * a list of degree of freedom indices along with the value this degree of
+ * freedom shall assume. To see how to get such a list, see below in the
+ * discussion of the #interpolate_boundary_values# function.
+ *
+ * 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.
+ *
+ * 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
+ * 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, 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
+ * 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.
+ *
+ *
+ * \subsection{Getting a list of boundary values}
+ *
+ * As discussed above, the #apply_boundary_values# function takes a list
+ * of boundary nodes and their values. You can get such a list by interpolation
+ * of a boundary function using the #interpolate_boundary_values# function.
+ * To use it, you have to
+ * 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 these into the process of assemblage therefore.
+ *
+ * Within this function, 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. For
+ * this purpose other functions exist in the #MatrixTools# library (or will
+ * exist at least).
+ *
+ * You should be aware that the boundary function may be evaluated at nodes
+ * on the interior of faces. These, however, need not be on the true
+ * boundary, but rather are on the approximation of the boundary represented
+ * by teh mapping of the unit cell to the real cell. Since this mapping will
+ * in most cases not be the exact one at the face, the boundary function is
+ * evaluated at points which are not on the boundary and you should make
+ * sure that the returned values are reasonable in some sense anyway.
+ * 
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MatrixTools : public MatrixCreator<dim> {
+  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>*> FunctionMap;
+
+                                    /**
+                                     * Apply dirichlet boundary conditions
+                                     * to the system matrix and vectors
+                                     * as described in the general
+                                     * documentation.
+                                     */
+    static void apply_boundary_values (const map<int,double> &boundary_values,
+                                      dSMatrix              &matrix,
+                                      dVector               &solution,
+                                      dVector               &right_hand_side);
+
+                                    /**
+                                     * Make up the list of node subject
+                                     * to Dirichlet boundary conditions
+                                     * and the values they are to be
+                                     * assigned, by interpolation around
+                                     * the boundary.
+                                     *
+                                     * See the general doc for more
+                                     * information.
+                                     */
+    static void interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                            const FunctionMap     &dirichlet_bc,
+                                            const FiniteElement<dim> &fe,
+                                            const Boundary<dim> &boundary,
+                                            map<int,double>     &boundary_values);
+    
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidBoundaryIndicator);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
+};
+
+
+
+
 /**
  * Equation class to be passed to the #Assembler# if you want to make up the
  * mass matrix for your problem. The mass matrix is the matrix with
@@ -377,6 +534,9 @@ class LaplaceMatrix :  public Equation<dim> {
 
 
 
+
+
+
 /*----------------------------   matrices.h     ---------------------------*/
 /* end of #ifndef __matrices_H */
 #endif
index c821ddd9d257827189430144ab0a0647ff4918d3..12641af42bbf3c333b73e12d9b0ac0ee9f4b956b 100644 (file)
@@ -5,18 +5,46 @@
 /*----------------------------   vectors.h     ---------------------------*/
 
 
+
+#include <base/exceptions.h>
+
 template <int dim> class DoFHandler;
 template <int dim> class Function;
 template <int dim> class Quadrature;
-template <int dim> class Boundary;
 template <int dim> class FiniteElement;
+template <int dim> class Boundary;
+template <int dim> class StraightBoundary;
 class ConstraintMatrix;
 class dVector;
 
 
+/**
+ *  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,
+      H1_seminorm,
+      H1_norm
+};
+
+
+
 
 /**
- * Provide a class which assembles some standard vectors. Among these are
+ * Provide a class which offers some operations on vectors. Amoung these are
+ * assemblage of standard vectors, integration of the difference of a
+ * finite element solution and a continuous function,
  * interpolations and projections of continuous functions to the finite
  * element space and other operations.
  *
@@ -57,10 +85,69 @@ class dVector;
  *   method without preconditioning and without multigrid. This is clearly not
  *   too efficient, but sufficient in many cases and simple to implement. This
  *   detail may change in the future.
+ *
+ * \item 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 of the #DoFHandler#
+ *   class 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#,
+ *   #H1_seminorm#.
+ *   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.
+ *   You should also take into account the superconvergence properties of finite
+ *   elements in some points: for example in 1D, the standard finite element
+ *   method is a collocation method and should return the exact value at nodal
+ *   points. Therefore, the trapezoidal rule should always return a vanishing
+ *   L-infinity error. Conversely, in 2D the maximum L-infinity error should
+ *   be located at the vertices or at the center of the cell, which would make
+ *   it plausible to use the Simpson quadrature rule. On the other hand, there
+ *   may be superconvergence at Gauss integration points. These examples are not
+ *   intended as a rule of thumb, rather they are though to illustrate that the
+ *   use of the wrong quadrature formula may show a significantly wrong result
+ *   a nd care should be taken to chose the right formula.
+ *
+ *   The $H_1$ seminorm is the $L_2$ norm of the gradient of the difference. The
+ *   full $H_1$ norm is the sum of the seminorm and the $L_2$ norm.
+ * 
+ *   To get the {\it global} L_1 error, you have to sum up the entries in
+ *   #difference#, e.g. using #dVector::l1_norm# function.
+ *   For the global L_2 difference, you have to sum up the squares of the
+ *   entries and take the root of the sum, e.g. using #dVector::l2_norm.
+ *   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 #dVector::linfty_norm# function.
+ *
+ *   For the global $H_1$ norm and seminorm, the same rule applies as for the
+ *   $L_2$ norm: compute the $l_2$ norm of the cell error vector.
  * \end{itemize}
  */
 template <int dim>
-class VectorCreator {
+class VectorTools {
   public:
                                     /**
                                      * Compute the interpolation of
@@ -90,6 +177,34 @@ class VectorCreator {
                         const Boundary<dim>      &boundary,
                         const Function<dim>      &function,
                         dVector                  &vec);
+
+                                    /**
+                                     * Integrate the difference between
+                                     * a finite element function and
+                                     * the reference function, which
+                                     * is given as a continuous function
+                                     * object.
+                                     *
+                                     * See the general documentation of this
+                                     * class for more information.
+                                     */
+    static void integrate_difference (const DoFHandler<dim>    &dof,
+                                     const dVector            &fe_function,
+                                     const Function<dim>      &exact_solution,
+                                     dVector                  &difference,
+                                     const Quadrature<dim>    &q,
+                                     const FiniteElement<dim> &fe,
+                                     const NormType           &norm,
+                                     const Boundary<dim> &boundary=StraightBoundary<dim>());
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidFE);
 };
 
 
index 1a0618c8e8d72ffdf180faeac4b0665101f2ef05..ebb649a0d0039f8412677c369aaf7e874971f997 100644 (file)
@@ -2,6 +2,7 @@
 
 #include <numerics/assembler.h>
 #include <numerics/base.h>
+#include <numerics/matrices.h>
 #include <grid/dof_constraints.h>
 #include <grid/tria_iterator.h>
 #include <basic/data_io.h>
 
 
 
-inline double sqr (const double x) {
-  return x*x;
-};
-
-
-template <int dim>
-inline double sqr_point (const Point<dim> &p) {
-  return p.square();
-};
-
-
-
 
 template <int dim>
 ProblemBase<dim>::ProblemBase () :
@@ -128,9 +117,13 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
 
                                   // apply Dirichlet bc as described
                                   // in the docs
-  apply_dirichlet_bc (system_matrix, solution,
-                     right_hand_side,
-                     dirichlet_bc, fe, boundary);
+  map<int, double> boundary_value_list;
+  MatrixTools<dim>::interpolate_boundary_values (*dof_handler,
+                                                dirichlet_bc, fe, boundary,
+                                                boundary_value_list);
+  MatrixTools<dim>::apply_boundary_values (boundary_value_list,
+                                          system_matrix, solution,
+                                          right_hand_side);
   
 };
 
@@ -156,196 +149,6 @@ void ProblemBase<dim>::solve () {
 
 
 
-template <int dim>
-void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_solution,
-                                            dVector                  &difference,
-                                            const Quadrature<dim>    &q,
-                                            const FiniteElement<dim> &fe,
-                                            const NormType           &norm,
-                                            const Boundary<dim>      &boundary) const {
-  Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());  
-  Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
-
-  difference.reinit (tria->n_active_cells());
-  
-  UpdateFlags update_flags = UpdateFlags (update_q_points  |
-                                         update_jacobians |
-                                         update_JxW_values);
-  if ((norm==H1_seminorm) || (norm==H1_norm))
-    update_flags = UpdateFlags (update_flags | update_gradients);
-  FEValues<dim> fe_values(fe, q, update_flags);
-  
-                                  // loop over all cells
-  DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active(),
-                                       endc = dof_handler->end();
-  for (unsigned int index=0; cell != endc; ++cell, ++index)
-    {
-      double diff=0;
-                                      // initialize for this cell
-      fe_values.reinit (cell, fe, boundary);
-
-      switch (norm) 
-       {
-         case mean:
-         case L1_norm:
-         case L2_norm:
-         case Linfty_norm:
-         case H1_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_q_points = q.n_quadrature_points;
-           vector<double>   psi (n_q_points);
-
-                                            // in praxi: first compute
-                                            // exact solution vector
-           exact_solution.value_list (fe_values.get_quadrature_points(),
-                                      psi);
-                                            // then subtract finite element
-                                            // solution
-           if (true) 
-             {
-               vector<double> function_values (n_q_points, 0);
-               fe_values.get_function_values (solution, function_values);
-
-               transform (psi.begin(), psi.end(),
-                          function_values.begin(),
-                          psi.begin(),
-                          minus<double>());
-             };
-           
-
-                                            // 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(fabs));
-                     break;
-               case L2_norm:
-               case H1_norm:
-                     transform (psi.begin(), psi.end(),
-                                psi.begin(), ptr_fun(sqr));
-                     break;
-               default:
-                     Assert (false, ExcNotImplemented());
-             };
-
-                                            // 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:
-                     diff = inner_product (psi.begin(), psi.end(),
-                                           fe_values.get_JxW_values().begin(),
-                                           0.0);
-                     break;
-               case L2_norm:
-               case H1_norm:
-                     diff = sqrt(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;
-               default:
-                     Assert (false, ExcNotImplemented());
-             };
-
-                                            // note: the H1_norm uses the result
-                                            // of the L2_norm and control goes
-                                            // over to the next case statement!
-           if (norm != H1_norm)
-             break;
-         };
-
-         case H1_seminorm:
-         {
-                                            // note: the computation of the
-                                            // H1_norm starts at the previous
-                                            // case statement, but continues
-                                            // here!
-
-                                            // for H1_norm: re-square L2_norm.
-           diff = sqr(diff);
-
-                                            // same procedure as above, but now
-                                            // psi is a vector of gradients
-           const unsigned int n_q_points = q.n_quadrature_points;
-           vector<Point<dim> >   psi (n_q_points);
-
-                                            // in praxi: first compute
-                                            // exact solution vector
-           exact_solution.gradient_list (fe_values.get_quadrature_points(),
-                                         psi);
-           
-                                            // then subtract finite element
-                                            // solution
-           if (true) 
-             {
-               vector<Point<dim> > function_grads (n_q_points, Point<dim>());
-               fe_values.get_function_grads (solution, function_grads);
-
-               transform (psi.begin(), psi.end(),
-                          function_grads.begin(),
-                          psi.begin(),
-                          minus<Point<dim> >());
-             };
-                                            // take square of integrand
-           vector<double> psi_square (psi.size(), 0.0);
-           for (unsigned int i=0; i<n_q_points; ++i)
-             psi_square[i] = sqr_point(psi[i]);
-
-                                            // add seminorm to L_2 norm or
-                                            // to zero
-           diff += inner_product (psi_square.begin(), psi_square.end(),
-                                  fe_values.get_JxW_values().begin(),
-                                  0.0);
-           diff = sqrt(diff);
-
-           break;
-         };
-                                            
-         default:
-               Assert (false, ExcNotImplemented());
-       };
-
-      
-                                      // append result of this cell
-                                      // to the end of the vector
-      difference(index) = diff;
-    };
-};
-
-
 
 template <int dim>
 void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
@@ -369,163 +172,6 @@ pair<char*,char*> ProblemBase<dim>::get_solution_name () const {
 
 
 
-template <int dim>
-void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
-                                          dVector  &solution,
-                                          dVector  &right_hand_side,
-                                          const FunctionMap &dirichlet_bc,
-                                          const FiniteElement<dim> &fe,
-                                          const Boundary<dim>      &boundary) {
-  Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
-  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, fe, boundary, boundary_values);
-
-  map<int,double>::const_iterator dof, endd;
-  const unsigned int n_dofs   = matrix.m();
-  const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern();
-  const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
-  const int          *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 (unsigned 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. Normally, however,
-                                      // the main diagonal entry should
-                                      // not be zero.
-                                      //
-                                      // store the new rhs entry to make
-                                      // the gauss step more efficient
-      double new_rhs;
-      if (matrix.diag_element(dof->first) != 0.0)
-       new_rhs = 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);
-         new_rhs = 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 (unsigned int j=sparsity_rowstart[row];
-            j<sparsity_rowstart[row+1]; ++j)
-         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 *
-                                     new_rhs;
-             
-                                              // 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 FunctionMap &,
-                                         const FiniteElement<1> &,
-                                         const Boundary<1> &,
-                                         map<int,double>   &) const {
-  Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());  
-  Assert (false, ExcNotImplemented());
-};
-
-
-
-
-template <int dim>
-void
-ProblemBase<dim>::make_boundary_value_list (const FunctionMap        &dirichlet_bc,
-                                           const FiniteElement<dim> &fe,
-                                           const Boundary<dim>      &boundary,
-                                           map<int,double>   &boundary_values) const {
-  Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
-
-                                  // use two face iterators, since we need
-                                  // a DoF-iterator for the dof indices, but
-                                  // a Tria-iterator for the fe object
-  DoFHandler<dim>::active_face_iterator face = dof_handler->begin_active_face(),
-                                       endf = dof_handler->end_face();
-  
-  FunctionMap::const_iterator function_ptr;
-
-                                  // field to store the indices of dofs
-  vector<int>         face_dofs (fe.dofs_per_face);
-  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
-  vector<double>      dof_values (fe.dofs_per_face);
-       
-  for (; face!=endf; ++face)
-    if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
-       dirichlet_bc.end()) 
-                                      // face is subject to one of the
-                                      // bc listed in #dirichlet_bc#
-      {
-                                        // get indices, physical location and
-                                        // boundary values of dofs on this
-                                        // face
-       face->get_dof_indices (face_dofs);
-       fe.get_face_ansatz_points (face, boundary, dof_locations);
-       function_ptr->second->value_list (dof_locations, dof_values);
-
-                                        // enter into list
-       for (unsigned int i=0; i<face_dofs.size(); ++i)
-         boundary_values[face_dofs[i]] = dof_values[i];
-      };
-};
-
   
 
 
index cad018acfe5128bda59f4795ed2b4a134ef6c62d..2ba113cd36b12d9641a327a43daf4e4c68b64e55 100644 (file)
@@ -5,6 +5,7 @@
 #include <grid/dof_accessor.h>
 #include <grid/tria_iterator.h>
 #include <fe/quadrature.h>
+#include <fe/fe.h>
 #include <fe/fe_values.h>
 #include <numerics/matrices.h>
 #include <numerics/assembler.h>
@@ -145,6 +146,156 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
+template <int dim>
+void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_values,
+                                             dSMatrix  &matrix,
+                                             dVector   &solution,
+                                             dVector   &right_hand_side) {
+
+  map<int,double>::const_iterator dof  = boundary_values.begin(),
+                                 endd = boundary_values.end();
+  const unsigned int n_dofs             = matrix.m();
+  const dSMatrixStruct &sparsity        = matrix.get_sparsity_pattern();
+  const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
+  const int          *sparsity_colnums  = sparsity.get_column_numbers();
+
+  for (; dof != endd; ++dof)
+    {
+                                      // for each boundary dof:
+
+                                      // set entries of this line
+                                      // to zero
+      for (unsigned 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. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      double new_rhs;
+      if (matrix.diag_element(dof->first) != 0.0)
+       new_rhs = 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);
+         new_rhs = 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 (unsigned int j=sparsity_rowstart[row];
+            j<sparsity_rowstart[row+1]; ++j)
+         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 *
+                                     new_rhs;
+             
+                                              // set matrix entry to zero
+             matrix.global_entry(j) = 0.;
+           };
+      
+      
+                                      // preset solution vector
+      solution(dof->first) = dof->second;
+    };
+};
+
+
+
+
+
+void
+MatrixTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
+                                            const FunctionMap &,
+                                            const FiniteElement<1> &,
+                                            const Boundary<1> &,
+                                            map<int,double>   &) {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
+template <int dim>
+void
+MatrixTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                              const FunctionMap     &dirichlet_bc,
+                                              const FiniteElement<dim> &fe,
+                                              const Boundary<dim>      &boundary,
+                                              map<int,double>   &boundary_values) {
+  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
+         ExcInvalidBoundaryIndicator());
+                                  // use two face iterators, since we need
+                                  // a DoF-iterator for the dof indices, but
+                                  // a Tria-iterator for the fe object
+  DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+                                       endf = dof.end_face();
+  
+  FunctionMap::const_iterator function_ptr;
+
+                                  // field to store the indices of dofs
+  vector<int>         face_dofs (fe.dofs_per_face);
+  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+  vector<double>      dof_values (fe.dofs_per_face);
+       
+  for (; face!=endf; ++face)
+    if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
+       dirichlet_bc.end()) 
+                                      // face is subject to one of the
+                                      // bc listed in #dirichlet_bc#
+      {
+                                        // get indices, physical location and
+                                        // boundary values of dofs on this
+                                        // face
+       face->get_dof_indices (face_dofs);
+       fe.get_face_ansatz_points (face, boundary, dof_locations);
+       function_ptr->second->value_list (dof_locations, dof_values);
+
+                                        // enter into list
+       for (unsigned int i=0; i<face_dofs.size(); ++i)
+         boundary_values[face_dofs[i]] = dof_values[i];
+      };
+};
+
+
+
+
+
+
 template <int dim>
 MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
                             const Function<dim> * const a) :
@@ -367,6 +518,8 @@ void LaplaceMatrix<dim>::assemble (dVector             &rhs,
 
 template class MatrixCreator<1>;
 template class MatrixCreator<2>;
+template class MatrixTools<1>;
+template class MatrixTools<2>;
 template class MassMatrix<1>;
 template class MassMatrix<2>;
 template class LaplaceMatrix<1>;
index 061aa87c067481d97869f6618f41861c58c9163d..a9136e875acd53c3a366b6c887db32bbcede0a47 100644 (file)
@@ -7,6 +7,8 @@
 #include <grid/tria_iterator.h>
 #include <grid/dof_constraints.h>
 #include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/quadrature.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <lac/dvector.h>
 #include "../../../mia/vectormemory.h"
 #include "../../../mia/cg.h"
 
+#include <numeric>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+inline double sqr (const double x) {
+  return x*x;
+};
+
+
+template <int dim>
+inline double sqr_point (const Point<dim> &p) {
+  return p.square();
+};
+
+
+
+
+
 
 template <int dim>
-void VectorCreator<dim>::interpolate (const DoFHandler<dim>    &dof,
+void VectorTools<dim>::interpolate (const DoFHandler<dim>    &dof,
                                      const FiniteElement<dim> &fe,
                                      const Boundary<dim>      &boundary,
                                      const Function<dim>      &function,
@@ -49,7 +72,7 @@ void VectorCreator<dim>::interpolate (const DoFHandler<dim>    &dof,
 
 
 template <int dim>
-void VectorCreator<dim>::project (const DoFHandler<dim>    &dof,
+void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
                                  const ConstraintMatrix   &constraints,
                                  const FiniteElement<dim> &fe,
                                  const Quadrature<dim>    &q,
@@ -87,5 +110,195 @@ void VectorCreator<dim>::project (const DoFHandler<dim>    &dof,
 
 
 
-template VectorCreator<1>;
-template VectorCreator<2>;
+template <int dim>
+void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
+                                            const dVector            &fe_function,
+                                            const Function<dim>      &exact_solution,
+                                            dVector                  &difference,
+                                            const Quadrature<dim>    &q,
+                                            const FiniteElement<dim> &fe,
+                                            const NormType           &norm,
+                                            const Boundary<dim>      &boundary) {
+  Assert (fe == dof.get_selected_fe(), ExcInvalidFE());
+  
+  difference.reinit (dof.get_tria().n_active_cells());
+  
+  UpdateFlags update_flags = UpdateFlags (update_q_points  |
+                                         update_jacobians |
+                                         update_JxW_values);
+  if ((norm==H1_seminorm) || (norm==H1_norm))
+    update_flags = UpdateFlags (update_flags | update_gradients);
+  FEValues<dim> fe_values(fe, q, update_flags);
+  
+                                  // loop over all cells
+  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+                                       endc = dof.end();
+  for (unsigned int index=0; cell != endc; ++cell, ++index)
+    {
+      double diff=0;
+                                      // initialize for this cell
+      fe_values.reinit (cell, fe, boundary);
+
+      switch (norm) 
+       {
+         case mean:
+         case L1_norm:
+         case L2_norm:
+         case Linfty_norm:
+         case H1_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
+                                            // fe_function 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_q_points = q.n_quadrature_points;
+           vector<double>   psi (n_q_points);
+
+                                            // in praxi: first compute
+                                            // exact fe_function vector
+           exact_solution.value_list (fe_values.get_quadrature_points(),
+                                      psi);
+                                            // then subtract finite element
+                                            // fe_function
+           if (true) 
+             {
+               vector<double> function_values (n_q_points, 0);
+               fe_values.get_function_values (fe_function, function_values);
+
+               transform (psi.begin(), psi.end(),
+                          function_values.begin(),
+                          psi.begin(),
+                          minus<double>());
+             };            
+
+                                            // 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(fabs));
+                     break;
+               case L2_norm:
+               case H1_norm:
+                     transform (psi.begin(), psi.end(),
+                                psi.begin(), ptr_fun(sqr));
+                     break;
+               default:
+                     Assert (false, ExcNotImplemented());
+             };
+
+                                            // 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:
+                     diff = inner_product (psi.begin(), psi.end(),
+                                           fe_values.get_JxW_values().begin(),
+                                           0.0);
+                     break;
+               case L2_norm:
+               case H1_norm:
+                     diff = sqrt(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;
+               default:
+                     Assert (false, ExcNotImplemented());
+             };
+
+                                            // note: the H1_norm uses the result
+                                            // of the L2_norm and control goes
+                                            // over to the next case statement!
+           if (norm != H1_norm)
+             break;
+         };
+
+         case H1_seminorm:
+         {
+                                            // note: the computation of the
+                                            // H1_norm starts at the previous
+                                            // case statement, but continues
+                                            // here!
+
+                                            // for H1_norm: re-square L2_norm.
+           diff = sqr(diff);
+
+                                            // same procedure as above, but now
+                                            // psi is a vector of gradients
+           const unsigned int n_q_points = q.n_quadrature_points;
+           vector<Point<dim> >   psi (n_q_points);
+
+                                            // in praxi: first compute
+                                            // exact fe_function vector
+           exact_solution.gradient_list (fe_values.get_quadrature_points(),
+                                         psi);
+           
+                                            // then subtract finite element
+                                            // fe_function
+           if (true) 
+             {
+               vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+               fe_values.get_function_grads (fe_function, function_grads);
+
+               transform (psi.begin(), psi.end(),
+                          function_grads.begin(),
+                          psi.begin(),
+                          minus<Point<dim> >());
+             };
+                                            // take square of integrand
+           vector<double> psi_square (psi.size(), 0.0);
+           for (unsigned int i=0; i<n_q_points; ++i)
+             psi_square[i] = sqr_point(psi[i]);
+
+                                            // add seminorm to L_2 norm or
+                                            // to zero
+           diff += inner_product (psi_square.begin(), psi_square.end(),
+                                  fe_values.get_JxW_values().begin(),
+                                  0.0);
+           diff = sqrt(diff);
+
+           break;
+         };
+                                            
+         default:
+               Assert (false, ExcNotImplemented());
+       };
+
+      
+                                      // append result of this cell
+                                      // to the end of the vector
+      difference(index) = diff;
+    };
+};
+
+
+template VectorTools<1>;
+template VectorTools<2>;

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.