]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add some standard matrices and code to build them.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 15 May 1998 15:59:56 +0000 (15:59 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 15 May 1998 15:59:56 +0000 (15:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@288 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/matrices.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h
new file mode 100644 (file)
index 0000000..18ca2af
--- /dev/null
@@ -0,0 +1,311 @@
+/*----------------------------   matrices.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __matrices_H
+#define __matrices_H
+/*----------------------------   matrices.h     ---------------------------*/
+
+
+
+#include <base/exceptions.h>
+
+template <int dim> class Triangulation;
+template <int dim> class DoFHandler;
+template <int dim> class FiniteElement;
+template <int dim> class FEValues;
+template <int dim> class Quadrature;
+template <int dim> class Function;
+template <int dim> class Boundary;
+template <int dim> 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<dim>::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<dim># 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 <int dim>
+class MatrixCreator {
+  public:
+                                    /**
+                                     * Assemble the mass matrix. 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);
+
+                                    /**
+                                     * 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<dim>    &dof,
+                                      const FiniteElement<dim> &fe,
+                                      const Quadrature<dim>    &q,
+                                      const Boundary<dim>      &boundary,
+                                      dSMatrix &matrix,
+                                      const Function<dim>      *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 <int dim>
+class MassMatrix :  public Equation<dim> {
+  public:
+                                    /**
+                                     * Constructor. Pass a function object if
+                                     * you want to create a right hand side
+                                     * vector, pass a function pointer (default
+                                     * is a NULL pointer). It is your duty to
+                                     * guarantee that the function object for
+                                     * the right hand side lives at least as
+                                     * long as this object does.
+                                     *
+                                     * You may also pass a function describing
+                                     * the weight to the integral (see the
+                                     * general docs for more information). The
+                                     * same applies for this object as said
+                                     * above.
+                                     */
+    MassMatrix (const Function<dim> * const rhs = 0,
+               const Function<dim> * const a = 0);
+
+                                    /**
+                                     * Assemble the cell matrix and right hand
+                                     * side vector for this cell. You need to
+                                     * give a right hand side object to the
+                                     * constructor to use this function. If
+                                     * a coefficient was given to the
+                                     * constructor, it is used.
+                                     */
+    virtual void assemble (dFMatrix            &cell_matrix,
+                          dVector             &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> * const right_hand_side;
+
+                                    /**
+                                     * Pointer to a function describing the
+                                     * coefficient to the integral for the
+                                     * matrix entries. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const coefficient;
+};
+
+
+
+
+
+/**
+ * Equation class to be passed to the #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 <int dim>
+class LaplaceMatrix :  public Equation<dim> {
+  public:
+                                    /**
+                                     * Constructor. Pass a function object if
+                                     * you want to create a right hand side
+                                     * vector, pass a function pointer (default
+                                     * is a NULL pointer). It is your duty to
+                                     * guarantee that the function object for
+                                     * the right hand side lives at least as
+                                     * long as this object does.
+                                     *
+                                     * You may also pass a function describing
+                                     * the weight to the integral (see the
+                                     * general docs for more information). The
+                                     * same applies for this object as said
+                                     * above.
+                                     */
+    LaplaceMatrix (const Function<dim> * const rhs = 0,
+                  const Function<dim> * const a = 0);
+
+                                    /**
+                                     * Assemble the cell matrix and right hand
+                                     * side vector for this cell. You need to
+                                     * give a right hand side object to the
+                                     * constructor to use this function. If
+                                     * a coefficient was given to the
+                                     * constructor, it is used.
+                                     */
+    virtual void assemble (dFMatrix            &cell_matrix,
+                          dVector             &rhs,
+                          const FEValues<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> &fe_values,
+                          const typename Triangulation<dim>::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<dim> * const right_hand_side;
+
+                                    /**
+                                     * Pointer to a function describing the
+                                     * coefficient to the integral for the
+                                     * matrix entries. Should
+                                     * be zero if not given to the constructor
+                                     * and should then not be used.
+                                     */
+    const Function<dim> * const coefficient;
+};
+
+
+
+
+
+/*----------------------------   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 (file)
index 0000000..f4b068a
--- /dev/null
@@ -0,0 +1,307 @@
+/*      $Id$                 */
+
+#include <basic/function.h>
+#include <grid/dof.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <fe/quadrature.h>
+#include <fe/fe_values.h>
+#include <numerics/matrices.h>
+#include <numerics/assembler.h>
+#include <lac/dsmatrix.h>
+
+
+
+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) {
+  dVector dummy;    // no entries, should give an error if accessed
+  const AssemblerData<dim> data (dof,
+                                true, false,  // assemble matrix but not rhs
+                                matrix, dummy,
+                                q, fe,
+                                UpdateFlags(update_jacobians |
+                                            update_JxW_values),
+                                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;
+  do 
+    {
+      assembler->assemble (equation);
+    }
+  while ((++assembler).state() == valid);
+};
+
+
+
+
+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> * const a) {
+  dVector dummy;   // no entries, should give an error if accessed
+  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);
+  TriaActiveIterator<dim, Assembler<dim> >
+    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
+              dof.get_tria().begin_active()->level(),
+              dof.get_tria().begin_active()->index(),
+              &data);
+  LaplaceMatrix<dim> equation (0, a);
+  do 
+    {
+      assembler->assemble (equation);
+    }
+  while ((++assembler).state() == valid);
+};
+
+
+
+
+
+
+
+
+template <int dim>
+MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
+                            const Function<dim> * const a) :
+               Equation<dim> (1),
+               right_hand_side (rhs),
+               coefficient (a)   {};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (dFMatrix            &cell_matrix,
+                               const FEValues<dim> &fe_values,
+                               const typename Triangulation<dim>::cell_iterator &) const {
+  const dFMatrix       &values    = fe_values.get_shape_values ();
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+
+  if (coefficient != 0)
+    {
+      vector<double> coefficient_values;
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+         for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+           cell_matrix(i,j) += (values(i,point) *
+                                values(j,point) *
+                                weights[point] *
+                                coefficient_values[point]);
+    }
+  else
+    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+      for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+       for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+         cell_matrix(i,j) += (values(i,point) *
+                              values(j,point) *
+                              weights[point]);
+};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (dFMatrix            &cell_matrix,
+                               dVector             &rhs,
+                               const FEValues<dim> &fe_values,
+                               const Triangulation<dim>::cell_iterator &) const {
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+  
+  const dFMatrix       &values    = fe_values.get_shape_values ();
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  vector<double>        rhs_values;
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  if (coefficient != 0)
+    {
+      vector<double> coefficient_values;
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+         {
+           for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+             cell_matrix(i,j) += (values(i,point) *
+                                  values(j,point) *
+                                  weights[point] *
+                                  coefficient_values[point]);
+           rhs(i) += values(i,point) *
+                     rhs_values[point] *
+                     weights[point];
+         };
+    }
+  else
+    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+      for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+       {
+         for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+           cell_matrix(i,j) += (values(i,point) *
+                                values(j,point) *
+                                weights[point]);
+         rhs(i) += values(i,point) *
+                   rhs_values[point] *
+                   weights[point];
+       };
+};
+
+
+
+template <int dim>
+void MassMatrix<dim>::assemble (dVector             &rhs,
+                               const FEValues<dim> &fe_values,
+                               const Triangulation<dim>::cell_iterator &) const {
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+  
+  const dFMatrix       &values    = fe_values.get_shape_values ();
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  vector<double>        rhs_values;
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+      rhs(i) += values(i,point) *
+               rhs_values[point] *
+               weights[point];
+};
+
+
+
+
+
+template <int dim>
+LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
+                                  const Function<dim> * const a) :
+               Equation<dim> (1),
+               right_hand_side (rhs),
+               coefficient (a) {};
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (dFMatrix            &cell_matrix,
+                                  dVector             &rhs,
+                                  const FEValues<dim> &fe_values,
+                                  const Triangulation<dim>::cell_iterator &) const {
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+  
+  const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+  const dFMatrix       &values    = fe_values.get_shape_values ();
+  vector<double>        rhs_values;
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+
+  if (coefficient != 0)
+    {
+      vector<double> coefficient_values;
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+         {
+           for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+             cell_matrix(i,j) += (gradients[i][point] *
+                                  gradients[j][point]) *
+                                 weights[point] *
+                                 coefficient_values[point];
+           rhs(i) += values(i,point) *
+                     rhs_values[point] *
+                     weights[point];
+         };
+    }
+  else
+    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+      for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+       {
+         for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+           cell_matrix(i,j) += (gradients[i][point] *
+                                gradients[j][point]) *
+                               weights[point];
+         rhs(i) += values(i,point) *
+                   rhs_values[point] *
+                   weights[point];
+       };
+
+};
+
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (dFMatrix            &cell_matrix,
+                                  const FEValues<dim> &fe_values,
+                                  const Triangulation<dim>::cell_iterator &) const {
+  const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+   
+  if (coefficient != 0)
+    {
+      vector<double> coefficient_values;
+      coefficient->value_list (fe_values.get_quadrature_points(),
+                              coefficient_values);
+      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+         for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+           cell_matrix(i,j) += (gradients[i][point] *
+                                gradients[j][point]) *
+                               weights[point] *
+                               coefficient_values[point];
+    }
+  else
+    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+      for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+       for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+         cell_matrix(i,j) += (gradients[i][point] *
+                              gradients[j][point]) *
+                             weights[point];
+};
+
+
+
+template <int dim>
+void LaplaceMatrix<dim>::assemble (dVector             &rhs,
+                                  const FEValues<dim> &fe_values,
+                                  const Triangulation<dim>::cell_iterator &) const {
+  Assert (right_hand_side != 0, ExcNoRHSSelected());
+  
+  const dFMatrix       &values    = fe_values.get_shape_values ();
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  vector<double>        rhs_values;
+  right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values);
+   
+  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+      rhs(i) += values(i,point) *
+               rhs_values[point] *
+               weights[point];
+};
+
+
+
+
+
+
+
+
+template class MatrixCreator<1>;
+template class MatrixCreator<2>;
+template class MassMatrix<1>;
+template class MassMatrix<2>;
+template class LaplaceMatrix<1>;
+template class LaplaceMatrix<2>;
+

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.