]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add creation of right hand sides to the MatrixCreator class.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 15 May 1998 17:11:36 +0000 (17:11 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 15 May 1998 17:11:36 +0000 (17:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@292 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/matrices.cc

index 18ca2af6bbcace42235605f6e21e6d5c11c9220c..bdfe0a2745c3ab14a64e377b1a268a6da197e663 100644 (file)
@@ -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<dim># 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 <int dim>
 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<dim>    &dof,
                                    const FiniteElement<dim> &fe,
                                    const Quadrature<dim>    &q,
                                    const Boundary<dim>      &boundary,
-                                   dSMatrix &matrix);
+                                   dSMatrix                 &matrix,
+                                   const Function<dim>      *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<dim>    &dof,
+                                   const FiniteElement<dim> &fe,
+                                   const Quadrature<dim>    &q,
+                                   const Boundary<dim>      &boundary,
+                                   dSMatrix                 &matrix,
+                                   const Function<dim>      &rhs,
+                                   dVector                  &rhs_vector,
+                                   const Function<dim>      *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<dim>    &dof,
+                                      const FiniteElement<dim> &fe,
+                                      const Quadrature<dim>    &q,
+                                      const Boundary<dim>      &boundary,
+                                      dSMatrix &matrix,
+                                      const Function<dim>      *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<dim>    &q,
                                       const Boundary<dim>      &boundary,
                                       dSMatrix &matrix,
+                                      const Function<dim>      &rhs,
+                                      dVector                  &rhs_vector,
                                       const Function<dim>      *a = 0);
 };
 
@@ -152,7 +212,7 @@ class MassMatrix :  public Equation<dim> {
     virtual void assemble (dFMatrix            &cell_matrix,
                           dVector             &rhs,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
 
                                     /**
                                      * Construct the cell matrix for this cell.
@@ -161,7 +221,7 @@ class MassMatrix :  public Equation<dim> {
                                      */
     virtual void assemble (dFMatrix            &cell_matrix,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
 
                                     /**
                                      * Only construct the right hand side
@@ -172,7 +232,7 @@ class MassMatrix :  public Equation<dim> {
                                      */
     virtual void assemble (dVector             &rhs,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
     
                                     /**
                                      * Exception
@@ -255,7 +315,7 @@ class LaplaceMatrix :  public Equation<dim> {
     virtual void assemble (dFMatrix            &cell_matrix,
                           dVector             &rhs,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
 
                                     /**
                                      * Construct the cell matrix for this cell.
@@ -264,7 +324,7 @@ class LaplaceMatrix :  public Equation<dim> {
                                      */
     virtual void assemble (dFMatrix            &cell_matrix,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
 
                                     /**
                                      * Only construct the right hand side
@@ -275,7 +335,7 @@ class LaplaceMatrix :  public Equation<dim> {
                                      */
     virtual void assemble (dVector             &rhs,
                           const FEValues<dim> &fe_values,
-                          const typename Triangulation<dim>::cell_iterator &cell) const;
+                          const typename Triangulation<dim>::cell_iterator &) const;
 
                                     /**
                                      * Exception
index f4b068a1bd53368697c63186b05ad545dde41ac0..53fc60946d381c4370a30df5b096dd29345b45ed 100644 (file)
@@ -17,21 +17,55 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             const FiniteElement<dim> &fe,
                                             const Quadrature<dim>    &q,
                                             const Boundary<dim>      &boundary,
-                                            dSMatrix &matrix) {
+                                            dSMatrix                 &matrix,
+                                            const Function<dim> * 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<dim> 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<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 
+    {
+      assembler->assemble (equation);
+    }
+  while ((++assembler).state() == valid);
+};
+
+
+
+
+template <int dim>
+void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
+                                            const FiniteElement<dim> &fe,
+                                            const Quadrature<dim>    &q,
+                                            const Boundary<dim>      &boundary,
+                                            dSMatrix                 &matrix,
+                                            const Function<dim>      &rhs,
+                                            dVector                  &rhs_vector,
+                                            const Function<dim> * const a) {
+  UpdateFlags update_flags = UpdateFlags(update_q_points |
+                                        update_jacobians |
+                                        update_JxW_values);
+  const AssemblerData<dim> data (dof,
+                                true, true,
+                                matrix, rhs_vector,
+                                q, fe,  update_flags, boundary);
   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;
+  MassMatrix<dim> equation(&rhs,a);
   do 
     {
       assembler->assemble (equation);
@@ -47,17 +81,18 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const FiniteElement<dim> &fe,
                                                const Quadrature<dim>    &q,
                                                const Boundary<dim>      &boundary,
-                                               dSMatrix &matrix,
+                                               dSMatrix                 &matrix,
                                                const Function<dim> * 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<dim> 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<dim, Assembler<dim> >
     assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
               dof.get_tria().begin_active()->level(),
@@ -73,6 +108,38 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
+template <int dim>
+void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
+                                               const FiniteElement<dim> &fe,
+                                               const Quadrature<dim>    &q,
+                                               const Boundary<dim>      &boundary,
+                                               dSMatrix                 &matrix,
+                                               const Function<dim>      &rhs,
+                                               dVector                  &rhs_vector,
+                                               const Function<dim> * const a) {
+  UpdateFlags update_flags = UpdateFlags(update_q_points  |
+                                        update_gradients |
+                                        update_jacobians |
+                                        update_JxW_values);
+  const AssemblerData<dim> data (dof,
+                                true, true,
+                                matrix, rhs_vector,
+                                q, fe,
+                                update_flags,
+                                boundary);
+  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 
+    {
+      assembler->assemble (equation);
+    }
+  while ((++assembler).state() == valid);
+};
+
 
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.