From b9a92342ea7db1d26a3751cfbf31e53b6877f362 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 18 Oct 2010 04:50:43 +0000 Subject: [PATCH] implement first integrators git-svn-id: https://svn.dealii.org/trunk@22355 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/integrators/laplace.h | 69 ++++++++++++++++++ .../include/integrators/local_integrators.h | 72 ++++++++++++++++++- 2 files changed, 140 insertions(+), 1 deletion(-) diff --git a/deal.II/deal.II/include/integrators/laplace.h b/deal.II/deal.II/include/integrators/laplace.h index 1863ab5713..0b05209326 100644 --- a/deal.II/deal.II/include/integrators/laplace.h +++ b/deal.II/deal.II/include/integrators/laplace.h @@ -33,6 +33,75 @@ namespace LocalIntegrators */ namespace Laplace { +/** + * Laplacian in weak form, namely on the cell Z the matrix + * \f[ + * \int_Z \nu \nabla u \cdot \nabla v \, dx. + * \f] + * + * The FiniteElement in fe 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 + void cell_matrix ( + FullMatrix& M, + const FEValuesBase& 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;kF the matrix + * @f[ + * \int_F \Bigl(\gamma u v - \partial_n u v - u \partial_n v\Bigr)\;ds. + * @f] + * + * Here, γ is the penalty parameter suitably computed + * with compute_penalty(). + */ + template + void nitsche_matrix ( + FullMatrix& M, + const FEValuesBase& 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& n = fe.normal_vector(k); + for (unsigned i=0;i + * void + * cell_matrix ( + * FullMatrix& M, + * const FEValuesBase& fe, + * const double factor = 1.); + * + * template + * void + * cell_residual ( + * BlockVector* v, + * const FEValuesBase& fe, + * const std::vector >& input, + * const double factor = 1.); + * \end{code} + * + * There is typically a pair of functions for the same operator, the + * function cell_residual implementing the mapping of the + * operator from the finite element space into its dual, and the + * function cell_matrix generating the bilinear form + * corresponding to the Frechet derivative of cell_residual. + * + * The first argument of these functions is the return type, which is + *
    + *
  • FullMatrix<double> for matrices + *
  • BlockVector<double> for vectors + *
  • + * + * 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 + *
      + *
    1. Data vectors from finite element functions + *
    2. Data vectors from other objects + *
    3. Additional data + *
    4. A factor which is multiplied with the whole result + *
    + * + *

    Usage

    + * + * 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 + *
      + *
    • 0: The vector Laplacian for the velocity (here with a vector valued element) + *
    • 1: The divergence matrix + *
    • 2: The pressure mass matrix used in the preconditioner + *
    + * + * With these matrices, the function called by MeshWorker::loop() + * could be written like + * \begin{code} +using namespace ::dealii:: LocalIntegrators; + +template +void MatrixIntegrator::cell( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationInfo& 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 { -- 2.39.5