You are here: Concepts>Concepts Web>Class documentation

Class documentation of Concepts

Searching...

No Matches

Introduction

Two dimensional FEM in Concepts. This code features high order Ansatz functions and irregular meshes (local, anisotropic refinements in h and p). All classes for the hp-FEM in two dimensions can be found in the namespace hp2D.

The equation solved is a reaction-diffusion equation (without reaction term)

$- \Delta u = 0$

on the L shaped domain (-1,1)² \ (0,1)x(-1,0) with mixed boundary conditions chosen such that the exact solution is $r^{2/3}\sin(2/3 \phi)$ . The linear system is solved using CG or a direct solver. The boundary conditions are homogeneous Dirichlet on the edges 0 and 9 and non-zero Neumann on the rest of the boundary.

The exact solution is $r^{2/3}\sin(2/3 \phi)$ .

The relative error in energy norm is computed.

Commented Program
Results
Complete Source Code

Commented Program

System include files.

#include <iostream>

Concepts include files. hp2D.hh is for the high order method.

 
#include "basics.hh"
#include "toolbox.hh"
#include "formula.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "operator.hh"
#include "space.hh"
 
#include "hp2D.hh"

With this using directive, we do not need to prepend concepts:: before Real everytime.

using concepts::Real;

concepts::Real

double Real

Definition typedefs.hh:39

Main Program

The mesh is read from three file containing the coordinates, the elements and the boundary conditions of the mesh.

int main(int argc, char** argv) {
 
  try {
    concepts::Import2dMesh msh(SOURCEDIR "/applications/lshape2DCoord.dat",
             SOURCEDIR "/applications/lshape2DElements.dat",
             SOURCEDIR "/applications/lshape2DBoundary.dat");
    std::cout << "msh = " << msh << std::endl;

The boundary conditions are different for each segment of the boundary.

    concepts::BoundaryConditions bc;
    bc.add(1, concepts::Boundary(concepts::Boundary::DIRICHLET));
    bc.add(2, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( (1/(1+y*y)^(2/3)) * (sin(2/3*atan(y)) - y*cos(2/3*atan(y)))*2/3 )")));
    bc.add(3, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*atan(1/x)) + x*cos(2/3*atan(1/x))) )")));
    bc.add(4, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(1/x)+pi)) + x*cos(2/3*(atan(1/x)+pi))) )")));
    bc.add(5, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )")));
    bc.add(6, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )")));
    bc.add(7, concepts::Boundary
     (concepts::Boundary::NEUMANN, concepts::ParsedFormula<>
      ("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(-1/x)+pi)) - x*cos(2/3*(atan(-1/x)+pi))) )")));
    bc.add(8, concepts::Boundary(concepts::Boundary::DIRICHLET));
    std::cout << bc << std::endl;
    std::cout << "--" << std::endl;

Using the mesh and the boundary conditions, the space can be built. Note that the elements are not yet created as this is an adaptive space. It builds its elements only on request via hp2D::Space::rebuild().

hp2D::Space space(msh, 0, 1, &bc);

hp2D::Space

Definition space.hh:88

Now, the loop over the different computational levels is initialized.

for (uint i = 0; i < 3; ++i) {

std::cout << "i = " << i << std::endl;

Now, it is time to build the elements of the space. The newly built space is printed.

space.rebuild();

std::cout << " space = " << space << std::endl;

The right hand side is computed. In our case, the vector rhs only contains the integrals of the Neumann boundary conditions as the right hand side f is 0.

hp2D::Riesz<> lform(concepts::ParsedFormula<>("(0)"), &bc);

concepts::Vector<Real> rhs(space, lform);

concepts::Vector

Definition vector.hh:39

hp2D::Riesz

Definition linearForm.hh:62

The stiffness matrix is computed from the bilinear form hp2D::Laplace.

      hp2D::Laplace<> la;
      concepts::SparseMatrix<Real> stiff(space, la);
      stiff.compress();
      std::cout << "  matrix: " << stiff << std::endl;

The next part is the linear solver. We use a diagonally preconditioned conjugate gradient solver.

      concepts::DiagonalMatrix<Real> diag(stiff);
      concepts::DiagonalSolver<Real> precond(diag);
      concepts::CG<Real> solver(stiff, precond, 1e-6, 2*stiff.dimX()+100);

Now, everything is prepared to solve the linear system.

sol is an empty vector to hold the solution of the linear system.

concepts::Vector<Real> sol(space);

Some immediate statistics of the solver.

The energy of the FE solution in sol can be computed by multiplying sol and rhs. It is compared to the exact energy.

      const Real exenergy = 1.8362264461350097;
      const Real FEenergy = sol*rhs;
      const Real relError = fabs(FEenergy - exenergy)/exenergy;
      std::cout << "  FE energy = " << FEenergy
                << ", rel. energy error = " << relError << std::endl;

If there are further computations (ie. i < 2), the space is refined using hp2D::APrioriRefinement towards the vertex carying the singularity (the origin). This information has to be given a-priori by the user.

      if (i < 2) {
        int pMax[2] = { 1, 1 };
        hp2D::APrioriRefinement refineSpace(space, 10, 0, pMax);
        concepts::GlobalPostprocess<Real> post(space);
        post(refineSpace);
      }
    }

Finally, exceptions are caught and a sane return value is given back.

  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }

Results

Output of the program:

msh = Import2dMesh(ncell = 3)
BoundaryConditions(1: Boundary(DIRICHLET, (0)), 2: Boundary(NEUMANN, ParsedFormula( (1/(1+y*y)^(2/3)) * (sin(2/3*atan(y)) - y*cos(2/3*atan(y)))*2/3 )), 3: Boundary(NEUMANN, ParsedFormula( 2/3/((1+x*x)^(2/3)) * (sin(2/3*atan(1/x)) + x*cos(2/3*atan(1/x))) )), 4: Boundary(NEUMANN, ParsedFormula( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(1/x)+pi)) + x*cos(2/3*(atan(1/x)+pi))) )), 5: Boundary(NEUMANN, ParsedFormula( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )), 6: Boundary(NEUMANN, ParsedFormula( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )), 7: Boundary(NEUMANN, ParsedFormula( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(-1/x)+pi)) - x*cos(2/3*(atan(-1/x)+pi))) )), 8: Boundary(DIRICHLET, (0)), 0: Boundary(FREE, (0)))
--
i = 0
  space = Space(dim = 5, nelm = 3)
  matrix: SparseMatrix(5x5, HashedSparseMatrix: 15 (60%) entries bound.)
  solver = CG(solves SparseMatrix(5x5, HashedSparseMatrix: 15 (60%) entries bound.), eps = 9.73773e-33, it = 2, relres = 0)
  FE energy = 1.7449, rel. energy error = 0.0497347
i = 1
  space = Space(dim = 16, nelm = 12)
  matrix: SparseMatrix(16x16, HashedSparseMatrix: 90 (35.1562%) entries bound.)
  solver = CG(solves SparseMatrix(16x16, HashedSparseMatrix: 90 (35.1562%) entries bound.), eps = 1.66559e-33, it = 6, relres = 0)
  FE energy = 1.79339, rel. energy error = 0.0233299
i = 2
  space = Space(dim = 50, nelm = 21)
  matrix: SparseMatrix(50x50, HashedSparseMatrix: 402 (16.08%) entries bound.)
  solver = CG(solves SparseMatrix(50x50, HashedSparseMatrix: 402 (16.08%) entries bound.), eps = 6.74911e-13, it = 17, relres = 0)
  FE energy = 1.81761, rel. energy error = 0.0101395 

Complete Source Code

Author: Philipp Frauenfelder, 2004