*/
namespace Laplace
{
+/**
+ * Laplacian in weak form, namely on the cell <i>Z</i> the matrix
+ * \f[
+ * \int_Z \nu \nabla u \cdot \nabla v \, dx.
+ * \f]
+ *
+ * The FiniteElement in <tt>fe</tt> may be scalar or vector valued. In
+ * the latter case, the Laplacian is applied to each component
+ * separately.
+ *
+ * @author Guido Kanschat
+ * @date 2008, 2009, 2010
+ */
+ template<int dim>
+ void cell_matrix (
+ FullMatrix<double>& M,
+ const FEValuesBase<dim>& fe,
+ const double factor = 1.)
+ {
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ const unsigned int n_components = fe.get_fe().n_components();
+
+ for (unsigned k=0;k<fe.n_quadrature_points;++k)
+ {
+ const double dx = fe.JxW(k) * factor;
+ for (unsigned i=0;i<n_dofs;++i)
+ {
+ for (unsigned j=0;j<n_dofs;++j)
+ for (unsigned int d=0;d<n_components;++d)
+ M(i,j) += dx *
+ (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
+ }
+ }
+ }
+/**
+ * Weak boundary condition of Nitsche type for the Laplacian, namely on the face <i>F<//i> the matrix
+ * @f[
+ * \int_F \Bigl(\gamma u v - \partial_n u v - u \partial_n v\Bigr)\;ds.
+ * @f]
+ *
+ * Here, γ is the <tt>penalty</tt> parameter suitably computed
+ * with compute_penalty().
+ */
+ template <int dim>
+ void nitsche_matrix (
+ FullMatrix<double>& M,
+ const FEValuesBase<dim>& fe,
+ double penalty,
+ double factor = 1.)
+ {
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ const unsigned int n_comp = fe.get_fe().n_components();
+
+ Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
+ Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
+
+ for (unsigned k=0;k<fe.n_quadrature_points;++k)
+ {
+ const double dx = fe.JxW(k) * factor;
+ const Point<dim>& n = fe.normal_vector(k);
+ for (unsigned i=0;i<n_dofs;++i)
+ for (unsigned j=0;j<n_dofs;++j)
+ for (unsigned int d=0;d<n_comp;++d)
+ M(i,j) += dx *
+ (fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
+ - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
+ - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
+ }
+ }
}
}
* Functions in this namespace follow a generic signature. In the
* simplest case, you have two related functions
* \begin{code}
- * \end{code}
+ * template <int dim>
+ * void
+ * cell_matrix (
+ * FullMatrix<double>& M,
+ * const FEValuesBase<dim>& fe,
+ * const double factor = 1.);
+ *
+ * template <int dim>
+ * void
+ * cell_residual (
+ * BlockVector<double>* v,
+ * const FEValuesBase<dim>& fe,
+ * const std::vector<Tensor<1,dim> >& input,
+ * const double factor = 1.);
+ * \end{code}
+ *
+ * There is typically a pair of functions for the same operator, the
+ * function <tt>cell_residual</tt> implementing the mapping of the
+ * operator from the finite element space into its dual, and the
+ * function <tt>cell_matrix</tt> generating the bilinear form
+ * corresponding to the Frechet derivative of <tt>cell_residual</tt>.
+ *
+ * The first argument of these functions is the return type, which is
+ * <ul>
+ * <li> FullMatrix<double> for matrices
+ * <li> BlockVector<double> for vectors
+ * </li>
+ *
+ * The next argument is the FEValuesBase object representing the
+ * finite element for integration. If the integrated operator maps
+ * from one finite element space into the dual of another (for
+ * instance an off-diagonal matrix in a block system), then first the
+ * FEValuesBase for the trial space and after this the one for the
+ * test space are specified.
+ *
+ * This list is followed by the set of required data in the order
+ * <ol>
+ * <li> Data vectors from finite element functions
+ * <li> Data vectors from other objects
+ * <li> Additional data
+ * <li> A factor which is multiplied with the whole result
+ * </ol>
+ *
+ * <h3>Usage</h3>
+ *
+ * The local integrators can be used wherever a local integration loop
+ * would have been implemented instead. The following example is from
+ * the implementation of a Stokes solver, using
+ * MeshWorker::Assembler::LocalBlocksToGlobalBlocks. The matrices are
+ * <ul>
+ * <li> 0: The vector Laplacian for the velocity (here with a vector valued element)
+ * <li> 1: The divergence matrix
+ * <li> 2: The pressure mass matrix used in the preconditioner
+ * </ul>
+ *
+ * With these matrices, the function called by MeshWorker::loop()
+ * could be written like
+ * \begin{code}
+using namespace ::dealii:: LocalIntegrators;
+
+template <int dim>
+void MatrixIntegrator<dim>::cell(
+ MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationInfo<dim>& info)
+{
+ Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0));
+ Differential::div_primal_matrix(dinfo.matrix(1,false).matrix, info.fe_values(0), info.fe_values(1));
+ L2::cell_matrix(dinfo.matrix(2,false).matrix, info.fe_values(1));
+}
+ * \end{code}
+ * See step-42 for a worked out example of this code.
*/
namespace LocalIntegrators
{