How to get started

Complete source code

#include "hp2D.hh"
#include "concepts.hh"

using concepts::Real;

int main(int argc, char** argv) {

  // build a mesh of a unit square with edges of attribute 1
  concepts::Square mesh(1, 1, 1);

  // set homogeneous Dirichlet boundary conditions for every edge of attribute 1 (i.e. all edges) 
  concepts::BoundaryConditions bc;
  bc.add(concepts::Attribute(1),concepts::Boundary(concepts::Boundary::DIRICHLET));

  // build the finite element space V_2_5
  uint levelOfRefinement = 2;
  uint polynomialDegree = 6;
  hp2D::hpAdaptiveSpaceH1 space(mesh,levelOfRefinement,polynomialDegree,&bc);
  space.rebuild();

  // define the linear form f of the right hand side
  concepts::ParsedFormula<> f("((2*pi*pi+1)*sin(pi*x)*sin(pi*y))");
  hp2D::Riesz<Real> lform(f);

  // compute the vector for the finite element discretization
  concepts::Vector<Real> rhs(space, lform);

  // define the bilinear forms
  hp2D::Laplace<Real> la;
  hp2D::Identity<Real> id;

  // compute the mass and stiffness matrices
  concepts::SparseMatrix<Real> S(space, la);
  concepts::SparseMatrix<Real> M(space, id);

  // add mass matrix to stiffness matrix
  M.addInto(S, 1.);

  // compute the solution of the linear problem using the CG-method
  std::unique_ptr<concepts::Operator<Real> > solver(nullptr);
  solver.reset(new concepts::CG<Real>(S, 1e-12, 400));
  concepts::Vector<Real> sol(space);
  (*solver)(rhs, sol);

  // set integration points to be equidistant using trapezoidal rule
  hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, polynomialDegree);
  space.recomputeShapefunctions();

  // save graphical output of the solution in a MAT-file
  graphics::MatlabBinaryGraphics mbg(space, "firstSolution");
  mbg.addSolution(space, "u", sol);

  return 0;
}