From: Wolfgang Bangerth Date: Fri, 15 May 1998 17:11:36 +0000 (+0000) Subject: Add creation of right hand sides to the MatrixCreator class. X-Git-Tag: v8.0.0~22962 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a6fa5a7f3d90d2e9f56fbc3e484da3819aef43a8;p=dealii.git Add creation of right hand sides to the MatrixCreator class. git-svn-id: https://svn.dealii.org/trunk@292 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 index 18ca2af6bb..bdfe0a2745 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -51,8 +51,9 @@ class dSMatrix; * 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. + * $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$. Here, the $\phi_i$ + * are the basis functions of the finite element space given. + * This function uses the #MassMatrix# class. * * \item #create_laplace_matrix#: there are two versions of this; the * one which takes the #Function# object creates @@ -63,26 +64,83 @@ class dSMatrix; * argument, which is a pointer to a function and defaults to zero. * This function uses the #LaplaceMatrix# class. * \end{itemize} + * + * + * \subsection{Right hand sides} + * + * In many cases, you will not only want to build the matrix, but also + * a right hand side, which will give a vector with + * $f_i = \int_\Omega f(x) \phi_i(x) dx$. For this purpose, each function + * exists in two versions, one only building the matrix and one also + * building the right hand side vector. + * + * Creation of the right hand side + * is the same for all operators and therefore for all of the functions + * below. It would be most orthogonal to write one single function which + * builds up the right hand side and not provide many functions doing + * the same thing. However, this may result in a heavy performance + * penalty, since then many values of a certain finite element have to + * be computed twice, so it is more economical to implement it more than + * once. If you only want to create a right hand side as above, there is + * a function in the #VectorCreator# class. The use of the latter may be + * useful if you want to create many right hand side vectors. */ template class MatrixCreator { public: /** - * Assemble the mass matrix. See - * the general doc of this class + * Assemble the mass matrix. If no + * coefficient is given, it is assumed + * to be constant one. + * + * 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); + dSMatrix &matrix, + const Function *a = 0); + + /** + * Assemble the mass matrix and a right + * hand side vector. If no + * coefficient is given, it is assumed + * to be constant one. + * + * 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, + const Function &rhs, + dVector &rhs_vector, + const Function *a = 0); + + /** + * Assemble the Laplace matrix. If no + * coefficient is given, it is assumed + * to be constant one. + * + * 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); /** - * Assemble the laplace matrix with a - * variable weight function. If no + * Assemble the Laplace matrix and a right + * hand side vector. If no * coefficient is given, it is assumed - * to be zero. + * to be constant one. * * See the general doc of this class * for more information. @@ -92,6 +150,8 @@ class MatrixCreator { const Quadrature &q, const Boundary &boundary, dSMatrix &matrix, + const Function &rhs, + dVector &rhs_vector, const Function *a = 0); }; @@ -152,7 +212,7 @@ class MassMatrix : public Equation { virtual void assemble (dFMatrix &cell_matrix, dVector &rhs, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Construct the cell matrix for this cell. @@ -161,7 +221,7 @@ class MassMatrix : public Equation { */ virtual void assemble (dFMatrix &cell_matrix, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Only construct the right hand side @@ -172,7 +232,7 @@ class MassMatrix : public Equation { */ virtual void assemble (dVector &rhs, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Exception @@ -255,7 +315,7 @@ class LaplaceMatrix : public Equation { virtual void assemble (dFMatrix &cell_matrix, dVector &rhs, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Construct the cell matrix for this cell. @@ -264,7 +324,7 @@ class LaplaceMatrix : public Equation { */ virtual void assemble (dFMatrix &cell_matrix, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Only construct the right hand side @@ -275,7 +335,7 @@ class LaplaceMatrix : public Equation { */ virtual void assemble (dVector &rhs, const FEValues &fe_values, - const typename Triangulation::cell_iterator &cell) const; + const typename Triangulation::cell_iterator &) const; /** * Exception diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index f4b068a1bd..53fc60946d 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -17,21 +17,55 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const FiniteElement &fe, const Quadrature &q, const Boundary &boundary, - dSMatrix &matrix) { + dSMatrix &matrix, + const Function * const a) { dVector dummy; // no entries, should give an error if accessed + UpdateFlags update_flags = UpdateFlags(update_jacobians | + update_JxW_values); + if (a != 0) + update_flags = UpdateFlags (update_flags | update_q_points); const AssemblerData data (dof, true, false, // assemble matrix but not rhs matrix, dummy, - q, fe, - UpdateFlags(update_jacobians | - update_JxW_values), - boundary); + q, fe, update_flags, boundary); + 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 + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); +}; + + + + +template +void MatrixCreator::create_mass_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const Function &rhs, + dVector &rhs_vector, + const Function * const a) { + UpdateFlags update_flags = UpdateFlags(update_q_points | + update_jacobians | + update_JxW_values); + const AssemblerData data (dof, + true, true, + matrix, rhs_vector, + q, fe, update_flags, boundary); TriaActiveIterator > assembler (const_cast*>(&dof.get_tria()), dof.get_tria().begin_active()->level(), dof.get_tria().begin_active()->index(), &data); - MassMatrix equation; + MassMatrix equation(&rhs,a); do { assembler->assemble (equation); @@ -47,17 +81,18 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const FiniteElement &fe, const Quadrature &q, const Boundary &boundary, - dSMatrix &matrix, + dSMatrix &matrix, const Function * const a) { dVector dummy; // no entries, should give an error if accessed + UpdateFlags update_flags = UpdateFlags(update_gradients | + update_jacobians | + update_JxW_values); + if (a != 0) + update_flags = UpdateFlags(update_flags | update_q_points); const AssemblerData data (dof, true, false, // assemble matrix but not rhs matrix, dummy, - q, fe, - UpdateFlags(update_gradients | - update_jacobians | - update_JxW_values), - boundary); + q, fe, update_flags, boundary); TriaActiveIterator > assembler (const_cast*>(&dof.get_tria()), dof.get_tria().begin_active()->level(), @@ -73,6 +108,38 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, +template +void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const Function &rhs, + dVector &rhs_vector, + const Function * const a) { + UpdateFlags update_flags = UpdateFlags(update_q_points | + update_gradients | + update_jacobians | + update_JxW_values); + const AssemblerData data (dof, + true, true, + matrix, rhs_vector, + q, fe, + update_flags, + boundary); + 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 + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); +}; +