How to get started

   1 #include "hp2D.hh"
   2 #include "concepts.hh"
   3 
   4 using concepts::Real;
   5 
   6 int main(int argc, char** argv) {
   7 
   8   // build a mesh of a unit square with edges of attribute 1
   9   concepts::Square mesh(1, 1, 1);
  10 
  11   // set homogeneous Dirichlet boundary conditions for every edge of attribute 1 (i.e. all edges) 
  12   concepts::BoundaryConditions bc;
  13   bc.add(concepts::Attribute(1),concepts::Boundary(concepts::Boundary::DIRICHLET));
  14 
  15   // build the finite element space V_2_5
  16   uint levelOfRefinement = 2;
  17   uint polynomialDegree = 6;
  18   hp2D::hpAdaptiveSpaceH1 space(mesh,levelOfRefinement,polynomialDegree,&bc);
  19   space.rebuild();
  20 
  21   // define the linear form f of the right hand side
  22   concepts::ParsedFormula<> f("((2*pi*pi+1)*sin(pi*x)*sin(pi*y))");
  23   hp2D::Riesz<Real> lform(f);
  24 
  25   // compute the vector for the finite element discretization
  26   concepts::Vector<Real> rhs(space, lform);
  27 
  28   // define the bilinear forms
  29   hp2D::Laplace<Real> la;
  30   hp2D::Identity<Real> id;
  31 
  32   // compute the mass and stiffness matrices
  33   concepts::SparseMatrix<Real> S(space, la);
  34   concepts::SparseMatrix<Real> M(space, id);
  35 
  36   // add mass matrix to stiffness matrix
  37   M.addInto(S, 1.);
  38 
  39   // compute the solution of the linear problem using the CG-method
  40   std::unique_ptr<concepts::Operator<Real> > solver(nullptr);
  41   solver.reset(new concepts::CG<Real>(S, 1e-12, 400));
  42   concepts::Vector<Real> sol(space);
  43   (*solver)(rhs, sol);
  44 
  45   // set integration points to be equidistant using trapezoidal rule
  46   hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, polynomialDegree);
  47   space.recomputeShapefunctions();
  48 
  49   // save graphical output of the solution in a MAT-file
  50   graphics::MatlabBinaryGraphics mbg(space, "firstSolution");
  51   mbg.addSolution(space, "u", sol);
  52 
  53   return 0;
  54 }