CoordinateTrafo causes segmentation fault with GridFunction-based CoefficientFunction (regression in 6.2.2503)

Hi,

I found a regression bug in NGSolve. Using CoordinateTrafo with a CoefficientFunction derived from GridFunction causes segmentation faults in versions 6.2.2503 through 6.2.2506.

Environment

  • OS: Windows 10/11

  • Python: 3.12

  • NGSolve versions tested: 6.2.2406 ~ 6.2.2506

Version Test Results

| Version | Result | Notes |

|---------|--------|-------|

| 6.2.2406 | OK | Works correctly |

| 6.2.2501 | OK | Works correctly |

| 6.2.2502 | N/A | Could not install (dependency issue with netgen-mesher) |

| 6.2.2503 | CRASH | Segmentation fault |

| 6.2.2504 | CRASH | Segmentation fault |

| 6.2.2505 | CRASH | Segmentation fault |

| 6.2.2506 | CRASH | Segmentation fault |

| 6.2.2506-139-ga14a97920 | OK | Works correctly (nightly build) |

Regression introduced in 6.2.2503, fixed in nightly build (6.2.2506-139-ga14a97920).

Minimal Reproducible Example


import ngsolve

print(ngsolve.__version__)

from numpy import *

from ngsolve import *

from ngsolve.webgui import Draw

from netgen.occ import *

from ngsolve.fem import CoordinateTrafo

# Create mesh

cyl = Cylinder(Pnt(0, 0, 0), Z, r=2, h=1)

geo = OCCGeometry(cyl)

mesh = Mesh(geo.GenerateMesh(maxh=0.3))

# Create HCurl GridFunction

fesHCurl = HCurl(mesh, order=1)

theta = atan2(y, x)

gfE = GridFunction(fesHCurl)

gfE.Set((-y * IfPos(theta, 1, 0) * IfPos(pi/4 - theta, 1, 0),

x * IfPos(theta, 1, 0) * IfPos(pi/4 - theta, 1, 0), 0))

# Create CoefficientFunction from GridFunction

Q = 0.5 * InnerProduct(gfE, gfE)

# Apply CoordinateTrafo - THIS CAUSES CRASH

phi = pi/4

trafo_cf = CF((cos(phi)*x + sin(phi)*y, -sin(phi)*x + cos(phi)*y, z))

trafo = CoordinateTrafo(trafo_cf, mesh.Materials(".*")) # <-- Crash here or next line

Q_rotated = Q(trafo)

# Set to GridFunction

fesH1 = H1(mesh, order=1)

gfQ = GridFunction(fesH1)

gfQ.Set(Q_rotated) # <-- Or crash here

Draw(gfQ, mesh, "Q", order=1)

Expected Behavior

The code should apply coordinate transformation to the CoefficientFunction and display the rotated field without crashing.

Actual Behavior

  • In versions 6.2.2503 ~ 6.2.2506: Segmentation fault occurs

  • Crash happens during script execution (not always at the same point)

  • Versions 6.2.2406 and 6.2.2501 work correctly

Working Versions

  • 6.2.2406: OK

  • 6.2.2501: OK

  • 6.2.2506-139-ga14a97920 (nightly/development build): OK

Conclusion

This is a regression bug introduced in version 6.2.2503. The bug has been fixed in the nightly build (6.2.2506-139-ga14a97920).

Workaround: Use version 6.2.2501 or earlier until the fix is released in a stable version.

Questions

  1. Is the fix in nightly build planned to be included in the next stable release (6.2.2507)?

  2. Is there a specific commit that fixed this issue?

Thank you for your help.


Appendix: Full Test Script


import sys

from numpy import pi, cos, sin

import ngsolve

version = ngsolve.__version__

print(f"NGSolve version: {version}")

from ngsolve import Mesh, HCurl, H1, GridFunction, CF, InnerProduct, atan2, IfPos, x, y, z

from netgen.occ import Cylinder, Pnt, Z, OCCGeometry

from ngsolve.fem import CoordinateTrafo

# Step 1: Create mesh

cyl = Cylinder(Pnt(0, 0, 0), Z, r=2, h=1)

geo = OCCGeometry(cyl)

mesh = Mesh(geo.GenerateMesh(maxh=0.3))

# Step 2: Create GridFunction

fesHCurl = HCurl(mesh, order=1)

theta = atan2(y, x)

gfE = GridFunction(fesHCurl)

gfE.Set((-y * IfPos(theta, 1, 0) * IfPos(pi/4 - theta, 1, 0),

x * IfPos(theta, 1, 0) * IfPos(pi/4 - theta, 1, 0), 0))

Q = 0.5 * InnerProduct(gfE, gfE)

# Step 3: Create CoordinateTrafo

phi = pi/4

trafo_cf = CF((cos(phi)*x + sin(phi)*y, -sin(phi)*x + cos(phi)*y, z))

trafo = CoordinateTrafo(trafo_cf, mesh.Materials(".*"))

# Step 4: Apply transformation

Q_rotated = Q(trafo)

# Step 5: Set to GridFunction

fesH1 = H1(mesh, order=1)

gfQ = GridFunction(fesH1)

gfQ.Set(Q_rotated)

# Step 6: Multiple rotations

for n in range(4):

phi = n * pi / 4

trafo_cf = CF((cos(phi)*x + sin(phi)*y, -sin(phi)*x + cos(phi)*y, z))

trafo = CoordinateTrafo(trafo_cf, mesh.Materials(".*"))

Q_rotated = Q(trafo)

gfQ.Set(Q_rotated)

print(f"[{version}] SUCCESS")