Torsion warping function not working

General
1

Hi,
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

#gfu.Set(555,definedon=mesh.BBoundaries("FixedPoint"))

a = BilinearForm(fes)

f = LinearForm(fes)

R=0

#Set Neumann boundary expression

n = specialcf.normal(mesh.dim)

G=InnerProduct(CF((y,-x)),n)

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

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

a.Assemble()

f.Assemble()

aInverse = a.mat.Inverse(freedofs=fes.FreeDofs())

gfu.vec.data = aInverse * f.vec

The ellipse gives right solution:
image

However, for circular section, the result is not right

image
image940×467 68.4 KB

2

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.

Joachim