Cylinder Coordinates / Integration rule

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

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

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