Inhomogeneous Dirichlet Boundary Conditions

In case of inhomogeneous Dirichlet boundary conditions, we combine the technique of Dirichlet data extension with the above static condensation principle in the following code.

U = x*x*(1-y)*(1-y)          # U = manufactured solution
DeltaU = 2*((1-y)*(1-y)+x*x) # Source: DeltaU = ∆U

f = LinearForm(fes)
f += -DeltaU * v * dx
f.Assemble()

u = GridFunction(fes)
u.Set(U, BND)               # Dirichlet b.c: u = U on boundary

r = f.vec.CreateVector()
r.data = f.vec - a.mat * u.vec
r.data += a.harmonic_extension_trans * r

u.vec.data += a.mat.Inverse(fes.FreeDofs(True)) * r
u.vec.data += a.harmonic_extension * u.vec
u.vec.data += a.inner_solve * r

sqrt(Integrate((U-u)*(U-u),mesh))  # Compute error

The code you provided is solving a Poisson equation with a manufactured solution and then calculates the error between the manufactured solution and the numerical solution obtained.

The first two lines define the manufactured solution U and the source term DeltaU for the Poisson equation.

The next few lines set up the right-hand side of the Poisson equation, by creating a linear form f and adding the term -DeltaU * v * dx to it. The Assemble() function is called to assemble the linear form.

Then, a grid function u is created and the manufactured solution U is set as the Dirichlet boundary condition using u.Set(U, BND)

Next, the residual is calculated as r = f.vec - a.mat * u.vec. The residual is then corrected to account for the effect of the eliminated degrees of freedom using the a.harmonic_extension_trans * r

The solution is then corrected by adding a.mat.Inverse(fes.FreeDofs(True)) * r to it. And the solution is then corrected again by adding a.harmonic_extension * u.vec and a.inner_solve * r

Finally, the last line of code calculates the error between the manufactured solution U and the numerical solution u by computing the L2-norm of the difference using Integrate((U-u)*(U-u),mesh) and taking the square root.

This script can be useful to evaluate the accuracy of the numerical solution and the effect of the static condensation technique on the solution.

If you find the above steps onerous, you can again use the automatic utility solvers.BVP, which will perform static condensation if condense is turned on in its input bilinear form.

Automatic utility BVP

u = GridFunction(fes)
u.Set(U, BND)

c = Preconditioner(a,"direct")
c.Update()
solvers.BVP(bf=a, lf=f, gf=u, pre=c)
sqrt(Integrate((U-u)*(U-u),mesh))

The code you provided is also solving a Poisson equation with a manufactured solution and then calculates the error between the manufactured solution and the numerical solution obtained, but it is using the solvers.BVP function of the ngsolve library to solve the system of equations.

The first two lines are the same as before, the manufactured solution U is set as the Dirichlet boundary condition using u.Set(U, BND)

Then, the solvers.BVP function is used to solve the Poisson equation. The function takes three arguments: the bilinear form a, the linear form f, and the grid function u to store the solution. It also takes a fourth argument, the preconditioner c. The preconditioner is used to accelerate the convergence of the iterative solver. The preconditioner is set to “direct” which uses a direct solver to invert the matrix.

The last line calculates the error between the manufactured solution U and the numerical solution u by computing the L2-norm of the difference using Integrate((U-u)*(U-u),mesh) and taking the square root.

The solvers.BVP function is a high-level function that automatically handles the solution of a boundary value problem, it takes care of assembling the matrix, applying the boundary conditions, and solving the linear system.

This script can be useful to evaluate the accuracy of the numerical solution and the effect of the preconditioner on the solution.