]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Purge more classes: Assembler, MassMatrix, LaplaceMatrix, Equation. big-tests do...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 27 Mar 2001 15:23:23 +0000 (15:23 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 27 Mar 2001 15:23:23 +0000 (15:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@4309 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/assembler.h [deleted file]
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/assembler.cc [deleted file]
deal.II/deal.II/source/numerics/equation.cc [deleted file]
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc
deal.II/doc/news/2001/c-3-1.html
tests/big-tests/poisson/poisson.h
tests/big-tests/problem_base.h

diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h
deleted file mode 100644 (file)
index 4c58e08..0000000
+++ /dev/null
@@ -1,311 +0,0 @@
-//----------------------------  assembler.h  ---------------------------
-//    Version: $Name$
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------  assembler.h  ---------------------------
-#ifndef __deal2__assembler_h
-#define __deal2__assembler_h
-
-
-/*----------------------------   problem_assembler.h     ---------------------------*/
-
-#include <base/exceptions.h>
-#include <lac/full_matrix.h>
-#include <lac/vector.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <fe/fe_values.h>
-#include <vector>
-
-
-/**
- * The use of this class is now deprecated!
- *
- * This is the base class for equation objects. Equations objects describe the
- * finite element discretisation of one or more equations.
- *
- * Equation objects need only provide functions which set up the cell
- * matrices and the cell right hand side. These are then automatically inserted
- * into the global matrices and vectors.
- *
- * @author Wolfgang Bangerth, 1998
- */
-template <int dim>
-class Equation
-{
-  public:
-                                    /**
-                                     * Constructor. You have to pass the number
-                                     * of equations you want to discretize, which
-                                     * equals the number of solution functions.
-                                     */
-    Equation (const unsigned int n_equations);
-
-                                    /**
-                                     * Virtual function which assembles the
-                                     * cell matrix and the right hand side
-                                     * on a given cell.
-                                     *
-                                     * This function assumes the cell matrix
-                                     * and right hand side to have the right
-                                     * size and to be empty. Functions of
-                                     * derived classes should check for
-                                     * this.
-                                     * For that purpose, the two exceptions
-                                     * @p{ExcWrongSize} and @p{ExcObjectNotEmpty}
-                                     * are declared.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &cell) const;
-
-                                    /**
-                                     * Virtual function which only assembles
-                                     * the cell matrix on a given cell.
-                                     *
-                                     * This function assumes the cell matrix
-                                     * and right hand side to have the right
-                                     * size and to be empty. Functions of
-                                     * derived classes should check for
-                                     * this.
-                                     * For that purpose, the two exceptions
-                                     * @p{ExcWrongSize} and @p{ExcObjectNotEmpty}
-                                     * are declared.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &cell) const;
-
-                                    /**
-                                     * Virtual function which only assembles
-                                     * the right hand side on a given cell.
-                                     *
-                                     * This function assumes the cell matrix
-                                     * and right hand side to have the right
-                                     * size and to be empty. Functions of
-                                     * derived classes should check for
-                                     * this.
-                                     * For that purpose, the two exceptions
-                                     * @p{ExcWrongSize} and @p{ExcObjectNotEmpty}
-                                     * are declared.
-                                     */
-    virtual void assemble (Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &cell) const;
-
-                                    /**
-                                     * Return number of equations for this
-                                     * equation object. This equals the number
-                                     * of solution functions.
-                                     */
-    unsigned int n_equations () const;
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcPureVirtualFunctionCalled);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcWrongSize,
-                   int, int,
-                   << "Object has wrong size " << arg1
-                   << ", but should have " << arg2 << ".");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcObjectNotEmpty);
-    
-  protected:
-                                    /**
-                                     * Store the number of solution functions,
-                                     * which is the same as the number of
-                                     * equations.
-                                     */
-    const unsigned int n_eq;
-};
-
-
-
-/**
- * The use of this class is now deprecated!
- *
- * An @p{Assembler} is a specialized version of a @p{DoFCellAccessor} which adds
- * functionality to assemble global matrices and vectors from cell base ones.
- *
- * @author Wolfgang Bangerth, 1998
- */
-template <int dim>
-class Assembler : public DoFCellAccessor<dim>
-{
-  public:
-
-                                    /**
-                                     * Structure to be passed upon
-                                     * construction of an assembler object
-                                     * through the iterator object. See
-                                     * @ref{TriaRawIterator} for a discussion
-                                     * of this mechanism.
-                                     */
-    struct AssemblerData {
-                                        /**
-                                         * Constructor.
-                                         */
-       AssemblerData (const DoFHandler<dim>    &dof,
-                      const bool                assemble_matrix,
-                      const bool                assemble_rhs,
-                      SparseMatrix<double>     &matrix,
-                      Vector<double>           &rhs_vector,
-                      const Quadrature<dim>    &quadrature,
-                      const UpdateFlags        &update_flags);
-       
-                                        /**
-                                         * Pointer to the dof handler object
-                                         * to be used to iterate on.
-                                         */
-       const DoFHandler<dim>  &dof;
-       
-                                        /**
-                                         * Flags to assemble the matrix.
-                                         */
-       const bool              assemble_matrix;
-
-                                        /**
-                                         * Flags whether to assemble the right hand sides.
-                                         */
-       const bool              assemble_rhs;
-       
-                                        /**
-                                         * Pointer to the matrix to be assembled
-                                         * by this object. Elements are summed
-                                         * up by the assembler, so you may want
-                                         * to clear this object (set all entries
-                                         * to zero) before use.
-                                         */
-       SparseMatrix<double>   &matrix;
-       
-                                        /**
-                                         * Pointer to the vector to be assembled
-                                         * by this object. Elements are summed
-                                         * up by the assembler, so you may want
-                                         * to clear this object (set all entries
-                                         * to zero) before use.
-                                         */
-       Vector<double>         &rhs_vector;
-       
-                                        /**
-                                         * Pointer to a quadrature object to be
-                                         * used for this assemblage process.
-                                         */
-       const Quadrature<dim>  &quadrature;
-       
-                                        /**
-                                         * Store which of the fields of the
-                                         * FEValues object need to be reinitialized
-                                         * on each cell.
-                                         */
-       const UpdateFlags       update_flags;
-    };
-
-
-                                    /**
-                                     * Declare the data type that this accessor
-                                     * class expects to get passed from the
-                                     * iterator classes.
-                                     */
-    typedef AssemblerData AccessorData;
-    
-                                    /**
-                                     * Default constructor, unused thus not
-                                     * implemented.
-                                     */
-    Assembler ();
-    
-                                    /**
-                                     * Constructor. The @p{local_data}
-                                     * argument is assumed to be a pointer
-                                     * to an @p{AssemblerData} object. The data
-                                     * is copied, so the object need not live
-                                     * longer than the constructor call.
-                                     */
-    Assembler (Triangulation<dim> *tria,
-              const int           level,
-              const int           index,
-              const AccessorData *local_data);
-    
-                                    /**
-                                     * Assemble on the present cell using
-                                     * the given equation objectand the data
-                                     * passed to the constructor. The elements
-                                     * of the local matrix and right hand side
-                                     * are added to the global matrix and
-                                     * vector so you may want to clear the
-                                     * matrix before use.
-                                     */
-    void assemble (const Equation<dim> &);
-
-                                    /**
-                                     * Exception.
-                                     */
-    DeclException0 (ExcNoAssemblingRequired);
-                                    /**
-                                     * Exception.
-                                     */
-    DeclException0 (ExcInvalidData);
-                                    /**
-                                     * Exception.
-                                     */
-                                    /**
-                                     * Exception.
-                                     */
-  private:
-                                    /**
-                                     * Store a local cell matrix.
-                                     */
-    FullMatrix<double>  cell_matrix;
-
-                                    /**
-                                     * Right hand side local to cell.
-                                     */
-    Vector<double>    cell_vector;
-
-                                    /**
-                                     * Store whether to assemble the
-                                     * global matrix.
-                                     */
-    bool              assemble_matrix;
-
-                                    /**
-                                     * Store whether to assemble the
-                                     * right hand side.
-                                     */
-    bool              assemble_rhs;
-
-                                    /**
-                                     * Pointer to the matrix to be assembled
-                                     * by this object.
-                                     */
-    SparseMatrix<double>    &matrix;
-
-                                    /**
-                                     * Pointer to the vector to be assembled
-                                     * by this object.
-                                     */
-    Vector<double>          &rhs_vector;
-
-                                    /**
-                                     * The finite element evaluated at the
-                                     * quadrature points.
-                                     */
-    FEValues<dim>     fe_values;
-};
-
-
-/*----------------------------   problem_assembler.h     ---------------------------*/
-
-#endif
-/*----------------------------   problem_assembler.h     ---------------------------*/
index 32e186486daee71caeabe925a4f774f15ef8ce9c..c3919beacf0cc83cd46289bc996208428f2df76c 100644 (file)
@@ -13,6 +13,8 @@
 #include <base/exceptions.h>
 #include <map>
 
+
+// forward declarations
 template <int dim> class Quadrature;
 
 template<typename number> class Vector;
@@ -25,7 +27,8 @@ template <typename Number> class BlockVector;
 template <int dim> class DoFHandler;
 template <int dim> class MGDoFHandler;
 template <int dim> class FEValues;
-template <int dim> class Equation;
+
+
 
 /**
  * Provide a class which assembles certain standard matrices for a given
@@ -85,8 +88,11 @@ template <int dim> class Equation;
  *   the different components. It will furthermore accept a single
  *   coefficient through the @ref{Function} parameter for all
  *   components. If you want different coefficients for the different
- *   parameters, you need to pass a function object representing the
- *   respective number of components.
+ *   parameters, you have to assemble the matrix yourself, sorry; the
+ *   implementation of the function will serve as a good starting
+ *   point, though. (You may also modify the implementation to accept
+ *   vector-valued functions and send this implementation to us -- we
+ *   will then include this implementation into the library.)
  *
  * @item @p{create_laplace_matrix}: there are two versions of this; the
  *   one which takes the @ref{Function} object creates
@@ -98,9 +104,10 @@ template <int dim> class Equation;
  *   This function uses the @ref{LaplaceMatrix} class.
  *
  *   If the finite element in use presently has more than only one
- *   component, this function may not be overly useful and presently
- *   throws an error.
- * @end{itemize}
+ *   component, this function may not be overly useful. It assembles a
+ *   Laplace matrix block for each component (with the same
+ *   coefficient for each component). These blocks are not coupled.  *
+ *   @end{itemize}
  *
  * All created matrices are `raw': they are not condensed, i.e. hanging
  * nodes are not eliminated. The reason is that you may want to add
@@ -294,6 +301,7 @@ class MatrixCreator
 };
 
 
+
 /**
  * Provide a collection of functions operating on matrices. These include
  * the application of boundary conditions to a linear system of equations
@@ -417,9 +425,9 @@ class MatrixTools : public MatrixCreator<dim>
                                      */
     static void
     apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                          BlockSparseMatrix<double> &matrix,
-                          BlockVector<double> &solution,
-                          BlockVector<double> &right_hand_side,
+                          BlockSparseMatrix<double>           &matrix,
+                          BlockVector<double>                 &solution,
+                          BlockVector<double>                 &right_hand_side,
                           const bool           eliminate_columns = true);
     
                                     /**
@@ -440,256 +448,5 @@ class MatrixTools : public MatrixCreator<dim>
 };
 
 
-/**
- * Equation class to be passed to the @ref{Assembler} if you want to make up the
- * mass matrix for your problem. The mass matrix is the matrix with
- * $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$.
- *
- * You may pass a coefficient function to the constructor. If you do so, the
- * assemble routines compute the matrix
- * $m_{ij} = \int_\Omega a(x) \phi_i(x) \phi_j(x) dx$
- * instead. The coefficient will in many cases be a strictly positive function.
- *
- * The class also has functions to create a right hand side
- * $f_i = \int_\Omega f(x) \phi_i(x) dx$. The function $f(x)$ has to be
- * given to the constructor; if none is given, an error is issued if you
- * try to create a right hand side vector. The function to create right
- * hand side vectors is the same for all the matrix class in this file,
- * since it does not depend on the operator.
- *
- * The defaults for both right hand side and coefficient function is a
- * @p{NULL} pointer. If you need a coefficient but no right hand side object,
- * simply pass a @p{NULL} pointer to the constructor for its first argument.
- *
- *
- * @sect3{Other possibilities}
- *
- * You will usually want to use this object only if you have coefficients
- * which vary over each cell. If you have coefficients which are constant
- * on each cell or even on the whole domain, you can get the local mass
- * matrix easier by calling the @ref{FiniteElement}@p{::get_local_mass_matrix} and
- * then scaling this one on each cell. This has the additional benefit that
- * the mass matrix is evaluated exactly, i.e. not using a quadrature formula
- * and is normally much faster since it can be precomputed and needs only
- * be scaled appropriately.
- *
- * The useful use of this object is therefore probable one of the following
- * cases:
- * @begin{itemize}
- * @item Mass lumping: use an @ref{Assembler} object and a special quadrature
- *   formula to voluntarily evaluate the mass matrix incorrect. For example
- *   by using the trapezoidal formula, the mass matrix will become a
- *   diagonal (at least if no hanging nodes are considered). However, there
- *   may be easier ways to set up the resulting matrix, for example by
- *   scaling the diagonal elements of the unit matrix by the area element
- *   of the respective cell.
- *
- * @item Nonconstant coefficient: if the coefficient varies considerably over
- *   each element, there is no way around this class. However, there are many
- *   cases where it is sufficient to assume that the function be constant on
- *   each cell (taking on its mean value throughout the cell for example, or
- *   more easily computed, its value at the center of mass of the element).
- *   A proper analysis of the error introduced by an assumed constant
- *   coefficient may be worth the effort.
- *
- *   Nonconstant coefficients to the mass matrix occur in mechanical problems
- *   if the density or other mechanical properties vary with the space
- *   coordinate.
- *
- * @item Simple plugging together of system matrices: if the system matrix has
- *    the form $s_{ij} = m_{ij} + \alpha a_{ij}$, for example, with $M$ and
- *    $A$ being the mass and laplace matrix, respectively (this matrix $S$
- *    occurs in the discretization of the heat and the wave equation, amoung
- *    others), once could conceive an equation object in which the @p{assemble}
- *    functions do nothing but sum up the contributions delivered by the
- *    @p{assemble} functions of the @ref{MassMatrix} and @ref{LaplaceMatrix} classes.
- *    Since numerical quadrature is necessary here anyway, this way is
- *    justifyable to quickly try something out. In the further process it
- *    may be useful to replace this behaviour by more sophisticated methods,
- *    however.
- * @end{itemize}
- */
-template <int dim>
-class MassMatrix :  public Equation<dim> {
-  public:
-                                    /**
-                                     * Constructor. Pass a function object if
-                                     * you want to create a right hand side
-                                     * vector, pass a function pointer (default
-                                     * is a NULL pointer). It is your duty to
-                                     * guarantee that the function object for
-                                     * the right hand side lives at least as
-                                     * long as this object does.
-                                     *
-                                     * You may also pass a function describing
-                                     * the weight to the integral (see the
-                                     * general docs for more information). The
-                                     * same applies for this object as said
-                                     * above.
-                                     */
-    MassMatrix (const Function<dim> * const rhs = 0,
-               const Function<dim> * const a = 0);
-
-                                    /**
-                                     * Assemble the cell matrix and right hand
-                                     * side vector for this cell. You need to
-                                     * give a right hand side object to the
-                                     * constructor to use this function. If
-                                     * a coefficient was given to the
-                                     * constructor, it is used.
-                                     *
-                                     * This function assumes the cell matrix
-                                     * and right hand side to have the right
-                                     * size and to be empty.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-
-                                    /**
-                                     * Construct the cell matrix for this cell.
-                                     * If a coefficient was given to the
-                                     * constructor, it is used.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-
-                                    /**
-                                     * Only construct the right hand side
-                                     * vector for this cell. You need to give
-                                     * a right hand side function to the
-                                     * constructor in order to call this
-                                     * function.
-                                     */
-    virtual void assemble (Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNoRHSSelected);
-    
-  protected:
-                                    /**
-                                     * Pointer to a function describing the
-                                     * right hand side of the problem. Should
-                                     * be zero if not given to the constructor
-                                     * and should then not be used.
-                                     */
-    const Function<dim> * const right_hand_side;
-
-                                    /**
-                                     * Pointer to a function describing the
-                                     * coefficient to the integral for the
-                                     * matrix entries. Should
-                                     * be zero if not given to the constructor
-                                     * and should then not be used.
-                                     */
-    const Function<dim> * const coefficient;
-};
-
-
-/**
- * Equation class to be passed to the @ref{Assembler} if you want to make up the
- * laplace matrix for your problem. The laplace matrix is the matrix with
- * $a_{ij} = \int_\Omega \nabla\phi_i(x) \cdot \nabla\phi_j(x) dx$.
- *
- * You may pass a coefficient function to the constructor. If you do so, the
- * assemble routines compute the matrix
- * $m_{ij} = \int_\Omega a(x) \nabla\phi_i(x) \cdot \nabla\phi_j(x) dx$
- * instead. The coefficient will in many cases be a strictly positive function.
- *
- * The class also has functions to create a right hand side
- * $f_i = \int_\Omega f(x) \phi_i(x) dx$. The function $f(x)$ has to be
- * given to the constructor; if none is given, an error is issued if you
- * try to create a right hand side vector. The function to create right
- * hand side vectors is the same for all the matrix class in this file,
- * since it does not depend on the operator.
- *
- * The defaults for both right hand side and coefficient function is a
- * @p{NULL} pointer. If you need a coefficient but no right hand side object,
- * simply pass a @p{NULL} pointer to the constructor for its first argument.
- */
-template <int dim>
-class LaplaceMatrix :  public Equation<dim> {
-  public:
-                                    /**
-                                     * Constructor. Pass a function object if
-                                     * you want to create a right hand side
-                                     * vector, pass a function pointer (default
-                                     * is a NULL pointer). It is your duty to
-                                     * guarantee that the function object for
-                                     * the right hand side lives at least as
-                                     * long as this object does.
-                                     *
-                                     * You may also pass a function describing
-                                     * the weight to the integral (see the
-                                     * general docs for more information). The
-                                     * same applies for this object as said
-                                     * above.
-                                     */
-    LaplaceMatrix (const Function<dim> * const rhs = 0,
-                  const Function<dim> * const a = 0);
-
-                                    /**
-                                     * Assemble the cell matrix and right hand
-                                     * side vector for this cell. You need to
-                                     * give a right hand side object to the
-                                     * constructor to use this function. If
-                                     * a coefficient was given to the
-                                     * constructor, it is used.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-
-                                    /**
-                                     * Construct the cell matrix for this cell.
-                                     * If a coefficient was given to the
-                                     * constructor, it is used.
-                                     */
-    virtual void assemble (FullMatrix<double>  &cell_matrix,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-
-                                    /**
-                                     * Only construct the right hand side
-                                     * vector for this cell. You need to give
-                                     * a right hand side function to the
-                                     * constructor in order to call this
-                                     * function.
-                                     */
-    virtual void assemble (Vector<double>      &rhs,
-                          const FEValues<dim> &fe_values,
-                          const typename DoFHandler<dim>::cell_iterator &) const;
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNoRHSSelected);
-    
-  protected:
-                                    /**
-                                     * Pointer to a function describing the
-                                     * right hand side of the problem. Should
-                                     * be zero if not given to the constructor
-                                     * and should then not be used.
-                                     */
-    const Function<dim> * const right_hand_side;
-
-                                    /**
-                                     * Pointer to a function describing the
-                                     * coefficient to the integral for the
-                                     * matrix entries. Should
-                                     * be zero if not given to the constructor
-                                     * and should then not be used.
-                                     */
-    const Function<dim> * const coefficient;
-};
-
 
 #endif
diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc
deleted file mode 100644 (file)
index 11b48d2..0000000
+++ /dev/null
@@ -1,132 +0,0 @@
-//----------------------------  assembler.cc  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------  assembler.cc  ---------------------------
-
-
-#include <numerics/assembler.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_iterator.templates.h>
-#include <fe/fe.h>
-#include <lac/full_matrix.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-#include <base/quadrature.h>
-#include <fe/mapping_q1.h>
-
-
-// if necessary try to work around a bug in the IBM xlC compiler
-#ifdef XLC_WORK_AROUND_STD_BUG
-using namespace std;
-#endif
-
-
-//TODO: purge this variable
-static const MappingQ1<deal_II_dimension> mapping;
-
-template <int dim>
-Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim>    &dof,
-                                             const bool                assemble_matrix,
-                                             const bool                assemble_rhs,
-                                             SparseMatrix<double>     &matrix,
-                                             Vector<double>           &rhs_vector,
-                                             const Quadrature<dim>    &quadrature,
-                                             const UpdateFlags        &update_flags) :
-               dof(dof),
-               assemble_matrix(assemble_matrix),
-               assemble_rhs(assemble_rhs),
-               matrix(matrix),
-               rhs_vector(rhs_vector),
-               quadrature(quadrature),
-               update_flags(update_flags)
-{};
-
-
-template <int dim>
-Assembler<dim>::Assembler (Triangulation<dim>  *tria,
-                          const int            level,
-                          const int            index,
-                          const AssemblerData *local_data) :
-               DoFCellAccessor<dim> (tria,level,index, &local_data->dof),
-               cell_matrix (dof_handler->get_fe().dofs_per_cell),
-               cell_vector (Vector<double>(dof_handler->get_fe().dofs_per_cell)),
-               assemble_matrix (local_data->assemble_matrix),
-               assemble_rhs (local_data->assemble_rhs),
-               matrix(local_data->matrix),
-               rhs_vector(local_data->rhs_vector),
-               fe_values (mapping, dof_handler->get_fe(),
-                          local_data->quadrature,
-                          local_data->update_flags)
-{
-  Assert (!assemble_matrix || (matrix.m() == dof_handler->n_dofs()),
-         ExcInvalidData());
-  Assert (!assemble_matrix || (matrix.n() == dof_handler->n_dofs()),
-         ExcInvalidData());
-  Assert (!assemble_rhs || (rhs_vector.size()==dof_handler->n_dofs()),
-         ExcInvalidData());
-};
-
-
-template <int dim>
-void Assembler<dim>::assemble (const Equation<dim> &equation) {
-                                  // re-init fe values for this cell
-  fe_values.reinit (DoFHandler<dim>::cell_iterator (*this));
-  const unsigned int n_dofs = dof_handler->get_fe().dofs_per_cell;
-
-  if (assemble_matrix)
-    cell_matrix.clear ();
-  if (assemble_rhs)
-    cell_vector.clear ();
-
-
-// fill cell matrix and vector if required
-  DoFHandler<dim>::cell_iterator this_cell (*this);
-  if (assemble_matrix && assemble_rhs) 
-    equation.assemble (cell_matrix, cell_vector, fe_values, this_cell);
-  else
-    if (assemble_matrix)
-      equation.assemble (cell_matrix, fe_values, this_cell);
-    else
-      if (assemble_rhs)
-       equation.assemble (cell_vector, fe_values, this_cell);
-      else
-       Assert (false, ExcNoAssemblingRequired());
-
-
-// get indices of dofs
-  std::vector<unsigned int> dofs (n_dofs);
-  get_dof_indices (dofs);
-
-                                  // one could use the
-                                  // @p{distribute_local_to_global} functions
-                                  // here, but they would require getting the
-                                  // dof indices twice, so we leave it the
-                                  // way it was originally programmed.
-  
-                                  // distribute cell matrix
-  if (assemble_matrix)
-    for (unsigned int i=0; i<n_dofs; ++i)
-      for (unsigned int j=0; j<n_dofs; ++j)
-       matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
-
-                                  // distribute cell vector
-  if (assemble_rhs)
-    for (unsigned int j=0; j<n_dofs; ++j)
-      rhs_vector(dofs[j]) += cell_vector(j);
-};
-
-
-// explicit instantiations
-template class Assembler<deal_II_dimension>;
-
-template class TriaRawIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
-template class TriaIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
-template class TriaActiveIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
diff --git a/deal.II/deal.II/source/numerics/equation.cc b/deal.II/deal.II/source/numerics/equation.cc
deleted file mode 100644 (file)
index af05025..0000000
+++ /dev/null
@@ -1,56 +0,0 @@
-//----------------------------  equation.cc  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------  equation.cc  ---------------------------
-
-
-#include <numerics/assembler.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_iterator.templates.h>
-#include <fe/fe.h>
-#include <lac/full_matrix.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
-template <int dim>
-Equation<dim>::Equation (const unsigned int n_equations) :
-               n_eq(n_equations) {};
-
-
-template <int dim>
-void Equation<dim>::assemble (FullMatrix<double>          &,
-                             Vector<double>           &,
-                             const FEValues<dim> &,
-                             const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-
-template <int dim>
-void Equation<dim>::assemble (FullMatrix<double>          &,
-                             const FEValues<dim> &,
-                             const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-
-template <int dim>
-void Equation<dim>::assemble (Vector<double>           &,
-                             const FEValues<dim> &,
-                             const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-
-template class Equation<deal_II_dimension>;
index d0f4c1abd8f076661775d7d5bd9e37dde454723e..004bac89b10077afce27f76c3d412c35d822bc1c 100644 (file)
@@ -21,7 +21,6 @@
 #include <fe/fe.h>
 #include <fe/fe_values.h>
 #include <numerics/matrices.h>
-#include <numerics/assembler.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
@@ -52,27 +51,68 @@ template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
-                                            const Function<dim> * const a)
+                                            const Function<dim> * const coefficient)
 {
-  Vector<double> dummy;    // no entries, should give an error if accessed
   UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
-  if (a != 0)
+  if (coefficient != 0)
     update_flags = UpdateFlags (update_flags | update_q_points);
-  const Assembler<dim>::AssemblerData data (dof,
-                                           true, false,  // assemble matrix but not rhs
-                                           matrix, dummy,
-                                           q, update_flags);
-  TriaActiveIterator<dim, Assembler<dim> >
-    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
-              dof.get_tria().begin_active()->level(),
-              dof.get_tria().begin_active()->index(),
-              &data);
-  MassMatrix<dim> equation(0,a);
-  do 
+
+  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+    
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<double> coefficient_values (n_q_points);
+  
+  std::vector<unsigned int> dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  for (; cell!=dof.end(); ++cell)
     {
-      assembler->assemble (equation);
-    }
-  while ((++assembler).state() == valid);
+      fe_values.reinit (cell);
+      
+      cell_matrix.clear ();
+      cell->get_dof_indices (dof_indices);
+      
+      const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+      const std::vector<double> &weights   = fe_values.get_JxW_values ();
+      
+      if (coefficient != 0)
+       {
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               if ((n_components==1) ||
+                   (fe.system_to_component_index(i).first ==
+                    fe.system_to_component_index(j).first))
+                 cell_matrix(i,j) += (values(i,point) *
+                                      values(j,point) *
+                                      weights[point] *
+                                      coefficient_values[point]);
+       }
+      else
+       for (unsigned int point=0; point<n_q_points; ++point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((n_components==1) ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               cell_matrix(i,j) += (values(i,point) *
+                                    values(j,point) *
+                                    weights[point]);
+
+                                      // transfer everything into the
+                                      // global object
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         matrix.add (dof_indices[i], dof_indices[j],
+                     cell_matrix(i,j));
+    };
 };
 
 
@@ -83,26 +123,86 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             SparseMatrix<double>     &matrix,
                                             const Function<dim>      &rhs,
                                             Vector<double>           &rhs_vector,
-                                            const Function<dim> * const a)
+                                            const Function<dim> * const coefficient)
 {
   UpdateFlags update_flags = UpdateFlags(update_values |
                                         update_q_points |
                                         update_JxW_values);
-  const Assembler<dim>::AssemblerData data (dof,
-                                           true, true,
-                                           matrix, rhs_vector,
-                                           q, update_flags);
-  TriaActiveIterator<dim, Assembler<dim> >
-    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
-              dof.get_tria().begin_active()->level(),
-              dof.get_tria().begin_active()->index(),
-              &data);
-  MassMatrix<dim> equation(&rhs,a);
-  do 
+  if (coefficient != 0)
+    update_flags = UpdateFlags (update_flags | update_q_points);
+
+  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+    
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>      local_rhs (dofs_per_cell);
+  std::vector<double> rhs_values (fe_values.n_quadrature_points);
+  std::vector<double> coefficient_values (n_q_points);
+  
+  std::vector<unsigned int> dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  for (; cell!=dof.end(); ++cell)
     {
-      assembler->assemble (equation);
-    }
-  while ((++assembler).state() == valid);
+      fe_values.reinit (cell);
+      
+      cell_matrix.clear ();
+      local_rhs.clear ();
+      cell->get_dof_indices (dof_indices);
+      
+      const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+      const std::vector<double> &weights   = fe_values.get_JxW_values ();
+      rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
+      
+      if (coefficient != 0)
+       {
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             {
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 if ((n_components==1) ||
+                     (fe.system_to_component_index(i).first ==
+                      fe.system_to_component_index(j).first))
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point] *
+                                        coefficient_values[point]);
+               local_rhs(i) += values(i,point) *
+                               rhs_values[point] *
+                               weights[point];
+             };
+       }
+      else
+       for (unsigned int point=0; point<n_q_points; ++point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           {
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               if ((n_components==1) ||
+                   (fe.system_to_component_index(i).first ==
+                    fe.system_to_component_index(j).first))
+                 cell_matrix(i,j) += (values(i,point) *
+                                      values(j,point) *
+                                      weights[point]);
+             local_rhs(i) += values(i,point) *
+                             rhs_values[point] *
+                             weights[point];
+           };
+
+                                      // transfer everything into the
+                                      // global object
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         matrix.add (dof_indices[i], dof_indices[j],
+                     cell_matrix(i,j));
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       rhs_vector(dof_indices[i]) += local_rhs(i);
+    };
 };
 
 
@@ -408,31 +508,70 @@ template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
                                                SparseMatrix<double>     &matrix,
-                                               const Function<dim> * const a)
+                                               const Function<dim> * const coefficient)
 {
-  const unsigned int n_components  = dof.get_fe().n_components();
-  Assert ((n_components==1) || (a==0), ExcNotImplemented());
+  UpdateFlags update_flags = UpdateFlags(update_JxW_values |
+                                        update_gradients);
+  if (coefficient != 0)
+    update_flags = UpdateFlags (update_flags | update_q_points);
 
-  Vector<double> dummy;   // no entries, should give an error if accessed
-  UpdateFlags update_flags = UpdateFlags(update_gradients |
-                                        update_JxW_values);
-  if (a != 0)
-    update_flags = UpdateFlags(update_flags | update_q_points);
-  const Assembler<dim>::AssemblerData data (dof,
-                                           true, false,  // assemble matrix but not rhs
-                                           matrix, dummy,
-                                           q, update_flags);
-  TriaActiveIterator<dim, Assembler<dim> >
-    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
-              dof.get_tria().begin_active()->level(),
-              dof.get_tria().begin_active()->index(),
-              &data);
-  LaplaceMatrix<dim> equation (0, a);
-  do 
+  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+    
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<double> coefficient_values (n_q_points);
+  
+  std::vector<unsigned int> dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  for (; cell!=dof.end(); ++cell)
     {
-      assembler->assemble (equation);
-    }
-  while ((++assembler).state() == valid);
+      fe_values.reinit (cell);
+      
+      cell_matrix.clear ();
+      cell->get_dof_indices (dof_indices);
+      
+      const std::vector<std::vector<Tensor<1,dim> > >
+       &grads   = fe_values.get_shape_grads ();
+      const std::vector<double> &weights = fe_values.get_JxW_values ();
+      
+      if (coefficient != 0)
+       {
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               if ((n_components==1) ||
+                   (fe.system_to_component_index(i).first ==
+                    fe.system_to_component_index(j).first))
+                 cell_matrix(i,j) += (grads[i][point] *
+                                      grads[j][point] *
+                                      weights[point] *
+                                      coefficient_values[point]);
+       }
+      else
+       for (unsigned int point=0; point<n_q_points; ++point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((n_components==1) ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               cell_matrix(i,j) += (grads[i][point] *
+                                    grads[j][point] *
+                                    weights[point]);
+
+                                      // transfer everything into the
+                                      // global object
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         matrix.add (dof_indices[i], dof_indices[j],
+                     cell_matrix(i,j));
+    };
 };
 
 
@@ -483,33 +622,94 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                SparseMatrix<double>     &matrix,
                                                const Function<dim>      &rhs,
                                                Vector<double>           &rhs_vector,
-                                               const Function<dim> * const a)
+                                               const Function<dim> * const coefficient)
 {
-  const unsigned int n_components  = dof.get_fe().n_components();
-  Assert ((n_components==1) || (a==0), ExcNotImplemented());
-
-  UpdateFlags update_flags = UpdateFlags(update_q_points  |
+  UpdateFlags update_flags = UpdateFlags(update_values    |
                                         update_gradients |
+                                        update_q_points  |
                                         update_JxW_values);
-  const Assembler<dim>::AssemblerData data (dof,
-                                           true, true,
-                                           matrix, rhs_vector,
-                                           q, update_flags);
-  TriaActiveIterator<dim, Assembler<dim> >
-    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
-              dof.get_tria().begin_active()->level(),
-              dof.get_tria().begin_active()->index(),
-              &data);
-  LaplaceMatrix<dim> equation (&rhs, a);
-  do 
+  if (coefficient != 0)
+    update_flags = UpdateFlags (update_flags | update_q_points);
+
+  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+    
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>      local_rhs (dofs_per_cell);
+  std::vector<double> rhs_values (fe_values.n_quadrature_points);
+  std::vector<double> coefficient_values (n_q_points);
+  
+  std::vector<unsigned int> dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  for (; cell!=dof.end(); ++cell)
     {
-      assembler->assemble (equation);
-    }
-  while ((++assembler).state() == valid);
+      fe_values.reinit (cell);
+      
+      cell_matrix.clear ();
+      local_rhs.clear ();
+      cell->get_dof_indices (dof_indices);
+      
+      const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+      const std::vector<std::vector<Tensor<1,dim> > >
+       &grads   = fe_values.get_shape_grads ();
+      const std::vector<double> &weights   = fe_values.get_JxW_values ();
+      rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
+      
+      if (coefficient != 0)
+       {
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             {
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 if ((n_components==1) ||
+                     (fe.system_to_component_index(i).first ==
+                      fe.system_to_component_index(j).first))
+                   cell_matrix(i,j) += (grads[i][point] *
+                                        grads[j][point] *
+                                        weights[point] *
+                                        coefficient_values[point]);
+               local_rhs(i) += values(i,point) *
+                               rhs_values[point] *
+                               weights[point];
+             };
+       }
+      else
+       for (unsigned int point=0; point<n_q_points; ++point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           {
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               if ((n_components==1) ||
+                   (fe.system_to_component_index(i).first ==
+                    fe.system_to_component_index(j).first))
+                 cell_matrix(i,j) += (grads[i][point] *
+                                      grads[j][point] *
+                                      weights[point]);
+             local_rhs(i) += values(i,point) *
+                             rhs_values[point] *
+                             weights[point];
+           };
+
+                                      // transfer everything into the
+                                      // global object
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         matrix.add (dof_indices[i], dof_indices[j],
+                     cell_matrix(i,j));
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       rhs_vector(dof_indices[i]) += local_rhs(i);
+    };
 };
 
 
 
+
 template <int dim>
 template <typename number>
 void
@@ -996,383 +1196,12 @@ MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &bo
 
 
 
-template <int dim>
-MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
-                            const Function<dim> * const a) :
-               Equation<dim> (1),
-               right_hand_side (rhs),
-               coefficient (a)
-{};
-
-
-
-template <int dim>
-void MassMatrix<dim>::assemble (FullMatrix<double>      &cell_matrix,
-                               const FEValues<dim>     &fe_values,
-                               const typename DoFHandler<dim>::cell_iterator &) const
-{
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-  Assert (cell_matrix.n() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
-  Assert (cell_matrix.m() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
-  Assert (cell_matrix.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-  
-  const FullMatrix<double> &values    = fe_values.get_shape_values ();
-  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
-
-
-  if (coefficient != 0)
-    {
-      if (coefficient->n_components == 1)
-                                        // scalar coefficient given
-       {
-         std::vector<double> coefficient_values (fe_values.n_quadrature_points);
-         coefficient->value_list (fe_values.get_quadrature_points(),
-                                  coefficient_values);
-         for (unsigned int i=0; i<dofs_per_cell; ++i) 
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             if ((n_components == 1)
-                 ||
-                 (fe.system_to_component_index(i).first ==
-                  fe.system_to_component_index(j).first))
-               {
-                 for (unsigned int point=0; point<n_q_points; ++point)
-                   cell_matrix(i,j) += (values(i,point) *
-                                        values(j,point) *
-                                        weights[point] *
-                                        coefficient_values[point]);
-               };
-       }
-      else
-                                        // vectorial coefficient
-                                        // given
-       {
-         std::vector<Vector<double> > coefficient_values (fe_values.n_quadrature_points,
-                                                     Vector<double>(n_components));
-         coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                         coefficient_values);
-         for (unsigned int i=0; i<dofs_per_cell; ++i) 
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             if ((n_components == 1)
-                 ||
-                 (fe.system_to_component_index(i).first ==
-                  fe.system_to_component_index(j).first))
-               {
-                 for (unsigned int point=0; point<n_q_points; ++point)
-                   cell_matrix(i,j) += (values(i,point) *
-                                        values(j,point) *
-                                        weights[point] *
-                                        coefficient_values[point](
-                                          fe.system_to_component_index(i).first));
-               };
-       };
-      
-    }
-  else
-                                    // no coefficient given
-    for (unsigned int i=0; i<dofs_per_cell; ++i) 
-      for (unsigned int j=0; j<dofs_per_cell; ++j)
-       if ((n_components == 1)
-           ||
-           (fe.system_to_component_index(i).first ==
-            fe.system_to_component_index(j).first))
-         {
-           for (unsigned int point=0; point<n_q_points; ++point)
-             cell_matrix(i,j) += (values(i,point) *
-                                  values(j,point) *
-                                  weights[point]);
-         };
-};
-
-
-
-template <int dim>
-void MassMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
-                               Vector<double>      &rhs,
-                               const FEValues<dim> &fe_values,
-                               const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (right_hand_side != 0, ExcNoRHSSelected());
-
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-                                  // for system elements: not
-                                  // implemented at present
-  Assert (n_components==1, ExcNotImplemented());
-  
-  Assert (cell_matrix.n() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
-  Assert (cell_matrix.m() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
-  Assert (rhs.size() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
-  Assert (cell_matrix.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-  Assert (rhs.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-
-  const FullMatrix<double> &values    = fe_values.get_shape_values ();
-  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
-  std::vector<double>            rhs_values (fe_values.n_quadrature_points);
-  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
-
-  if (coefficient != 0)
-    {
-      std::vector<double> coefficient_values (n_q_points);
-      coefficient->value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-      for (unsigned int point=0; point<n_q_points; ++point)
-       for (unsigned int i=0; i<dofs_per_cell; ++i) 
-         {
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (values(i,point) *
-                                  values(j,point) *
-                                  weights[point] *
-                                  coefficient_values[point]);
-           rhs(i) += values(i,point) *
-                     rhs_values[point] *
-                     weights[point];
-         };
-    }
-  else
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<dofs_per_cell; ++i) 
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           cell_matrix(i,j) += (values(i,point) *
-                                values(j,point) *
-                                weights[point]);
-         rhs(i) += values(i,point) *
-                   rhs_values[point] *
-                   weights[point];
-       };
-};
-
-
-
-template <int dim>
-void MassMatrix<dim>::assemble (Vector<double>      &rhs,
-                               const FEValues<dim> &fe_values,
-                               const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (right_hand_side != 0, ExcNoRHSSelected());
-
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-                                  // for system elements: not
-                                  // implemented at present
-  Assert (n_components==1, ExcNotImplemented());
-
-  Assert (rhs.size() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
-  Assert (rhs.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-
-  const FullMatrix<double> &values    = fe_values.get_shape_values ();
-  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
-  std::vector<double>            rhs_values(fe_values.n_quadrature_points);
-  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
-
-  for (unsigned int point=0; point<n_q_points; ++point)
-    for (unsigned int i=0; i<dofs_per_cell; ++i) 
-      rhs(i) += values(i,point) *
-               rhs_values[point] *
-               weights[point];
-};
-
-
-
-
-
-template <int dim>
-LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
-                                  const Function<dim> * const a) :
-               Equation<dim> (1),
-               right_hand_side (rhs),
-               coefficient (a) {};
-
-
-
-template <int dim>
-void LaplaceMatrix<dim>::assemble (FullMatrix<double>         &cell_matrix,
-                                  Vector<double>             &rhs,
-                                  const FEValues<dim>        &fe_values,
-                                  const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (right_hand_side != 0, ExcNoRHSSelected());
-  
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-                                  // for system elements: might be
-                                  // not so useful, not implemented
-                                  // at present
-  Assert (n_components==1, ExcNotImplemented());
-
-  Assert (cell_matrix.n() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
-  Assert (cell_matrix.m() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
-  Assert (rhs.size() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
-  Assert (cell_matrix.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-  Assert (rhs.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-
-  const std::vector<std::vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
-  const FullMatrix<double>             &values    = fe_values.get_shape_values ();
-  std::vector<double>        rhs_values(fe_values.n_quadrature_points);
-  const std::vector<double> &weights   = fe_values.get_JxW_values ();
-  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
-
-  if (coefficient != 0)
-    {
-      std::vector<double> coefficient_values(n_q_points);
-      coefficient->value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-      for (unsigned int point=0; point<n_q_points; ++point)
-       for (unsigned int i=0; i<dofs_per_cell; ++i) 
-         {
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (gradients[i][point] *
-                                  gradients[j][point]) *
-                                 weights[point] *
-                                 coefficient_values[point];
-           rhs(i) += values(i,point) *
-                     rhs_values[point] *
-                     weights[point];
-         };
-    }
-  else
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<dofs_per_cell; ++i) 
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           cell_matrix(i,j) += (gradients[i][point] *
-                                gradients[j][point]) *
-                               weights[point];
-         rhs(i) += values(i,point) *
-                   rhs_values[point] *
-                   weights[point];
-       };
-
-};
-
-
-
-template <int dim>
-void LaplaceMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
-                                  const FEValues<dim> &fe_values,
-                                  const typename DoFHandler<dim>::cell_iterator &) const
-{
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-                                  // for system elements: might be
-                                  // not so useful, not implemented
-                                  // at present
-  Assert ((n_components==1) || (coefficient==0), ExcNotImplemented());
-
-  Assert (cell_matrix.n() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
-  Assert (cell_matrix.m() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
-  Assert (cell_matrix.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-  
-  const std::vector<std::vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
-  const std::vector<double> &weights   = fe_values.get_JxW_values ();
-   
-  if (coefficient != 0)
-    {
-      std::vector<double> coefficient_values(n_q_points);
-      coefficient->value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-      for (unsigned int point=0; point<n_q_points; ++point)
-       for (unsigned int i=0; i<dofs_per_cell; ++i) 
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           cell_matrix(i,j) += (gradients[i][point] *
-                                gradients[j][point]) *
-                               weights[point] *
-                               coefficient_values[point];
-    }
-  else
-    for (unsigned int i=0; i<dofs_per_cell; ++i) 
-      for (unsigned int j=0; j<dofs_per_cell; ++j)
-       if ((n_components==1)
-           ||
-           (fe.system_to_component_index(i).first ==
-            fe.system_to_component_index(j).first))
-         {
-           for (unsigned int point=0; point<n_q_points; ++point)
-             cell_matrix(i,j) += (gradients[i][point] *
-                                  gradients[j][point]) *
-                                 weights[point];
-         };
-};
-
-
-
-template <int dim>
-void LaplaceMatrix<dim>::assemble (Vector<double>      &rhs,
-                                  const FEValues<dim> &fe_values,
-                                  const typename DoFHandler<dim>::cell_iterator &) const
-{
-  Assert (right_hand_side != 0, ExcNoRHSSelected());
-
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-                                  // for system elements: might be
-                                  // not so useful, not implemented
-                                  // at present
-  Assert (n_components==1, ExcNotImplemented());
-
-  Assert (rhs.size() == dofs_per_cell,
-         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
-  Assert (rhs.all_zero(),
-         Equation<dim>::ExcObjectNotEmpty());
-
-  const FullMatrix<double> &values    = fe_values.get_shape_values ();
-  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
-  std::vector<double>        rhs_values(fe_values.n_quadrature_points);
-  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
-   
-  for (unsigned int point=0; point<n_q_points; ++point)
-    for (unsigned int i=0; i<dofs_per_cell; ++i) 
-      rhs(i) += values(i,point) *
-               rhs_values[point] *
-               weights[point];
-};
-
 
 
 // explicit instantiations
 
 template class MatrixCreator<deal_II_dimension>;
 template class MatrixTools<deal_II_dimension>;
-template class MassMatrix<deal_II_dimension>;
-template class LaplaceMatrix<deal_II_dimension>;
 
 
 template
index 95cb3a9f8a7fe08711678fbd635108de2d2a2f06..c41ef605e24c4b273d6d7c043befb5ffb7758fe5 100644 (file)
@@ -20,7 +20,6 @@
 #include <fe/fe.h>
 #include <fe/fe_values.h>
 #include <base/quadrature.h>
-#include <numerics/assembler.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <dofs/dof_tools.h>
index 886e1892038be6b13ef9cd41e4643ea931b036ec..995dc8eefcff7dedd5d77218c0be29fea2a1528e 100644 (file)
@@ -283,7 +283,12 @@ documentation, etc</a>.
   <li> <p>
        Removed: the <code class="class">ProblemBase</code> class,
        which has been deprecated since before the release of 
-       <acronym>deal.II</acronym> 3.0 has finally been removed.
+       <acronym>deal.II</acronym> 3.0 has finally been removed. The
+       same applied for the classes
+       <code class="class">Assembler</code>,
+       <code class="class">Equation</code>,
+       <code class="class">MassMatrix</code>, and
+       <code class="class">LaplaceMatrix</code>.
        <br>
        (WB 2001/03/27)
        </p>
index ce4ea0c407fca1a376be3a057d6cbcb072125356..021381e3943328d823ff514e607803a027adbf7a 100644 (file)
@@ -19,7 +19,6 @@
 #include <fe/fe_lib.lagrange.h>
 #include <base/quadrature_lib.h>
 #include "../problem_base.h"
-#include <numerics/assembler.h>
 #include <lac/sparse_matrix.h>
 
 
index fa5b32c7418605a05db012daf2f8c574f2f00a91..ddcfed4313e3013130ac98683c144a3009610f84 100644 (file)
@@ -249,8 +249,752 @@ class ProblemBase
 
 
 
+/**
+ * Equation class to be passed to the @ref{Assembler} if you want to make up the
+ * mass matrix for your problem. The mass matrix is the matrix with
+ * $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$.
+ *
+ * You may pass a coefficient function to the constructor. If you do so, the
+ * assemble routines compute the matrix
+ * $m_{ij} = \int_\Omega a(x) \phi_i(x) \phi_j(x) dx$
+ * instead. The coefficient will in many cases be a strictly positive function.
+ *
+ * The class also has functions to create a right hand side
+ * $f_i = \int_\Omega f(x) \phi_i(x) dx$. The function $f(x)$ has to be
+ * given to the constructor; if none is given, an error is issued if you
+ * try to create a right hand side vector. The function to create right
+ * hand side vectors is the same for all the matrix class in this file,
+ * since it does not depend on the operator.
+ *
+ * The defaults for both right hand side and coefficient function is a
+ * @p{NULL} pointer. If you need a coefficient but no right hand side object,
+ * simply pass a @p{NULL} pointer to the constructor for its first argument.
+ *
+ *
+ * @sect3{Other possibilities}
+ *
+ * You will usually want to use this object only if you have coefficients
+ * which vary over each cell. If you have coefficients which are constant
+ * on each cell or even on the whole domain, you can get the local mass
+ * matrix easier by calling the @ref{FiniteElement}@p{::get_local_mass_matrix} and
+ * then scaling this one on each cell. This has the additional benefit that
+ * the mass matrix is evaluated exactly, i.e. not using a quadrature formula
+ * and is normally much faster since it can be precomputed and needs only
+ * be scaled appropriately.
+ *
+ * The useful use of this object is therefore probable one of the following
+ * cases:
+ * @begin{itemize}
+ * @item Mass lumping: use an @ref{Assembler} object and a special quadrature
+ *   formula to voluntarily evaluate the mass matrix incorrect. For example
+ *   by using the trapezoidal formula, the mass matrix will become a
+ *   diagonal (at least if no hanging nodes are considered). However, there
+ *   may be easier ways to set up the resulting matrix, for example by
+ *   scaling the diagonal elements of the unit matrix by the area element
+ *   of the respective cell.
+ *
+ * @item Nonconstant coefficient: if the coefficient varies considerably over
+ *   each element, there is no way around this class. However, there are many
+ *   cases where it is sufficient to assume that the function be constant on
+ *   each cell (taking on its mean value throughout the cell for example, or
+ *   more easily computed, its value at the center of mass of the element).
+ *   A proper analysis of the error introduced by an assumed constant
+ *   coefficient may be worth the effort.
+ *
+ *   Nonconstant coefficients to the mass matrix occur in mechanical problems
+ *   if the density or other mechanical properties vary with the space
+ *   coordinate.
+ *
+ * @item Simple plugging together of system matrices: if the system matrix has
+ *    the form $s_{ij} = m_{ij} + \alpha a_{ij}$, for example, with $M$ and
+ *    $A$ being the mass and laplace matrix, respectively (this matrix $S$
+ *    occurs in the discretization of the heat and the wave equation, amoung
+ *    others), once could conceive an equation object in which the @p{assemble}
+ *    functions do nothing but sum up the contributions delivered by the
+ *    @p{assemble} functions of the @ref{MassMatrix} and @ref{LaplaceMatrix} classes.
+ *    Since numerical quadrature is necessary here anyway, this way is
+ *    justifyable to quickly try something out. In the further process it
+ *    may be useful to replace this behaviour by more sophisticated methods,
+ *    however.
+ * @end{itemize}
+ */
+template <int dim>
+class MassMatrix :  public Equation<dim> {
+  public:
+                                    /**
+                                     * Constructor. Pass a function object if
+                                     * you want to create a right hand side
+                                     * vector, pass a function pointer (default
+                                     * is a NULL pointer). It is your duty to
+                                     * guarantee that the function object for
+                                     * the right hand side lives at least as
+                                     * long as this object does.
+                                     *
+                                     * You may also pass a function describing
+                                     * the weight to the integral (see the
+                                     * general docs for more information). The
+                                     * same applies for this object as said
+                                     * above.
+                                     */
+    MassMatrix (const Function<dim> * const rhs = 0,
+               const Function<dim> * const a = 0);
+
+                                    /**
+                                     * Assemble the cell matrix and right hand
+                                     * side vector for this cell. You need to
+                                     * give a right hand side object to the
+                                     * constructor to use this function. If
+                                     * a coefficient was given to the
+                                     * constructor, it is used.
+                                     *
+                                     * This function assumes the cell matrix
+                                     * and right hand side to have the right
+                                     * size and to be empty.
+                                     */
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+
+                                    /**
+                                     * Construct the cell matrix for this cell.
+                                     * If a coefficient was given to the
+                                     * constructor, it is used.
+                                     */
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+
+                                    /**
+                                     * Only construct the right hand side
+                                     * vector for this cell. You need to give
+                                     * a right hand side function to the
+                                     * constructor in order to call this
+                                     * function.
+                                     */
+    virtual void assemble (Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+    
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNoRHSSelected);
+    
+  protected:
+                                    /**
+                                     * Pointer to a function describing the
+                                     * right hand side of the problem. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const right_hand_side;
+
+                                    /**
+                                     * Pointer to a function describing the
+                                     * coefficient to the integral for the
+                                     * matrix entries. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const coefficient;
+};
+
+
+/**
+ * Equation class to be passed to the @ref{Assembler} if you want to make up the
+ * laplace matrix for your problem. The laplace matrix is the matrix with
+ * $a_{ij} = \int_\Omega \nabla\phi_i(x) \cdot \nabla\phi_j(x) dx$.
+ *
+ * You may pass a coefficient function to the constructor. If you do so, the
+ * assemble routines compute the matrix
+ * $m_{ij} = \int_\Omega a(x) \nabla\phi_i(x) \cdot \nabla\phi_j(x) dx$
+ * instead. The coefficient will in many cases be a strictly positive function.
+ *
+ * The class also has functions to create a right hand side
+ * $f_i = \int_\Omega f(x) \phi_i(x) dx$. The function $f(x)$ has to be
+ * given to the constructor; if none is given, an error is issued if you
+ * try to create a right hand side vector. The function to create right
+ * hand side vectors is the same for all the matrix class in this file,
+ * since it does not depend on the operator.
+ *
+ * The defaults for both right hand side and coefficient function is a
+ * @p{NULL} pointer. If you need a coefficient but no right hand side object,
+ * simply pass a @p{NULL} pointer to the constructor for its first argument.
+ */
+template <int dim>
+class LaplaceMatrix :  public Equation<dim> {
+  public:
+                                    /**
+                                     * Constructor. Pass a function object if
+                                     * you want to create a right hand side
+                                     * vector, pass a function pointer (default
+                                     * is a NULL pointer). It is your duty to
+                                     * guarantee that the function object for
+                                     * the right hand side lives at least as
+                                     * long as this object does.
+                                     *
+                                     * You may also pass a function describing
+                                     * the weight to the integral (see the
+                                     * general docs for more information). The
+                                     * same applies for this object as said
+                                     * above.
+                                     */
+    LaplaceMatrix (const Function<dim> * const rhs = 0,
+                  const Function<dim> * const a = 0);
+
+                                    /**
+                                     * Assemble the cell matrix and right hand
+                                     * side vector for this cell. You need to
+                                     * give a right hand side object to the
+                                     * constructor to use this function. If
+                                     * a coefficient was given to the
+                                     * constructor, it is used.
+                                     */
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+
+                                    /**
+                                     * Construct the cell matrix for this cell.
+                                     * If a coefficient was given to the
+                                     * constructor, it is used.
+                                     */
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+
+                                    /**
+                                     * Only construct the right hand side
+                                     * vector for this cell. You need to give
+                                     * a right hand side function to the
+                                     * constructor in order to call this
+                                     * function.
+                                     */
+    virtual void assemble (Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename DoFHandler<dim>::cell_iterator &) const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNoRHSSelected);
+    
+  protected:
+                                    /**
+                                     * Pointer to a function describing the
+                                     * right hand side of the problem. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const right_hand_side;
+
+                                    /**
+                                     * Pointer to a function describing the
+                                     * coefficient to the integral for the
+                                     * matrix entries. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const coefficient;
+};
+
+
+
+
+
+template <int dim>
+MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
+                            const Function<dim> * const a) :
+               Equation<dim> (1),
+               right_hand_side (rhs),
+               coefficient (a)
+{};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (FullMatrix<double>      &cell_matrix,
+                               const FEValues<dim>     &fe_values,
+                               const typename DoFHandler<dim>::cell_iterator &) const
+{
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  Assert (cell_matrix.n() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
+  Assert (cell_matrix.m() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
+  Assert (cell_matrix.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+  
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+
+
+  if (coefficient != 0)
+    {
+      if (coefficient->n_components == 1)
+                                        // scalar coefficient given
+       {
+         std::vector<double> coefficient_values (fe_values.n_quadrature_points);
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((n_components == 1)
+                 ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               {
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point] *
+                                        coefficient_values[point]);
+               };
+       }
+      else
+                                        // vectorial coefficient
+                                        // given
+       {
+         std::vector<Vector<double> > coefficient_values (fe_values.n_quadrature_points,
+                                                     Vector<double>(n_components));
+         coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                         coefficient_values);
+         for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((n_components == 1)
+                 ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               {
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point] *
+                                        coefficient_values[point](
+                                          fe.system_to_component_index(i).first));
+               };
+       };
+      
+    }
+  else
+                                    // no coefficient given
+    for (unsigned int i=0; i<dofs_per_cell; ++i) 
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+       if ((n_components == 1)
+           ||
+           (fe.system_to_component_index(i).first ==
+            fe.system_to_component_index(j).first))
+         {
+           for (unsigned int point=0; point<n_q_points; ++point)
+             cell_matrix(i,j) += (values(i,point) *
+                                  values(j,point) *
+                                  weights[point]);
+         };
+};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
+                               Vector<double>      &rhs,
+                               const FEValues<dim> &fe_values,
+                               const typename DoFHandler<dim>::cell_iterator &) const
+{
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+                                  // for system elements: not
+                                  // implemented at present
+  Assert (n_components==1, ExcNotImplemented());
+  
+  Assert (cell_matrix.n() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
+  Assert (cell_matrix.m() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
+  Assert (rhs.size() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
+  Assert (cell_matrix.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+  Assert (rhs.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+  std::vector<double>            rhs_values (fe_values.n_quadrature_points);
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  if (coefficient != 0)
+    {
+      std::vector<double> coefficient_values (n_q_points);
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<n_q_points; ++point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i) 
+         {
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             cell_matrix(i,j) += (values(i,point) *
+                                  values(j,point) *
+                                  weights[point] *
+                                  coefficient_values[point]);
+           rhs(i) += values(i,point) *
+                     rhs_values[point] *
+                     weights[point];
+         };
+    }
+  else
+    for (unsigned int point=0; point<n_q_points; ++point)
+      for (unsigned int i=0; i<dofs_per_cell; ++i) 
+       {
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           cell_matrix(i,j) += (values(i,point) *
+                                values(j,point) *
+                                weights[point]);
+         rhs(i) += values(i,point) *
+                   rhs_values[point] *
+                   weights[point];
+       };
+};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (Vector<double>      &rhs,
+                               const FEValues<dim> &fe_values,
+                               const typename DoFHandler<dim>::cell_iterator &) const
+{
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+                                  // for system elements: not
+                                  // implemented at present
+  Assert (n_components==1, ExcNotImplemented());
+
+  Assert (rhs.size() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
+  Assert (rhs.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+  std::vector<double>            rhs_values(fe_values.n_quadrature_points);
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  for (unsigned int point=0; point<n_q_points; ++point)
+    for (unsigned int i=0; i<dofs_per_cell; ++i) 
+      rhs(i) += values(i,point) *
+               rhs_values[point] *
+               weights[point];
+};
+
+
+
+
+
+template <int dim>
+LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
+                                  const Function<dim> * const a) :
+               Equation<dim> (1),
+               right_hand_side (rhs),
+               coefficient (a) {};
+
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (FullMatrix<double>         &cell_matrix,
+                                  Vector<double>             &rhs,
+                                  const FEValues<dim>        &fe_values,
+                                  const typename DoFHandler<dim>::cell_iterator &) const
+{
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+  
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+                                  // for system elements: might be
+                                  // not so useful, not implemented
+                                  // at present
+  Assert (n_components==1, ExcNotImplemented());
+
+  Assert (cell_matrix.n() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
+  Assert (cell_matrix.m() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
+  Assert (rhs.size() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
+  Assert (cell_matrix.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+  Assert (rhs.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+
+  const std::vector<std::vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
+  const FullMatrix<double>             &values    = fe_values.get_shape_values ();
+  std::vector<double>        rhs_values(fe_values.n_quadrature_points);
+  const std::vector<double> &weights   = fe_values.get_JxW_values ();
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  if (coefficient != 0)
+    {
+      std::vector<double> coefficient_values(n_q_points);
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<n_q_points; ++point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i) 
+         {
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             cell_matrix(i,j) += (gradients[i][point] *
+                                  gradients[j][point]) *
+                                 weights[point] *
+                                 coefficient_values[point];
+           rhs(i) += values(i,point) *
+                     rhs_values[point] *
+                     weights[point];
+         };
+    }
+  else
+    for (unsigned int point=0; point<n_q_points; ++point)
+      for (unsigned int i=0; i<dofs_per_cell; ++i) 
+       {
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           cell_matrix(i,j) += (gradients[i][point] *
+                                gradients[j][point]) *
+                               weights[point];
+         rhs(i) += values(i,point) *
+                   rhs_values[point] *
+                   weights[point];
+       };
+
+};
+
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
+                                  const FEValues<dim> &fe_values,
+                                  const typename DoFHandler<dim>::cell_iterator &) const
+{
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+                                  // for system elements: might be
+                                  // not so useful, not implemented
+                                  // at present
+  Assert ((n_components==1) || (coefficient==0), ExcNotImplemented());
+
+  Assert (cell_matrix.n() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.n(), dofs_per_cell));
+  Assert (cell_matrix.m() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(cell_matrix.m(), dofs_per_cell));
+  Assert (cell_matrix.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+  
+  const std::vector<std::vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
+  const std::vector<double> &weights   = fe_values.get_JxW_values ();
+   
+  if (coefficient != 0)
+    {
+      std::vector<double> coefficient_values(n_q_points);
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<n_q_points; ++point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i) 
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           cell_matrix(i,j) += (gradients[i][point] *
+                                gradients[j][point]) *
+                               weights[point] *
+                               coefficient_values[point];
+    }
+  else
+    for (unsigned int i=0; i<dofs_per_cell; ++i) 
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+       if ((n_components==1)
+           ||
+           (fe.system_to_component_index(i).first ==
+            fe.system_to_component_index(j).first))
+         {
+           for (unsigned int point=0; point<n_q_points; ++point)
+             cell_matrix(i,j) += (gradients[i][point] *
+                                  gradients[j][point]) *
+                                 weights[point];
+         };
+};
+
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (Vector<double>      &rhs,
+                                  const FEValues<dim> &fe_values,
+                                  const typename DoFHandler<dim>::cell_iterator &) const
+{
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+                                  // for system elements: might be
+                                  // not so useful, not implemented
+                                  // at present
+  Assert (n_components==1, ExcNotImplemented());
+
+  Assert (rhs.size() == dofs_per_cell,
+         Equation<dim>::ExcWrongSize(rhs.size(), dofs_per_cell));
+  Assert (rhs.all_zero(),
+         Equation<dim>::ExcObjectNotEmpty());
+
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+  std::vector<double>        rhs_values(fe_values.n_quadrature_points);
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+   
+  for (unsigned int point=0; point<n_q_points; ++point)
+    for (unsigned int i=0; i<dofs_per_cell; ++i) 
+      rhs(i) += values(i,point) *
+               rhs_values[point] *
+               weights[point];
+};
+
+
+
+#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
+#include <fe/fe.h>
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <base/quadrature.h>
+#include <fe/mapping_q1.h>
+
+
+// if necessary try to work around a bug in the IBM xlC compiler
+#ifdef XLC_WORK_AROUND_STD_BUG
+using namespace std;
+#endif
+
+
+//TODO: purge this variable
+static const MappingQ1<deal_II_dimension> mapping;
+
+template <int dim>
+Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim>    &dof,
+                                             const bool                assemble_matrix,
+                                             const bool                assemble_rhs,
+                                             SparseMatrix<double>     &matrix,
+                                             Vector<double>           &rhs_vector,
+                                             const Quadrature<dim>    &quadrature,
+                                             const UpdateFlags        &update_flags) :
+               dof(dof),
+               assemble_matrix(assemble_matrix),
+               assemble_rhs(assemble_rhs),
+               matrix(matrix),
+               rhs_vector(rhs_vector),
+               quadrature(quadrature),
+               update_flags(update_flags)
+{};
+
+
+template <int dim>
+Assembler<dim>::Assembler (Triangulation<dim>  *tria,
+                          const int            level,
+                          const int            index,
+                          const AssemblerData *local_data) :
+               DoFCellAccessor<dim> (tria,level,index, &local_data->dof),
+               cell_matrix (dof_handler->get_fe().dofs_per_cell),
+               cell_vector (Vector<double>(dof_handler->get_fe().dofs_per_cell)),
+               assemble_matrix (local_data->assemble_matrix),
+               assemble_rhs (local_data->assemble_rhs),
+               matrix(local_data->matrix),
+               rhs_vector(local_data->rhs_vector),
+               fe_values (mapping, dof_handler->get_fe(),
+                          local_data->quadrature,
+                          local_data->update_flags)
+{
+  Assert (!assemble_matrix || (matrix.m() == dof_handler->n_dofs()),
+         ExcInvalidData());
+  Assert (!assemble_matrix || (matrix.n() == dof_handler->n_dofs()),
+         ExcInvalidData());
+  Assert (!assemble_rhs || (rhs_vector.size()==dof_handler->n_dofs()),
+         ExcInvalidData());
+};
+
+
+template <int dim>
+void Assembler<dim>::assemble (const Equation<dim> &equation) {
+                                  // re-init fe values for this cell
+  fe_values.reinit (DoFHandler<dim>::cell_iterator (*this));
+  const unsigned int n_dofs = dof_handler->get_fe().dofs_per_cell;
+
+  if (assemble_matrix)
+    cell_matrix.clear ();
+  if (assemble_rhs)
+    cell_vector.clear ();
+
+
+// fill cell matrix and vector if required
+  DoFHandler<dim>::cell_iterator this_cell (*this);
+  if (assemble_matrix && assemble_rhs) 
+    equation.assemble (cell_matrix, cell_vector, fe_values, this_cell);
+  else
+    if (assemble_matrix)
+      equation.assemble (cell_matrix, fe_values, this_cell);
+    else
+      if (assemble_rhs)
+       equation.assemble (cell_vector, fe_values, this_cell);
+      else
+       Assert (false, ExcNoAssemblingRequired());
+
+
+// get indices of dofs
+  std::vector<unsigned int> dofs (n_dofs);
+  get_dof_indices (dofs);
+
+                                  // one could use the
+                                  // @p{distribute_local_to_global} functions
+                                  // here, but they would require getting the
+                                  // dof indices twice, so we leave it the
+                                  // way it was originally programmed.
+  
+                                  // distribute cell matrix
+  if (assemble_matrix)
+    for (unsigned int i=0; i<n_dofs; ++i)
+      for (unsigned int j=0; j<n_dofs; ++j)
+       matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
+
+                                  // distribute cell vector
+  if (assemble_rhs)
+    for (unsigned int j=0; j<n_dofs; ++j)
+      rhs_vector(dofs[j]) += cell_vector(j);
+};
+
+
+// explicit instantiations
+template class Assembler<deal_II_dimension>;
+
+template class TriaRawIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
+template class TriaIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
+template class TriaActiveIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
+
+
 
-#include <numerics/assembler.h>
 #include <numerics/matrices.h>
 #include <numerics/vectors.h>
 #include <dofs/dof_constraints.h>
@@ -422,5 +1166,7 @@ std::string ProblemBase<dim>::get_solution_name () const
 
 
 
+
+
 #endif
 

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.