* 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
* 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.
const Quadrature<dim> &q,
const Boundary<dim> &boundary,
dSMatrix &matrix,
+ const Function<dim> &rhs,
+ dVector &rhs_vector,
const Function<dim> *a = 0);
};
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.
*/
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
*/
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
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.
*/
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
*/
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
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);
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(),
+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);
+};
+