From 256989d9fe22d0812ff933b73c1c503021ba03d0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 27 Mar 2001 15:23:23 +0000 Subject: [PATCH] Purge more classes: Assembler, MassMatrix, LaplaceMatrix, Equation. big-tests do not compile at present, will look into that tomorrow. git-svn-id: https://svn.dealii.org/trunk@4309 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/assembler.h | 311 -------- deal.II/deal.II/include/numerics/matrices.h | 277 +------ deal.II/deal.II/source/numerics/assembler.cc | 132 ---- deal.II/deal.II/source/numerics/equation.cc | 56 -- deal.II/deal.II/source/numerics/matrices.cc | 719 +++++++----------- deal.II/deal.II/source/numerics/vectors.cc | 1 - deal.II/doc/news/2001/c-3-1.html | 7 +- tests/big-tests/poisson/poisson.h | 1 - tests/big-tests/problem_base.h | 748 ++++++++++++++++++- 9 files changed, 1044 insertions(+), 1208 deletions(-) delete mode 100644 deal.II/deal.II/include/numerics/assembler.h delete mode 100644 deal.II/deal.II/source/numerics/assembler.cc delete mode 100644 deal.II/deal.II/source/numerics/equation.cc diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h deleted file mode 100644 index 4c58e08aa3..0000000000 --- a/deal.II/deal.II/include/numerics/assembler.h +++ /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 -#include -#include -#include -#include -#include -#include - - -/** - * 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 -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 &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::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 &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 -class Assembler : public DoFCellAccessor -{ - 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 &dof, - const bool assemble_matrix, - const bool assemble_rhs, - SparseMatrix &matrix, - Vector &rhs_vector, - const Quadrature &quadrature, - const UpdateFlags &update_flags); - - /** - * Pointer to the dof handler object - * to be used to iterate on. - */ - const DoFHandler &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 &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 &rhs_vector; - - /** - * Pointer to a quadrature object to be - * used for this assemblage process. - */ - const Quadrature &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 *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 &); - - /** - * Exception. - */ - DeclException0 (ExcNoAssemblingRequired); - /** - * Exception. - */ - DeclException0 (ExcInvalidData); - /** - * Exception. - */ - /** - * Exception. - */ - private: - /** - * Store a local cell matrix. - */ - FullMatrix cell_matrix; - - /** - * Right hand side local to cell. - */ - Vector 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 &matrix; - - /** - * Pointer to the vector to be assembled - * by this object. - */ - Vector &rhs_vector; - - /** - * The finite element evaluated at the - * quadrature points. - */ - FEValues fe_values; -}; - - -/*---------------------------- problem_assembler.h ---------------------------*/ - -#endif -/*---------------------------- problem_assembler.h ---------------------------*/ diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 32e186486d..c3919beacf 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -13,6 +13,8 @@ #include #include + +// forward declarations template class Quadrature; template class Vector; @@ -25,7 +27,8 @@ template class BlockVector; template class DoFHandler; template class MGDoFHandler; template class FEValues; -template class Equation; + + /** * Provide a class which assembles certain standard matrices for a given @@ -85,8 +88,11 @@ template 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 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 */ static void apply_boundary_values (const std::map &boundary_values, - BlockSparseMatrix &matrix, - BlockVector &solution, - BlockVector &right_hand_side, + BlockSparseMatrix &matrix, + BlockVector &solution, + BlockVector &right_hand_side, const bool eliminate_columns = true); /** @@ -440,256 +448,5 @@ class MatrixTools : public MatrixCreator }; -/** - * 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 -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. - * - * This function assumes the cell matrix - * and right hand side to have the right - * size and to be empty. - */ - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::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 &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 * 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 @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 -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 (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::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 &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 * 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; -}; - #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 index 11b48d2918..0000000000 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ /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 -#include -#include -#include -#include -#include -#include -#include -#include - - -// 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 mapping; - -template -Assembler::AssemblerData::AssemblerData (const DoFHandler &dof, - const bool assemble_matrix, - const bool assemble_rhs, - SparseMatrix &matrix, - Vector &rhs_vector, - const Quadrature &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 -Assembler::Assembler (Triangulation *tria, - const int level, - const int index, - const AssemblerData *local_data) : - DoFCellAccessor (tria,level,index, &local_data->dof), - cell_matrix (dof_handler->get_fe().dofs_per_cell), - cell_vector (Vector(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 -void Assembler::assemble (const Equation &equation) { - // re-init fe values for this cell - fe_values.reinit (DoFHandler::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::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 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; - -template class TriaRawIterator >; -template class TriaIterator >; -template class TriaActiveIterator >; diff --git a/deal.II/deal.II/source/numerics/equation.cc b/deal.II/deal.II/source/numerics/equation.cc deleted file mode 100644 index af05025e27..0000000000 --- a/deal.II/deal.II/source/numerics/equation.cc +++ /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 -#include -#include -#include -#include -#include -#include - -template -Equation::Equation (const unsigned int n_equations) : - n_eq(n_equations) {}; - - -template -void Equation::assemble (FullMatrix &, - Vector &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - - -template -void Equation::assemble (FullMatrix &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - - -template -void Equation::assemble (Vector &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - - -template class Equation; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index d0f4c1abd8..004bac89b1 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -21,7 +21,6 @@ #include #include #include -#include #include #include #include @@ -52,27 +51,68 @@ template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const a) + const Function * const coefficient) { - Vector 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::AssemblerData data (dof, - true, false, // assemble matrix but not rhs - matrix, dummy, - q, update_flags); - TriaActiveIterator > - assembler (const_cast*>(&dof.get_tria()), - dof.get_tria().begin_active()->level(), - dof.get_tria().begin_active()->index(), - &data); - MassMatrix equation(0,a); - do + + FEValues 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 &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector coefficient_values (n_q_points); + + std::vector dof_indices (dofs_per_cell); + + typename DoFHandler::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 &values = fe_values.get_shape_values (); + const std::vector &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::create_mass_matrix (const DoFHandler &dof, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const a) + const Function * const coefficient) { UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | update_JxW_values); - const Assembler::AssemblerData data (dof, - true, true, - matrix, rhs_vector, - q, update_flags); - TriaActiveIterator > - assembler (const_cast*>(&dof.get_tria()), - dof.get_tria().begin_active()->level(), - dof.get_tria().begin_active()->index(), - &data); - MassMatrix equation(&rhs,a); - do + if (coefficient != 0) + update_flags = UpdateFlags (update_flags | update_q_points); + + FEValues 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 &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + std::vector rhs_values (fe_values.n_quadrature_points); + std::vector coefficient_values (n_q_points); + + std::vector dof_indices (dofs_per_cell); + + typename DoFHandler::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 &values = fe_values.get_shape_values (); + const std::vector &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 void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const a) + const Function * 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 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::AssemblerData data (dof, - true, false, // assemble matrix but not rhs - matrix, dummy, - q, update_flags); - 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 + FEValues 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 &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector coefficient_values (n_q_points); + + std::vector dof_indices (dofs_per_cell); + + typename DoFHandler::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 > > + &grads = fe_values.get_shape_grads (); + const std::vector &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::create_laplace_matrix (const DoFHandler &dof, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const a) + const Function * 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::AssemblerData data (dof, - true, true, - matrix, rhs_vector, - q, update_flags); - TriaActiveIterator > - assembler (const_cast*>(&dof.get_tria()), - dof.get_tria().begin_active()->level(), - dof.get_tria().begin_active()->index(), - &data); - LaplaceMatrix equation (&rhs, a); - do + if (coefficient != 0) + update_flags = UpdateFlags (update_flags | update_q_points); + + FEValues 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 &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + std::vector rhs_values (fe_values.n_quadrature_points); + std::vector coefficient_values (n_q_points); + + std::vector dof_indices (dofs_per_cell); + + typename DoFHandler::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 &values = fe_values.get_shape_values (); + const std::vector > > + &grads = fe_values.get_shape_grads (); + const std::vector &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 template void @@ -996,383 +1196,12 @@ MatrixTools::apply_boundary_values (const std::map &bo -template -MassMatrix::MassMatrix (const Function * const rhs, - const Function * const a) : - Equation (1), - right_hand_side (rhs), - coefficient (a) -{}; - - - -template -void MassMatrix::assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::cell_iterator &) const -{ - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert (cell_matrix.n() == dofs_per_cell, - Equation::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); - Assert (cell_matrix.m() == dofs_per_cell, - Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); - Assert (cell_matrix.all_zero(), - Equation::ExcObjectNotEmpty()); - - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); - - - if (coefficient != 0) - { - if (coefficient->n_components == 1) - // scalar coefficient given - { - std::vector coefficient_values (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int i=0; i > coefficient_values (fe_values.n_quadrature_points, - Vector(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int i=0; i -void MassMatrix::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); - Assert (cell_matrix.m() == dofs_per_cell, - Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); - Assert (rhs.size() == dofs_per_cell, - Equation::ExcWrongSize(rhs.size(), dofs_per_cell)); - Assert (cell_matrix.all_zero(), - Equation::ExcObjectNotEmpty()); - Assert (rhs.all_zero(), - Equation::ExcObjectNotEmpty()); - - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); - std::vector rhs_values (fe_values.n_quadrature_points); - right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - - if (coefficient != 0) - { - std::vector coefficient_values (n_q_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point -void MassMatrix::assemble (Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &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::ExcWrongSize(rhs.size(), dofs_per_cell)); - Assert (rhs.all_zero(), - Equation::ExcObjectNotEmpty()); - - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); - std::vector 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 -LaplaceMatrix::LaplaceMatrix (const Function * const rhs, - const Function * const a) : - Equation (1), - right_hand_side (rhs), - coefficient (a) {}; - - - -template -void LaplaceMatrix::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); - Assert (cell_matrix.m() == dofs_per_cell, - Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); - Assert (rhs.size() == dofs_per_cell, - Equation::ExcWrongSize(rhs.size(), dofs_per_cell)); - Assert (cell_matrix.all_zero(), - Equation::ExcObjectNotEmpty()); - Assert (rhs.all_zero(), - Equation::ExcObjectNotEmpty()); - - const std::vector > >&gradients = fe_values.get_shape_grads (); - const FullMatrix &values = fe_values.get_shape_values (); - std::vector rhs_values(fe_values.n_quadrature_points); - const std::vector &weights = fe_values.get_JxW_values (); - right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - - if (coefficient != 0) - { - std::vector coefficient_values(n_q_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point -void LaplaceMatrix::assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::cell_iterator &) const -{ - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - - const FiniteElement &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); - Assert (cell_matrix.m() == dofs_per_cell, - Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); - Assert (cell_matrix.all_zero(), - Equation::ExcObjectNotEmpty()); - - const std::vector > >&gradients = fe_values.get_shape_grads (); - const std::vector &weights = fe_values.get_JxW_values (); - - if (coefficient != 0) - { - std::vector coefficient_values(n_q_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point -void LaplaceMatrix::assemble (Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::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 &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::ExcWrongSize(rhs.size(), dofs_per_cell)); - Assert (rhs.all_zero(), - Equation::ExcObjectNotEmpty()); - - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); - std::vector 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; template class MatrixTools; -template class MassMatrix; -template class LaplaceMatrix; template diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 95cb3a9f8a..c41ef605e2 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -20,7 +20,6 @@ #include #include #include -#include #include #include #include diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 886e189203..995dc8eefc 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -283,7 +283,12 @@ documentation, etc.
  • Removed: the ProblemBase class, which has been deprecated since before the release of - deal.II 3.0 has finally been removed. + deal.II 3.0 has finally been removed. The + same applied for the classes + Assembler, + Equation, + MassMatrix, and + LaplaceMatrix.
    (WB 2001/03/27)

    diff --git a/tests/big-tests/poisson/poisson.h b/tests/big-tests/poisson/poisson.h index ce4ea0c407..021381e394 100644 --- a/tests/big-tests/poisson/poisson.h +++ b/tests/big-tests/poisson/poisson.h @@ -19,7 +19,6 @@ #include #include #include "../problem_base.h" -#include #include diff --git a/tests/big-tests/problem_base.h b/tests/big-tests/problem_base.h index fa5b32c741..ddcfed4313 100644 --- a/tests/big-tests/problem_base.h +++ b/tests/big-tests/problem_base.h @@ -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 +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. + * + * This function assumes the cell matrix + * and right hand side to have the right + * size and to be empty. + */ + virtual void assemble (FullMatrix &cell_matrix, + Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &cell_matrix, + const FEValues &fe_values, + const typename DoFHandler::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 &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 * 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 @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 +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 (FullMatrix &cell_matrix, + Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &cell_matrix, + const FEValues &fe_values, + const typename DoFHandler::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 &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 * 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; +}; + + + + + +template +MassMatrix::MassMatrix (const Function * const rhs, + const Function * const a) : + Equation (1), + right_hand_side (rhs), + coefficient (a) +{}; + + + +template +void MassMatrix::assemble (FullMatrix &cell_matrix, + const FEValues &fe_values, + const typename DoFHandler::cell_iterator &) const +{ + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + const FiniteElement &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert (cell_matrix.n() == dofs_per_cell, + Equation::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); + Assert (cell_matrix.m() == dofs_per_cell, + Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); + Assert (cell_matrix.all_zero(), + Equation::ExcObjectNotEmpty()); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + + + if (coefficient != 0) + { + if (coefficient->n_components == 1) + // scalar coefficient given + { + std::vector coefficient_values (fe_values.n_quadrature_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int i=0; i > coefficient_values (fe_values.n_quadrature_points, + Vector(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int i=0; i +void MassMatrix::assemble (FullMatrix &cell_matrix, + Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); + Assert (cell_matrix.m() == dofs_per_cell, + Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); + Assert (rhs.size() == dofs_per_cell, + Equation::ExcWrongSize(rhs.size(), dofs_per_cell)); + Assert (cell_matrix.all_zero(), + Equation::ExcObjectNotEmpty()); + Assert (rhs.all_zero(), + Equation::ExcObjectNotEmpty()); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + std::vector rhs_values (fe_values.n_quadrature_points); + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + if (coefficient != 0) + { + std::vector coefficient_values (n_q_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void MassMatrix::assemble (Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &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::ExcWrongSize(rhs.size(), dofs_per_cell)); + Assert (rhs.all_zero(), + Equation::ExcObjectNotEmpty()); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + std::vector 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 +LaplaceMatrix::LaplaceMatrix (const Function * const rhs, + const Function * const a) : + Equation (1), + right_hand_side (rhs), + coefficient (a) {}; + + + +template +void LaplaceMatrix::assemble (FullMatrix &cell_matrix, + Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); + Assert (cell_matrix.m() == dofs_per_cell, + Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); + Assert (rhs.size() == dofs_per_cell, + Equation::ExcWrongSize(rhs.size(), dofs_per_cell)); + Assert (cell_matrix.all_zero(), + Equation::ExcObjectNotEmpty()); + Assert (rhs.all_zero(), + Equation::ExcObjectNotEmpty()); + + const std::vector > >&gradients = fe_values.get_shape_grads (); + const FullMatrix &values = fe_values.get_shape_values (); + std::vector rhs_values(fe_values.n_quadrature_points); + const std::vector &weights = fe_values.get_JxW_values (); + right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); + + if (coefficient != 0) + { + std::vector coefficient_values(n_q_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void LaplaceMatrix::assemble (FullMatrix &cell_matrix, + const FEValues &fe_values, + const typename DoFHandler::cell_iterator &) const +{ + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + + const FiniteElement &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::ExcWrongSize(cell_matrix.n(), dofs_per_cell)); + Assert (cell_matrix.m() == dofs_per_cell, + Equation::ExcWrongSize(cell_matrix.m(), dofs_per_cell)); + Assert (cell_matrix.all_zero(), + Equation::ExcObjectNotEmpty()); + + const std::vector > >&gradients = fe_values.get_shape_grads (); + const std::vector &weights = fe_values.get_JxW_values (); + + if (coefficient != 0) + { + std::vector coefficient_values(n_q_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point +void LaplaceMatrix::assemble (Vector &rhs, + const FEValues &fe_values, + const typename DoFHandler::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 &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::ExcWrongSize(rhs.size(), dofs_per_cell)); + Assert (rhs.all_zero(), + Equation::ExcObjectNotEmpty()); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + std::vector 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 +#include +#include +#include +#include +#include +#include +#include + + +// 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 mapping; + +template +Assembler::AssemblerData::AssemblerData (const DoFHandler &dof, + const bool assemble_matrix, + const bool assemble_rhs, + SparseMatrix &matrix, + Vector &rhs_vector, + const Quadrature &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 +Assembler::Assembler (Triangulation *tria, + const int level, + const int index, + const AssemblerData *local_data) : + DoFCellAccessor (tria,level,index, &local_data->dof), + cell_matrix (dof_handler->get_fe().dofs_per_cell), + cell_vector (Vector(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 +void Assembler::assemble (const Equation &equation) { + // re-init fe values for this cell + fe_values.reinit (DoFHandler::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::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 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; + +template class TriaRawIterator >; +template class TriaIterator >; +template class TriaActiveIterator >; + + -#include #include #include #include @@ -422,5 +1166,7 @@ std::string ProblemBase::get_solution_name () const + + #endif -- 2.39.5