+++ /dev/null
-//---------------------------- 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 ---------------------------*/
#include <base/exceptions.h>
#include <map>
+
+// forward declarations
template <int dim> class Quadrature;
template<typename number> class Vector;
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
* 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
* 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
};
+
/**
* Provide a collection of functions operating on matrices. These include
* the application of boundary conditions to a linear system of equations
*/
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);
/**
};
-/**
- * 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
+++ /dev/null
-//---------------------------- 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> >;
+++ /dev/null
-//---------------------------- 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>;
#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>
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));
+ };
};
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);
+ };
};
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));
+ };
};
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
-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
#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>
<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>
#include <fe/fe_lib.lagrange.h>
#include <base/quadrature_lib.h>
#include "../problem_base.h"
-#include <numerics/assembler.h>
#include <lac/sparse_matrix.h>
+/**
+ * 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>
+
+
#endif