Ctrl+K

A first example in NGSove

2.2. A first example in NGSove#

This first notebook introduces in few basic steps how to solve a simple PDE problem using NGSolve; each step contains different objects that will be explained in detail in the following notebooks.

To solve the problem one need to consider the following steps:

  1. Define the domain \(\Omega\), name the boundaries and mesh it.

  2. Define the finite element space and reserve space for the solution.

  3. Define the linear and bilinear forms and assemble respective matrix and a vector.

  4. Solve the linear system. (We may need homogeneization techniques)

2.2.1. Poisson problem on the unit square#

The first problem we are going to solve a problem on the unit square \(\Omega = (0,1)^2\), we call the boundaries \(\Gamma_b, \Gamma_t, \Gamma_l, \Gamma_r\) the bottom, top, left and right boundaries respectively.

Find \(u\) such that

\[\begin{align*} -\Delta u &= 20\,xy^2 \quad \text{in } \Omega, \\ u &= 0 \quad \text{on } \Gamma_b \cup \Gamma_l,\\ n\cdot \nabla u &= 0 \quad \text{on } \Gamma_r \cup \Gamma_t, \end{align*}\]

2.2.1.1. Step 0: Import the necessary libraries#

  • ngsolve is the library that contains the finite element methods and the solvers.

  • Draw is a function that allows us to visualize the results in the notebook.

from ngsolve import *
from ngsolve.webgui import Draw # Drawing in Jupyter notebook. More about this later.

2.2.1.2. Step 1: Define the domain#

netgen_mesh = unit_square.GenerateMesh(maxh=0.2) # generate a netgen mesh from the OCCGeometry class
mesh = Mesh(netgen_mesh)                         # create a NGSolve mesh class

Draw(mesh); # names are predefined for this case

2.2.1.3. Step 2: Define the finite element space and solution vector#

We need to create the Lagrange finite element space of order 1 on the mesh, to do so me need to pass a mesh and the essential boundary conditions.

The essential boundaries are selected using the dirichlet flag. In particular we select the bottom and (united with) the left boundary using the regex bottom|left.

fes = H1(mesh, order = 1, dirichlet="bottom|left") # define the finite element space

With the use of a finite element space one can create a GridFunction that will be used to store the solution.

For now you can understand the GridFunction as a vector that stores the values of the solution at the nodes of the mesh.

gfu = GridFunction(fes) # define the grid function: store a function on the finite element space

# help(gfu) # query the grid function

# print(gfu.vec) # print the vector of the grid function

2.2.1.4. Step 3: Define the symbolic forms#

Our problem is to find \(u \) such that

\[\begin{align*} \int_{\Omega} \nabla u \cdot \nabla v \,dx &= \int_{\Omega} 20\,xy^2 v \,dx \quad \forall v \end{align*}\]

For this specific problem we need to define the bilinear form \(a(u,v)\) and the linear form \(f(v)\). To write it in a math-like way we make use of the ProxyFunctions:

  • u, trial function.

  • v, test function.

u = fes.TrialFunction() # define the trial function
v = fes.TestFunction() # define the test function

# u, v = fes.TnT() # one-liner to define the trial and test functions

a = BilinearForm(fes) # define the bilinear form
a += grad(u) * grad(v) * dx
a += u*v * dx 

# print(a.mat)  # print it ... does not work! 

a.Assemble();

#print(a.mat) # print the matrix of the bilinear form

in the biliear form we have the dx that represents the integration over the domain \(\Omega\), it is short hand for DifferentialSymbol(VOL).

and grad that represents the gradient of the symbolic test/trial functions.

We want to declare the right hand side as follows: $\( f(v) = \int_{\Omega} 20 \,xy^2\, v \, dx. \)$ To do so we use the CoefficientFunctions x and y that represent the first two coordinate

f = LinearForm(fes)
f += 20 * x * y ** 2 * v * dx
f.Assemble();

# print(f.vec) # print the vector of the linear form. RHS is zero so far.

2.2.1.5. Step 4: Solve the linear system#

So far our problem was defined only symbolically, the only computation that was done was in the creation of the mesh. To actively solve the problem we need to assemble the linear system.

We can now invert the matrix and solve the linear system.

inv = a.mat.Inverse(freedofs=fes.FreeDofs()) # compute the inverse of the matrix

gfu.vec.data = inv * f.vec  # solve the system

What freedofs does is to select the degrees of freedom that are not in the essential boundary conditions.

Draw(gfu); # draw the solution

2.2.2. Interact with NGSolve#

  • A jupyter notebook (like this one) gives you one way to interact with NGSolve. When you have a complex sequence of tasks to perform, the notebook may not be adequate.

  • You can write an entire python module in a text editor and call python on the command line. (A script of the above is provided in poisson.py.)

    python3 poisson.py

  • If you want the Netgen GUI, then use netgen on the command line:

    netgen poisson.py

    You can then ask for a python shell from the GUI’s menu options (Solve -> Python shell).

  • One can use the Netgen GUI within a jupyter notebook: import netgen.gui.

    import netgen.gui
from ngsolve import *

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))  # generate a mesh
Draw(mesh);  # draw the mesh

# import netgen.gui
# from ngsolve import *

# mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))  # generate a mesh
# Draw(mesh);  # draw the mesh
Contents