Attaching geom2d to mesh

Netgen
1

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 :slight_smile:

disk.py (485 Bytes)

disk_csg.py (523 Bytes)

2 
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?

3

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")
4

Hi @christopher! Indeed @pbrubeck script reproduces what ngsPETSc does under the hood :slight_smile:

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.

5

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.

6

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