# 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) on the L shaped domain (-1,1)2 \ (0,1)x(-1,0) with mixed boundary conditions chosen such that the exact solution is . 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 .

The relative error in energy norm is computed.

# 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.

## 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.

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().

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

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

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.

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

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

Now, everything is prepared to solve the linear system.

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

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.

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.

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

# 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

#include <iostream>
#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"
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;
// ** boundary conditions **
("( (1/(1+y*y)^(2/3)) * (sin(2/3*atan(y)) - y*cos(2/3*atan(y)))*2/3 )")));
("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*atan(1/x)) + x*cos(2/3*atan(1/x))) )")));
("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(1/x)+pi)) + x*cos(2/3*(atan(1/x)+pi))) )")));
("( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )")));
("( 2/3/((1+y*y)^(2/3)) * (sin(2/3*(atan(-y)+pi)) + y*cos(2/3*(atan(-y)+pi))) )")));
("( 2/3/((1+x*x)^(2/3)) * (sin(2/3*(atan(-1/x)+pi)) - x*cos(2/3*(atan(-1/x)+pi))) )")));
std::cout << bc << std::endl;
std::cout << "--" << std::endl;
// ** space **
hp2D::Space space(msh, 0, 1, &bc);
for (uint i = 0; i < 3; ++i) {
std::cout << "i = " << i << std::endl;
space.rebuild();
std::cout << " space = " << space << std::endl;
// ** right hand side **
concepts::Vector<Real> rhs(space, lform);
// ** matrices **
concepts::SparseMatrix<Real> stiff(space, la);
stiff.compress();
std::cout << " matrix: " << stiff << std::endl;
// ** set up solver **
concepts::CG<Real> solver(stiff, precond, 1e-6, 2*stiff.dimX()+100);
// ** solve **
solver(rhs, sol);
std::cout << " solver = " << solver << std::endl;
// ** error **
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;
// ** refinement **
if (i < 2) {
int pMax = { 1, 1 };
hp2D::APrioriRefinement refineSpace(space, 10, 0, pMax);
post(refineSpace);
}
}
}
std::cout << e << std::endl;
return 1;
}
return 0;
}
Solves a symmetric system of linear equations with conjugate gradients (CG).
Definition: cg.hh:39
Base class for exceptions.
Definition: exceptions.hh:86
void add(const Set< Attribute > &attrib, const Boundary &bcObject)
Adds a boundary condition with this attribute to the list of boundary conditions.
A solver for diagonal matrices.
Definition: diagonal.hh:142
void compress(Real threshold=EPS)
Compresses the matrix by dropping small entries.
Space rebuild()
Rebuilds the space.
Class to describe an element of the boundary.
Definition: boundary.hh:35
Linear form in 2D.
Definition: linearForm.hh:62
Diagonal matrix.
Definition: diagonal.hh:24
double pi()
A function class to calculate element matrices for the Laplacian.
Definition: bf_laplace.hh:91
A 2D hp FEM space with continuous, piecewise polynomial basis functions.
Definition: space.hh:88
Carries out a-priori given refinements.
Definition: aprioriRef2D.hh:40
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich