]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Extend each function of MatrixCreator to use an arbitrary mapping. Keep second versio...
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Thu, 5 Apr 2001 11:26:06 +0000 (11:26 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Thu, 5 Apr 2001 11:26:06 +0000 (11:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@4373 0785d39b-7218-0410-832d-ea1e28bc413d

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

index c3919beacf0cc83cd46289bc996208428f2df76c..f9e1200666772af692cf03bb785cc46ea95fcd1e 100644 (file)
@@ -24,6 +24,7 @@ template<typename number> class SparseMatrix;
 template <typename number> class BlockSparseMatrix;
 template <typename Number> class BlockVector;
 
+template <int dim> class Mapping;
 template <int dim> class DoFHandler;
 template <int dim> class MGDoFHandler;
 template <int dim> class FEValues;
@@ -31,14 +32,30 @@ template <int dim> class FEValues;
 
 
 /**
- * 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.
+ * Provide a class which assembles certain standard matrices for a
+ * given triangulation, using a given finite element, a given mapping
+ * 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.
  *
  *
  * @sect3{Conventions for all functions}
  *
+ * There exist two versions of each function. One with a @ref{Mapping}
+ * argument and one without. If a code uses a mapping different from
+ * @ref{MappingQ1} the functions @em{with} mapping argument should be
+ * used. Code that uses only @ref{MappingQ1} may also use the
+ * functions @em{without} @ref{Mapping} argument. Each of these latter
+ * functions create a @ref{MappingQ1} object and just call the
+ * respective functions with that object as mapping argument. The
+ * functions without @ref{Mapping} argument still exist to ensure
+ * backward compatibility. Nevertheless it is advised to change the
+ * user's codes to store a specific @ref{Mapping} object and to use
+ * the functions that take this @p{Mapping} object as argument. This
+ * gives the possibility to easily extend the user codes to work also
+ * on mappings of higher degree, this just by exchanging
+ * @ref{MappingQ1} by, for example, a @ref{MappingQ} or another
+ * @ref{Mapping} object of interest.
+ *
  * 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 (@ref{SparsityPattern}) corresponding to the given degree
@@ -48,9 +65,24 @@ template <int dim> class FEValues;
  *
  * 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 @p{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.
+ * are not touched and in special are not set to zero.
+ * In all cases, the elements of the matrix to be assembled are simply
+ * summed up from the contributions of each cell. Therefore you may want
+ * to clear the matrix before assemblage.
+ *
+ * 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 several matrices and could then condense afterwards
+ * only once, instead of for every matrix. To actually do computations
+ * with these matrices, you have to condense the matrix using the
+ * @ref{ConstraintMatrix}@p{::condense} function; you also have to
+ * condense the right hand side accordingly and distribute the
+ * solution afterwards.
+ *
+ * If you want to use boundary conditions with the matrices generated
+ * by the functions of this class, you have to use a function like
+ * @p{ProblemBase<>::apply_dirichlet_bc} to matrix and right hand
+ * side.
  *
  *
  * @sect3{Supported matrices}
@@ -58,88 +90,56 @@ template <int dim> class FEValues;
  * At present there are functions to create the following matrices:
  * @begin{itemize}
  * @item @p{create_mass_matrix}: create the matrix with entries
- *   $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 @ref{MassMatrix} class.
+ *   $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$ by numerical
+ *   quadrature. Here, the $\phi_i$ are the basis functions of the
+ *   finite element space given.
  *
- *   Two ways to create this matrix are offered. The first one uses
- *   numerical quadrature and the @ref{MassMatrix} class. In this case,
- *   a coefficient may be given to evaluate
+ *   A coefficient may be given to evaluate
  *   $m_{ij} = \int_\Omega a(x) \phi_i(x) \phi_j(x) dx$ instead.
- *   This way of setting up the mass matrix is quite general, but has
- *   some drawbacks, see the documentation of the @ref{MassMatrix} class.
- *
- *   The other way uses exact integration, as offered by the finite
- *   element class used. This way you can avoid quadrature errors and
- *   the assemblage is much faster. However, no coefficient can be
- *   given.
- *
- *   Note that the effect of the two ways of setting up the mass
- *   matrix is not the same if you use finite elements which are
- *   composed of several subelements. In this case, using the
- *   quadrature free way (without coefficient) results in a matrix
- *   which does not couple the subelements, as described in the
- *   @ref{FESystem}@p{::get_local_mass_matrix} documentation, while the
- *   way using quadrature sets up the full matrix, i.e. with the
- *   cross coupling of shape functions belonging to different subelements.
- *
- *   If the finite element for which the mass matrix is to be built
- *   has more than one component, the resulting matrix will not couple
- *   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 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
- *   $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 optimizations. 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 @ref{LaplaceMatrix} class.
- *
- *   If the finite element in use presently has more than only one
- *   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
- * several matrices and could then condense afterwards only once,
- * instead of for every matrix. To actually do computations with these
- * matrices, you have to condense the matrix using the
- * @ref{ConstraintMatrix}@p{::condense} function; you also have to condense the
- * right hand side accordingly and distribute the solution afterwards.
  *
- * In all cases, the elements of the matrix to be assembled are simply
- * summed up from the contributions of each cell. Therefore you may want
- * to clear the matrix before assemblage.
+ * @item @p{create_laplace_matrix}: create the matrix with entries
+ *   $m_{ij} = \int_\Omega \nabla\phi_i(x) \nabla\phi_j(x) dx$ by
+ *   numerical quadrature.
  *
- * If you want to use boundary conditions, you have to use a function
- * like @p{ProblemBase<>::apply_dirichlet_bc} to matrix and right hand
- * side.
+ *   Again, a coefficient may be given to evaluate
+ *   $m_{ij} = \int_\Omega a(x) \nabla\phi_i(x) \phi_j(x) dx$ instead.
+ * @end{itemize}
+ *
+ * Make sure that the order of the @ref{Quadrature} formula given to these
+ * functions is sufficiently high to compute the matrices with the
+ * required accuracy. For the choice of this quadrature rule you need
+ * to take into account the polynomial degree of the @ref{FiniteElement}
+ * basis functions, the roughness of the coefficient @p{a}, as well as
+ * the degree of the given @p{Mapping}.
+ *
+ * Note, that for system elements the mass matrix and the laplace
+ * matrix is implemented such that each components couples only with
+ * itself. I.e. there is no coupling of shape functions belonging to
+ * different components.
+ *
+ * If the finite element for which the mass matrix or the laplace
+ * matrix is to be built has more than one component, this function
+ * accepts a single coefficient as well as a vector valued coefficient
+ * function. For the latter case make sure that the number of
+ * components coincides with the number of components of the system
+ * finite element.
  *
  *
  * @sect3{Matrices on the boundary}
  *
  * The @p{create_boundary_mass_matrix} creates the matrix with entries
- * $m_{ij} = \int_{\Gamma} \phi_i \phi_j dx$, where $\Gamma$ is the union
- * of boundary parts with indicators contained in a set passed to the
- * function (i.e. if you want to set up the mass matrix for the parts of
- * the boundary with indicators zero and 2, you pass the function a set
- * of @p{unsigned char}s as parameter @p{boundary_parts} containing the elements
- * zero and 2). The $\phi_i$ are the basis functions which have at least
- * part of their support om $\Gamma$. The mapping between row and column
- * indices in the mass matrix and the right hand side and the global degree
- * of freedom numbers of the respective basis functions on the whole domain
- * is returned as a vector of numbers which has the same size as the dimension
- * of matrix and right hand side.
+ * $m_{ij} = \int_{\Gamma} \phi_i \phi_j dx$, where $\Gamma$ is the
+ * union of boundary parts with indicators contained in a
+ * @p{FunctionMap} passed to the function (i.e. if you want to set up
+ * the mass matrix for the parts of the boundary with indicators zero
+ * and 2, you pass the function a map of @p{unsigned char}s as
+ * parameter @p{boundary_functions} containing the keys zero and
+ * 2). The $\phi_i$ are the basis functions which have at least part
+ * of their support on $\Gamma$. The mapping
+ * @p{dof_to_boundary_mapping} required by this function maps global
+ * DoF numbers to a numbering of the degrees of freedom located on the
+ * boundary, and can be obtained using the function
+ * @p{DoFTools::map_dof_to_boundary_indices}.
  *
  * Since in most cases we are not interested in the pure mass matrix on the
  * boundary, but rather need it to compute the projection of a function to
@@ -153,35 +153,33 @@ template <int dim> class FEValues;
  * The object describing the exact form of the boundary is obtained from the
  * triangulation object.
  *
+ *
  * @sect3{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. (The @p{create_mass_matrix} function
- * which does not use quadrature does not offer a version to evaluate a right
- * hand side also, since this needs quadrature anyway. Take look at the
- * @ref{VectorTools} class to find a function to set up a right hand side vector
- * only.)
- *
- * 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 @p{VectorCreator} class. The use of the latter may be
- * useful if you want to create many right hand side vectors.
- *
- *
- * All functions in this collection use the finite elemen given to the
- * @ref{DoFHandler} object the last time that the degrees of freedom were
- * distributed on the triangulation.
- *
- * @author Wolfgang Bangerth, 1998
+ * building the right hand side vector. If you want to create a right
+ * hand side vector without creating a matrix, you can use the
+ * @ref{VectorTools::create_right_hand_side} function. The use of the
+ * latter may be useful if you want to create many right hand side
+ * vectors.
+ *
+ * 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.
+ *
+ * All functions in this collection use the finite element given to
+ * the @ref{DoFHandler} object the last time that the degrees of
+ * freedom were distributed on the triangulation.
+ *
+ * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
  */
 template <int dim>
 class MatrixCreator
@@ -209,6 +207,17 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
+    static void create_mass_matrix (const Mapping<dim>       &mapping,
+                                   const DoFHandler<dim>    &dof,
+                                   const Quadrature<dim>    &q,
+                                   SparseMatrix<double>     &matrix,
+                                   const Function<dim>      *a = 0);
+
+                                    /**
+                                     * Calls the @p{create_mass_matrix}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
     static void create_mass_matrix (const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
@@ -223,13 +232,26 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
-    static void create_mass_matrix (const DoFHandler<dim>    &dof,
+    static void create_mass_matrix (const Mapping<dim>       &mapping,
+                                   const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
                                    const Function<dim>      &rhs,
                                    Vector<double>           &rhs_vector,
                                    const Function<dim>      *a = 0);
 
+                                    /**
+                                     * Calls the @p{create_mass_matrix}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
+    static void create_mass_matrix (const DoFHandler<dim>    &dof,
+                                   const Quadrature<dim>    &q,
+                                   SparseMatrix<double>     &matrix,
+                                   const Function<dim>      &rhs,
+                                   Vector<double>           &rhs_vector,
+                                   const Function<dim>      *a = 0);
+    
                                     /**
                                      * Assemble the mass matrix and a right
                                      * hand side vector along the boundary.
@@ -245,12 +267,26 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
+    static void create_boundary_mass_matrix (const Mapping<dim>       &mapping,
+                                            const DoFHandler<dim>    &dof,
+                                            const Quadrature<dim-1>  &q,
+                                            SparseMatrix<double>     &matrix,
+                                            const FunctionMap        &boundary_functions,
+                                            Vector<double>           &rhs_vector,
+                                            std::vector<unsigned int>&dof_to_boundary_mapping,
+                                            const Function<dim>      *a = 0);
+
+                                    /**
+                                     * Calls the @p{create_boundary_mass_matrix}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
     static void create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
                                             const Quadrature<dim-1>  &q,
                                             SparseMatrix<double>     &matrix,
-                                            const FunctionMap        &rhs,
+                                            const FunctionMap        &boundary_functions,
                                             Vector<double>           &rhs_vector,
-                                            std::vector<unsigned int>&vec_to_dof_mapping,
+                                            std::vector<unsigned int>&dof_to_boundary_mapping,
                                             const Function<dim>      *a = 0);
 
                                     /**
@@ -261,6 +297,17 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
+    static void create_laplace_matrix (const Mapping<dim>       &mapping,
+                                      const DoFHandler<dim>    &dof,
+                                      const Quadrature<dim>    &q,
+                                      SparseMatrix<double>     &matrix,
+                                      const Function<dim>      *a = 0);
+    
+                                    /**
+                                     * Calls the @p{create_laplace_matrix}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
     static void create_laplace_matrix (const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
                                       SparseMatrix<double>     &matrix,
@@ -287,6 +334,19 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
+    static void create_laplace_matrix (const Mapping<dim>       &mapping,
+                                      const DoFHandler<dim>    &dof,
+                                      const Quadrature<dim>    &q,
+                                      SparseMatrix<double>     &matrix,
+                                      const Function<dim>      &rhs,
+                                      Vector<double>           &rhs_vector,
+                                      const Function<dim>      *a = 0);
+
+                                    /**
+                                     * Calls the @p{create_laplace_matrix}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
     static void create_laplace_matrix (const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
                                       SparseMatrix<double>     &matrix,
index 3c81a0d92f31a421aac11c601f5285afced13613..f286436a7f37a4efd510e92fbad397adeeceb4a3 100644 (file)
@@ -44,11 +44,11 @@ using namespace std;
 static const MappingQ1<deal_II_dimension> mapping_q1;
 
 
-
-//TODO:[RH,GK] maybe re-create the create_mass_matrix function with 2 args
-
+// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
+// TODO:[WB] implement multithreading for this function.
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
+void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
+                                            const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
                                             const Function<dim> * const coefficient)
@@ -57,7 +57,7 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
   if (coefficient != 0)
     update_flags = UpdateFlags (update_flags | update_q_points);
 
-  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+  FEValues<dim> fe_values (mapping, 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;
@@ -116,9 +116,22 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 };
 
 
-
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
+                                            const Quadrature<dim>    &q,
+                                            SparseMatrix<double>     &matrix,
+                                            const Function<dim> * const coefficient)
+{
+  static const MappingQ1<dim> mapping;
+  create_mass_matrix(mapping, dof, q, matrix, coefficient);
+}
+
+
+// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
+// TODO:[WB] implement multithreading for this function.
+template <int dim>
+void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
+                                            const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
                                             const Function<dim>      &rhs,
@@ -131,7 +144,7 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
   if (coefficient != 0)
     update_flags = UpdateFlags (update_flags | update_q_points);
 
-  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+  FEValues<dim> fe_values (mapping, 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;
@@ -206,6 +219,18 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 };
 
 
+template <int dim>
+void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
+                                            const Quadrature<dim>    &q,
+                                            SparseMatrix<double>     &matrix,
+                                            const Function<dim>      &rhs,
+                                            Vector<double>           &rhs_vector,
+                                            const Function<dim> * const coefficient)
+{
+  static const MappingQ1<dim> mapping;
+  create_mass_matrix(mapping, dof, q, matrix, rhs, rhs_vector, coefficient);
+}
+
 
 #if deal_II_dimension == 1
 
@@ -221,14 +246,32 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>       &,
   Assert (false, ExcNotImplemented());
 };
 
+
+template <>
+void MatrixCreator<1>::create_boundary_mass_matrix (const Mapping<1>          &,
+                                                   const DoFHandler<1>       &,
+                                                   const Quadrature<0>       &,
+                                                   SparseMatrix<double>      &,
+                                                   const FunctionMap         &,
+                                                   Vector<double>            &,
+                                                   std::vector<unsigned int> &,
+                                                   const Function<1>         *)
+{
+  Assert (false, ExcNotImplemented());
+};
+
+
 #endif
 
 
+// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
+// TODO:[WB] implement multithreading for this function.
 template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>     &dof,
+void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapping,
+                                                     const DoFHandler<dim>     &dof,
                                                      const Quadrature<dim-1>   &q,
                                                      SparseMatrix<double>      &matrix,
-                                                     const FunctionMap         &rhs,
+                                                     const FunctionMap         &boundary_functions,
                                                      Vector<double>            &rhs_vector,
                                                      std::vector<unsigned int> &dof_to_boundary_mapping,
                                                      const Function<dim>       *a)
@@ -237,14 +280,14 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
   const unsigned int n_components  = fe.n_components();
   const bool         fe_is_system  = (n_components != 1);
   
-  Assert (matrix.n() == dof.n_boundary_dofs(rhs), ExcInternalError());
+  Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), ExcInternalError());
   Assert (matrix.n() == matrix.m(), ExcInternalError());
   Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
-  Assert (rhs.size() != 0, ExcInternalError());
+  Assert (boundary_functions.size() != 0, ExcInternalError());
   Assert (dof.get_fe() == fe, ExcInternalError());
   Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
          ExcInternalError());
-  Assert (n_components == rhs.begin()->second->n_components,
+  Assert (n_components == boundary_functions.begin()->second->n_components,
          ExcComponentMismatch());
 #ifdef DEBUG
   if (true)
@@ -269,7 +312,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
   UpdateFlags update_flags = UpdateFlags (update_values     |
                                          update_JxW_values |
                                          update_q_points);
-  FEFaceValues<dim> fe_values (mapping_q1, fe, q, update_flags);
+  FEFaceValues<dim> fe_values (mapping, fe, q, update_flags);
 
                                   // two variables for the coefficient,
                                   // one for the two cases indicated in
@@ -292,7 +335,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
     for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                                       // check if this face is on that part of
                                       // the boundary we are interested in
-      if (rhs.find(cell->face(face)->boundary_indicator()) != rhs.end())
+      if (boundary_functions.find(cell->face(face)->boundary_indicator()) != boundary_functions.end())
        {
          cell_matrix.clear ();
          cell_vector.clear ();
@@ -305,7 +348,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
          if (fe_is_system)
                                             // FE has several components
            {
-             rhs.find(cell->face(face)->boundary_indicator())
+             boundary_functions.find(cell->face(face)->boundary_indicator())
                ->second->vector_value_list (fe_values.get_quadrature_points(),
                                             rhs_values_system);
 
@@ -356,7 +399,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
          else
                                             // FE is a scalar one
            {
-             rhs.find(cell->face(face)->boundary_indicator())
+             boundary_functions.find(cell->face(face)->boundary_indicator())
                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
 
              if (a != 0)
@@ -473,7 +516,11 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
            if (fabs(cell_matrix(i,i)) > max_diag_entry)
              max_diag_entry = fabs(cell_matrix(i,i));
 #endif  
-         
+//TODO: [WB] check this place for more efficient alternatives
+//in the innermost loop, we traverse a set twice in "find", but
+//this could be made much faster, e.g. by first building a vector<bool>
+//that already stores the result of find beforehand
+// (then remove dofs_on_face altogether)         
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            for (unsigned int j=0; j<dofs_per_cell; ++j)
              if ((dofs_on_face.find(dofs[i]) != dofs_on_face.end()) &&
@@ -503,9 +550,26 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
 };
 
 
+template <int dim>
+void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>     &dof,
+                                                     const Quadrature<dim-1>   &q,
+                                                     SparseMatrix<double>      &matrix,
+                                                     const FunctionMap         &rhs,
+                                                     Vector<double>            &rhs_vector,
+                                                     std::vector<unsigned int> &dof_to_boundary_mapping,
+                                                     const Function<dim>       *a)
+{
+  static const MappingQ1<dim> mapping;
+  create_boundary_mass_matrix(mapping, dof, q, matrix, rhs, rhs_vector,
+                             dof_to_boundary_mapping, a);
+}
+
 
+// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
+// TODO:[WB] implement multithreading for this function.
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
+void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
+                                               const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
                                                SparseMatrix<double>     &matrix,
                                                const Function<dim> * const coefficient)
@@ -515,7 +579,7 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
   if (coefficient != 0)
     update_flags = UpdateFlags (update_flags | update_q_points);
 
-  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+  FEValues<dim> fe_values (mapping, 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;
@@ -576,6 +640,16 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
+template <int dim>
+void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
+                                               const Quadrature<dim>    &q,
+                                               SparseMatrix<double>     &matrix,
+                                               const Function<dim> * const coefficient)
+{
+  static const MappingQ1<dim> mapping;
+  create_laplace_matrix(mapping, dof, q, matrix, coefficient);
+}
+
 
 //TODO:[GK,RH] maybe recreate this function
 /*
@@ -614,10 +688,11 @@ void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
 
 
 
-
-
+// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
+// TODO:[WB] implement multithreading for this function.
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
+void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
+                                               const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
                                                SparseMatrix<double>     &matrix,
                                                const Function<dim>      &rhs,
@@ -631,7 +706,7 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
   if (coefficient != 0)
     update_flags = UpdateFlags (update_flags | update_q_points);
 
-  FEValues<dim> fe_values (dof.get_fe(), q, update_flags);
+  FEValues<dim> fe_values (mapping, 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;
@@ -709,6 +784,18 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
+template <int dim>
+void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
+                                               const Quadrature<dim>    &q,
+                                               SparseMatrix<double>     &matrix,
+                                               const Function<dim>      &rhs,
+                                               Vector<double>           &rhs_vector,
+                                               const Function<dim> * const coefficient)
+{
+  static const MappingQ1<dim> mapping;  
+  create_laplace_matrix(mapping, dof, q, matrix, rhs, rhs_vector, coefficient);
+}
+
 
 template <int dim>
 template <typename number>

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.