Torsion warping function not working

I am trying to solve torsion Laplace equation with Neuman boundary
dw/dn = ynx-xny

Following code works well for a ellipse geometry, but it doesn’t work for a circle one. Not sure what’s the difference in between.

#Construct ngsolve FE space with direchlet boundary at out of section

fes = H1(mesh, order=3)#,dirichlet = "bcOuter")

#Declare test function, trial function, and grid function

u = fes.TrialFunction() # symbolic object

v = fes.TestFunction() # symbolic object

gfu = GridFunction(fes) # solution


a = BilinearForm(fes)

f = LinearForm(fes)


#Set Neumann boundary expression

n = specialcf.normal(mesh.dim)


a += grad(u)*grad(v)*dx

f += R*v*dx + G*v*ds(definedon="bcOuter")



aInverse = a.mat.Inverse(freedofs=fes.FreeDofs()) = aInverse * f.vec

The ellipse gives right solution:

However, for circular section, the result is not right

for the circle we expect the zero-function, right ?

Your matrix is singular, you are solving a pure Neumann problem. We depend heavily how the direct solver treats roundoff errors.

Simple way out: add a small regularization term to the BilineraForm (something like 1e-6uv*dx), or use a Lagrange multiplier for the mean value.

Try also curved elements (mesh.Curve(5)) to get the orthogonality of ynx-xny on the circle at high accuracy.