Cylinder Coordinates / Integration rule

Main Forum Old-Forum - do not post here
1

Hi everyone,

currently I am trying to compute a vector valued Poisson problem in cylindrical coordinates, under the assumption of rotational symmetry which allows the problem to be reduced into a 2-dimensional one. The weak formulation of this is the following:
[tex]
\int_{\Omega^{2d}} x \nabla u : \nabla v + \frac{1}{x} u_x v_x dx dy
[/tex]

  • The term 1/x * u_x * v_x term is of course singular on the rotational axis (r → x= 0). However I have u_x = 0 on the boundary x = 0 so I do not have to solve on this boundary.
  • In order to not evaluate anything on the boundary, I have manually changed the integration rule to only use interior points.

However, the stiffness-matrix still contains ‘nan’ values for the dofs (see attached MWE).
https://ngsolve.org/media/kunena/attachments/830/Poission_CylCoord_RotSym.py

Can anyone explain/solve this problem?

Best wishes,
Henry

Attachment: Poission_CylCoord_RotSym.py

2

Hi Henry,

it is fixed in the coming nightly.
The integration with simd was using IntegrationPoints = (0,0,0), with 0 weights, for the overhead, which caused the division 0/0.

it works with your current version if you disable simd by SymbolicBFI(…, simd_evaluate=False)

Joachim

3

Hi Joachim,

thank you for the quick reply and fix!

Is there an internal way to get integration rules which do not use points on the element boundary, or do I always have to specify them manually?

Best wishes,
Henry