From: Wolfgang Bangerth Date: Fri, 15 May 1998 15:59:56 +0000 (+0000) Subject: Add some standard matrices and code to build them. X-Git-Tag: v8.0.0~22966 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=556c63d1fce6ebffbf1e169a8fda826626b0cd17;p=dealii.git Add some standard matrices and code to build them. git-svn-id: https://svn.dealii.org/trunk@288 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h new file mode 100644 index 0000000000..18ca2af6bb --- /dev/null +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -0,0 +1,311 @@ +/*---------------------------- matrices.h ---------------------------*/ +/* $Id$ */ +#ifndef __matrices_H +#define __matrices_H +/*---------------------------- matrices.h ---------------------------*/ + + + +#include + +template class Triangulation; +template class DoFHandler; +template class FiniteElement; +template class FEValues; +template class Quadrature; +template class Function; +template class Boundary; +template class Equation; + +class dVector; +class dFMatrix; +class dSMatrix; + + + +/** + * Provide a class which assembles certain standard matrices for a given + * triangulation, using a given finite element and a quadrature formula. + * All functions are static, so it is not necessary to create an object + * of this type, though you may do so. + * + * + * \subsection{Conventions for all functions} + * + * All functions take a sparse matrix object to hold the matrix to be + * created. The functions assume that the matrix is initialized with a + * sparsity pattern (#dSMatrixStruct#) corresponding to the given degree + * of freedom handler, i.e. the sparsity structure is already as needed. + * You can do this by calling the #DoFHandler::make_sparsity_pattern# + * function. + * + * Furthermore it is assumed that no relevant data is in the matrix. All + * entries will be overwritten. Entries which are not needed by the matrix + * (and were thus added 'by hand' after #make_sparsity_pattern# was called) + * are not touched and in special are not set to zero, so you have to care + * yourself about that if you really need these entries. + * + * + * \subsection{Supported matrices} + * + * At present there are functions to create the following matrices: + * \begin{itemize} + * \item #create_mass_matrix#: create the matrix with entries + * $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$. Uses the + * #MassMatrix# class. + * + * \item #create_laplace_matrix#: there are two versions of this; the + * one which takes the #Function# object creates + * $a_{ij} = \int_\Omega a(x) \nabla\phi_i(x) \nabla\phi_j(x) dx$, + * $a$ being the given function, while the other one assumes that + * $a=1$ which enables some optimzations. In fact the two versions + * are in one function, the coefficient being given as a defaulted + * argument, which is a pointer to a function and defaults to zero. + * This function uses the #LaplaceMatrix# class. + * \end{itemize} + */ +template +class MatrixCreator { + public: + /** + * Assemble the mass matrix. See + * the general doc of this class + * for more information. + */ + static void create_mass_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix); + + /** + * Assemble the laplace matrix with a + * variable weight function. If no + * coefficient is given, it is assumed + * to be zero. + * + * See the general doc of this class + * for more information. + */ + static void create_laplace_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const Function *a = 0); +}; + + + + + +/** + * Equation class to be passed to the #Assembler# if you want to make up the + * mass matrix for your problem. The mass matrix is the matrix with + * $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 + * #NULL# pointer. If you need a coefficient but no right hand side object, + * simply pass a #NULL# pointer to the constructor for its first argument. + */ +template +class MassMatrix : public Equation { + 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 * const rhs = 0, + const Function * 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 (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) const; + + /** + * Construct the cell matrix for this cell. + * If a coefficient was given to the + * constructor, it is used. + */ + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) 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 (dVector &rhs, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) 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 * 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 * const coefficient; +}; + + + + + +/** + * Equation class to be passed to the #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 + * #NULL# pointer. If you need a coefficient but no right hand side object, + * simply pass a #NULL# pointer to the constructor for its first argument. + */ +template +class LaplaceMatrix : public Equation { + 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 * const rhs = 0, + const Function * 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 (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) const; + + /** + * Construct the cell matrix for this cell. + * If a coefficient was given to the + * constructor, it is used. + */ + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) 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 (dVector &rhs, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &cell) 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 * 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 * const coefficient; +}; + + + + + +/*---------------------------- matrices.h ---------------------------*/ +/* end of #ifndef __matrices_H */ +#endif +/*---------------------------- matrices.h ---------------------------*/ diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc new file mode 100644 index 0000000000..f4b068a1bd --- /dev/null +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -0,0 +1,307 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +template +void MatrixCreator::create_mass_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix) { + dVector dummy; // no entries, should give an error if accessed + const AssemblerData data (dof, + true, false, // assemble matrix but not rhs + matrix, dummy, + q, fe, + UpdateFlags(update_jacobians | + update_JxW_values), + boundary); + TriaActiveIterator > + assembler (const_cast*>(&dof.get_tria()), + dof.get_tria().begin_active()->level(), + dof.get_tria().begin_active()->index(), + &data); + MassMatrix equation; + do + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); +}; + + + + +template +void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const Function * const a) { + dVector dummy; // no entries, should give an error if accessed + const AssemblerData data (dof, + true, false, // assemble matrix but not rhs + matrix, dummy, + q, fe, + UpdateFlags(update_gradients | + update_jacobians | + update_JxW_values), + boundary); + TriaActiveIterator > + assembler (const_cast*>(&dof.get_tria()), + dof.get_tria().begin_active()->level(), + dof.get_tria().begin_active()->index(), + &data); + LaplaceMatrix equation (0, a); + do + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); +}; + + + + + + + + +template +MassMatrix::MassMatrix (const Function * const rhs, + const Function * const a) : + Equation (1), + right_hand_side (rhs), + coefficient (a) {}; + + + +template +void MassMatrix::assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const typename Triangulation::cell_iterator &) const { + const dFMatrix &values = fe_values.get_shape_values (); + const vector &weights = fe_values.get_JxW_values (); + + if (coefficient != 0) + { + vector coefficient_values; + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void MassMatrix::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &) const { + Assert (right_hand_side != 0, ExcNoRHSSelected()); + + const dFMatrix &values = fe_values.get_shape_values (); + const vector &weights = fe_values.get_JxW_values (); + vector rhs_values; + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + if (coefficient != 0) + { + vector coefficient_values; + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void MassMatrix::assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &) const { + Assert (right_hand_side != 0, ExcNoRHSSelected()); + + const dFMatrix &values = fe_values.get_shape_values (); + const vector &weights = fe_values.get_JxW_values (); + vector rhs_values; + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + for (unsigned int point=0; point +LaplaceMatrix::LaplaceMatrix (const Function * const rhs, + const Function * const a) : + Equation (1), + right_hand_side (rhs), + coefficient (a) {}; + + +template +void LaplaceMatrix::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &) const { + Assert (right_hand_side != 0, ExcNoRHSSelected()); + + const vector > >&gradients = fe_values.get_shape_grads (); + const dFMatrix &values = fe_values.get_shape_values (); + vector rhs_values; + const vector &weights = fe_values.get_JxW_values (); + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + if (coefficient != 0) + { + vector coefficient_values; + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void LaplaceMatrix::assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &) const { + const vector > >&gradients = fe_values.get_shape_grads (); + const vector &weights = fe_values.get_JxW_values (); + + if (coefficient != 0) + { + vector coefficient_values; + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void LaplaceMatrix::assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &) const { + Assert (right_hand_side != 0, ExcNoRHSSelected()); + + const dFMatrix &values = fe_values.get_shape_values (); + const vector &weights = fe_values.get_JxW_values (); + vector rhs_values; + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + for (unsigned int point=0; point; +template class MatrixCreator<2>; +template class MassMatrix<1>; +template class MassMatrix<2>; +template class LaplaceMatrix<1>; +template class LaplaceMatrix<2>; +