We are developing support for multigrid starting from a Netgen mesh in Firedrake. Doing this I realised that attaching a geom2d to a mesh seg faults, while it works perfectly fine with occ. See attached, scripts 
disk.py (485 Bytes)
disk_csg.py (523 Bytes)
from netgen.geom2d import *
geo = SplineGeometry()
geo.AddCircle(c=(0, 0),
r=1.0,
bc="circle")
ngmesh = geo.GenerateMesh(maxh=0.5)
ngmesh.SetGeometry(None)
ngmesh.Save("disk_original.vol")
from netgen.meshing import Mesh
from ngsolve import *
mesh = Mesh("disk_original.vol")
mesh.ngmesh.SetGeometry(geo)
mesh.Refine()
Draw(mesh)
Without ngspetsc it works for me on master. Does this work for you?
The issue is not to do with refining, but creating a mesh from a list of vertices. This should reproduce the issue without the ngsPETSc dependency:
import netgen.meshing as ngm
# occ works with shift = 0 and -1
# csg works with shift = 0 but fails with -1
#backend = "occ"
backend = "csg"
# shift = 0
shift = -1
if backend == "occ":
from netgen.occ import Circle, OCCGeometry
disk = Circle((0,0), 1).Face()
geo = OCCGeometry(disk, dim=2)
elif backend == "csg":
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddCircle(c=(0, 0), r=1.0)
else:
raise ValueError("Unexpected backend")
ngmesh = geo.GenerateMesh(maxh=0.5)
ngmesh.Save("disk_original.vol")
dim = ngmesh.dim
geo = ngmesh.GetGeometry()
coordinates = ngmesh.Coordinates()
rngmesh = ngm.Mesh(dim=dim)
rngmesh.SetGeometry(geo)
rngmesh.AddPoints(coordinates)
cells_np = ngmesh.Elements2D().NumPy()
nodes = cells_np["nodes"][:, :dim+1]
index = cells_np["index"][0]
rngmesh.Add(ngm.FaceDescriptor(bc=index))
rngmesh.AddElements(dim=dim,
index=index,
data=nodes, base=1,
project_geometry=True)
faces = ngmesh.Elements1D()
for f in faces:
edge = ngm.Element1D(f.vertices, index=f.index, edgenr=f.index+shift)
rngmesh.Add(edge, project_geominfo=True)
rngmesh.Refine()
rngmesh.Curve(2)
rngmesh.Save("disk.vol")
Hi @christopher! Indeed @pbrubeck script reproduces what ngsPETSc does under the hood
With one caveat, we use base=0 (ngsPETSc/ngsPETSc/plex.py at 582948328382130037daae80167fd4abb6cd3ad1 · NGSolve/ngsPETSc · GitHub) and pass the edgenr ( ngsPETSc/ngsPETSc/plex.py at 582948328382130037daae80167fd4abb6cd3ad1 · NGSolve/ngsPETSc · GitHub ) to the boundary edges.
Hi @christopher, the ngsPETSc code needs edgenr and index to differ by a -1 shift. Could you help us understand why is this needed?
Also I tried to define 1D boundary elements via AddElements(dim=1, ...) but this will not project them to the geometry, whereas this works in higher dimensions.
by historic reasons, there are different indices, which are either 0 or 1 based, for different geometry formats. Is not a problem for Netgen internal use, but it’s now an issue when interfacing other packages.
There is certainly this quick fix by adding +/-1 for certain geometry formats, but there is no promise that this will not change in future.
There is one robust recommendation:
```
fd_idx = rngmesh.Add(ngm.FaceDescriptor(bc=index))
rngmesh.AddElements(dim=dim,
index=fd_idx,
data=nodes, base=1,
project_geometry=True)
```
A FaceDescrtiptor stores data for a face (e.g. string label, a geometry index, …), and each individual element indexes the facedescriptor. The user shall not have to care if it’s 0-based or 1-based indexing.
In future, I like to have a similar EdgeDescriptor for curves, and each 1D element is indexing into that.
Recommendation is to add many elements using AddElements(dim=1, np.array, ..), We will have a look into that.
A strong recommendation is to focus on OCC-geometry.
SplineGeometry (2d), CSGeoemtry (3d) and csg2d are currently maintained, but may be dropped in future.
Joachim
Pablo’s example is now working, without any index-shifts, for occ and spline2d geometries.
Now in gitlab, if tests run through I’ll make a prerelease later today.
Now trucks are driving on the right, in the UK
Joachim
Thanks Joachim. I wanted to try your fix, but I am having issues to install netgen from the source code (skbuild ModuleNotFoundError, possibly because I am installing through a python venv).
Please let us know when the pre-release becomes available through pip.
It’d be really nice if you could enable EdgeDescriptor for us. It would greatly simplify the logic in ngsPETSc.
Hi Pablo,
the pre-release from yesterday (6.2.2602.post11) should work for you.
I restructured some indexing, which needs now testing. Introducing the EdgeDesriptor will be the next step, maybe in a month or so.
For the skbuild error it needs Matthias.
Joachim
I did pip install netgen-mesher==6.2.2602.post11.dev0. As you pointed out, the segfault can by avoided by correctly setting edgenr depending on the geometry backend. Nevertheless, the edge points are not projected to the boundary. In our case, we are refining via DMPlex, which is geometry-agnostic. I modified the example to reproduce the failure:
import netgen.meshing as ngm
#backend = "occ"
backend = "spline"
#backend = "csg2d"
if backend == "occ":
from netgen.occ import Circle, OCCGeometry
disk = Circle((0,0), 1).Face()
geo = OCCGeometry(disk, dim=2)
shift = -1
elif backend == "spline":
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddCircle(c=(0, 0), r=1.0)
shift = 0
elif backend == "csg2d":
from netgen.geom2d import CSG2d, Circle
geo = CSG2d()
geo.Add(Circle(center=(0, 0), radius=1, bc="surface"))
shift = 0
else:
raise ValueError("Unexpected backend")
ngmesh = geo.GenerateMesh(maxh=0.5)
dim = ngmesh.dim
geo = ngmesh.GetGeometry()
# This refines the mesh without projecting the points to the boundary
ngmesh.SetGeometry(None)
ngmesh.Refine()
coordinates = ngmesh.Coordinates()
# Construct a new mesh and project the points on the boundary
rngmesh = ngm.Mesh(dim=dim)
rngmesh.SetGeometry(geo)
rngmesh.AddPoints(coordinates)
cells_np = ngmesh.Elements2D().NumPy()
nodes = cells_np["nodes"][:, :dim+1]
index = cells_np["index"][0]
rngmesh.Add(ngm.FaceDescriptor(bc=index))
rngmesh.AddElements(dim=dim,
index=index,
data=nodes, base=1,
project_geometry=True)
faces = ngmesh.Elements1D()
for f in faces:
edge = ngm.Element1D(f.vertices, index=f.index, edgenr=f.index+shift)
rngmesh.Add(edge, project_geominfo=True)
rngmesh.Save("disk.vol")
you need
edge = ngm.Element1D(f.vertices, index=f.index, edgenr=f.edgenr)
then your testcase is working for me for all 3 backends
If I replace that line I get this error:
Traceback (most recent call last):
File "/home/brubeckmarti/git/dphil_thesis/static_condensation/snippets/bug_netgen.py", line 53, in <module>
edge = ngm.Element1D(f.vertices, index=f.index, edgenr=f.edgenr)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
1. netgen.libngpy._meshing.Element1D(vertices: list, surfaces: list = [], index: typing.SupportsInt = 1, edgenr: typing.SupportsInt = 1, trignums: list = [])
Invoked with: [22, 2]; kwargs: index=1, edgenr=18446744073709551615
Could this mean that your fix is not there?