Assembling a linear system from a PDE

This recipe shows how to assemble a linear system of equations discretizing a PDE.

In particular, we start from a LocalOperator representing the PDE's bilinear form. Rather than defining one by ourselves, we use one already provided in PDELab. Further, we assume that a GridFunctionSpace and Constraints have already been set up.

  typedef Dune::PDELab::ConvectionDiffusionFEM<Problem,FEM> LOP;
  LOP lop(problem);

This per-element LocalOperator can now be plugged into a GridOperator, which is capable of assembling a matrix for the entire domain.

  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,NumberType,NumberType,NumberType,CC,CC> GO;
  auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));

Next, we set up our degree of freedom vector.

  typedef Dune::PDELab::Backend::Vector<GFS,NumberType> X;
  X x(gfs,0.0);

and ensure it matches the Dirichlet boundary conditions at constrained degrees of freedom. In addition to specifying Dirichlet constrained degrees of freedom, it also serves as initial guess at unconstrained ones.

  typedef Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem> G;
  G g(grid.leafGridView(),problem);
  Dune::PDELab::interpolate(g,gfs,x);

Now, from the GridOperator, we can assemble our residual vector depending on our initial guess

  X d(gfs,0.0);
  go.residual(x,d);

and print it to console.

  Dune::printvector(std::cout, native(d), "Residual vector", "");

Likewise, we can assemble our matrix which again depends on the initial guess

  typedef GO::Jacobian M;
  M A(go);
  go.jacobian(x,A);

and print it as well.

  Dune::printmatrix(std::cout, native(A), "Stiffness matrix", "");

Full example code: recipe-linear-system-assembly.cc