The Poisson equation:

\begin{eqnarray*} -\Delta\mathbf{u} &=& \mathbf{f} \quad \text{in} \quad \Omega, \\ \mathbf{u} &=& \mathbf{g} \quad \text{on} \quad \Gamma_D, \\ \mathbf{n} \cdot \nabla \mathbf{u} &=& \mathbf{h} \quad \text{on} \quad \Gamma_N, \end{eqnarray*}

where \(\mathbf{n}\) is the outward pointing unit normal vector, \(\Gamma_D\) and \(\Gamma_N\) are the parts of the boundary where Dirichlet (essential) and Neumann (natural) boundary conditions are prescribed, respectively.

We can also create a own .cpp and solve more complicated Poisson problems. In this tutorial, we explain how to do this. As an example we will use the poisson_example.cpp file in the example folder. We will go through essentials of the code in details (For more detail of the solver see the documentation for gsPoissonAssembler, gsAssembler, and gsVisitorPoisson.)

First we include the library and use the gismo namespace.

We add a command-line argument which allows the user to create a Paraview visualization file.

Our Poisson problem is defined by source (right-hand side) function and Dirichlet function data:

Geometry

Inside the main function we declare the geometry. Note that we use the gsMultiPatch class even if we only have a single patch geometry. The gsNurbsCreator class creates several simple geometries:

The chosen geometry represents the unit square consisting of 2 x 2 squares, which all have side-length 0.5, i.e., the computational domain is the unit square \(\Omega = (0,1)^2\) .

In this example, the four patches are numbered as follows:

1	3
0	2

Boundary Conditions

The boundary conditions are stored in an object of type gsBoundaryConditions :

In the 2D case, each patch has 4 boundary sides: west, east, south and north (in addition there is front and back for 3D geometries). The boundary conditions must be specified for each patch.

In this example we have 4 patches forming a unit square, hence we need to specify 8 boundary conditions:

Note that the first argument in the addCondition is the patch number.

Discretization Basis

We now need a discretization basis. We use the basis from the patches. We can h or/and p refine this basis (in any order).

Here, we perform h- and k-refinement.

Setting up and using the Assembler

We use gsPoissonAssembler class to assemble and the Poisson equation. First we initialize the class object by giving the geometry (patches), discretization basis, boundary condition, source function, and the strategies for treating Dirichlet boundary conditions and patch interfaces:

With dirichlet::elimination, we specify that the degrees of freedom (DOF) on the Dirichlet bondary should be eliminated already during the assembling process.

By iFace::glue, we specify that the DOFs at interfaces should be "glued" together in the sense that corresponding DOFs will be identified with each other and treated as one global DOF (which results in \(C^0\)-continuity at the interfaces).

Executing the assembling process is then done by calling

assembler.assemble();

The assembled system matrix and the right-hand-side can be accessed via

assembler.matrix();

assembler.rhs();

where the system matrix is stored as a gsSparseMatrix and the right-hand-side as gsMatrix.

Computing and Post processing the solution

Solving the linear system of equations is done in the main file. In this example of the Poisson problem, where we have a symmetric and positive definite system matrix, we can, for example, use the CG-solver of the Eigen library, with a diagonal preconditioner.

The solVector contains the computed coefficients without the Dirichlet DOFs (which were eliminated). In order to obtain the discrete solution as a gsField, we call

Finally we want to evaluate the (exact and approximate) solutions and write to a file:

the solutions can now opened and viewed in paraview.

Annotated source file

Here is the full file examples/poisson_example.cpp. Clicking on a function or class name will lead you to its reference documentation.