Mesh refinement with Netgen

General
1

Hello everyone!

I’m moving my code for adaptive refinement algorithm from Firedrake with Netgen to FEniCSx.

The initial meshes for both platforms are the same:

initial
initial1840×1152 78.7 KB

which is generated by the following

from mpi4py import MPI
from netgen.csg import OrthoBrick, Pnt, CSGeometry
import ngsPETSc.utils.fenicsx as ngfx
from dolfinx.fem import (
    functionspace, Function
)
import numpy as np

cube1 = OrthoBrick(Pnt(-1,-1,0), Pnt(0,0,1))
cube2 = OrthoBrick(Pnt(0,-1,0), Pnt(1,0,1))
cube3 = OrthoBrick(Pnt(-1,0,0), Pnt(0,1,1))
geo = CSGeometry()
geo.Add(cube1 + cube2 + cube3)
geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD)
domain, (ct, ft), region_map = geoModel.model_to_mesh(gdim=3, hmax=1.0)

The uniform refinement in Firedrake is implemented with the following code

# Firedrake code
msh = msh.refine_marked_elements(
        Function(FunctionSpace(msh, 'DG', 0)).interpolate(1.0))

which generates the following finer mesh

firedrake
firedrake1840×1152 142 KB

I tried the Netgen refine in FEniCSx but it produces an error:

DG0 = functionspace(domain, ('DG', 0))
markers = Function(DG0)
markers.x.array[:] = 1.0
markers = np.flatnonzero(np.isclose(markers.x.array.astype(np.int32), 1))
domain, (_, ft) = geoModel.refineMarkedElements(domain.topology.dim, markers)

with the error message:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[8], line 24
     21 markers.x.array[:] = 1.0
     23 markers = np.flatnonzero(np.isclose(markers.x.array.astype(np.int32), 1))
---> 24 domain, (_, ft) = geoModel.refineMarkedElements(domain.topology.dim, markers)

File /dolfinx-env/lib/python3.12/site-packages/ngsPETSc/utils/fenicsx.py:619, in GeometricModel.refineMarkedElements(self, dim, elements, netgen_flags)
    617         el.refine = False
    618 if not refine_faces and dim == 3:
--> 619     _dim_to_element_wrapper(self.ngmesh)[2]().Numpy()["refine"] = 0
    620 self.ngmesh.Refine(adaptive=True)
    621 self.ngmesh.Curve(1)  # Reset mesh to be linear

AttributeError: 'netgen.libngpy._meshing.Array_N6netgen9Element2dE_N6netgen19SurfaceElementIndexE' object has no attribute 'Numpy'

I also try the uniform refinement imported as from dolfinx.mesh import uniform_refine but it gives a different result:

fenicsx
fenicsx1840×1152 147 KB

I’m trying to get the same refined mesh in FEniCSx that I got in Firedrake. I’m using the FEniCSx environment described in the tutorial: Adaptive mesh refinement with NetGen and DOLFINx — FEniCSx tutorial.

I’d appreciate any guidance on this issue. Thank you!

2

I realized that’s not an intelligent question. After double-checking the ngsPETSc documentation (see ngsPETSc/ngsPETSc/utils/fenicsx.py at main · NGSolve/ngsPETSc), it becomes clear that I need to set the “refine_faces” = True, like this:

domain, (_, ft) = geoModel.refineMarkedElements(
        domain.topology.dim,
        markers,
        netgen_flags={'refine_faces': True})

This now produces the result I want.

3

In addition, I think there might be a typo in line 619 ngsPETSc/ngsPETSc/utils/fenicsx.py at main · NGSolve/ngsPETSc, where Numpy() should be NumPy().

For the above uniform refinement in 3D, that typo doesn’t affect the results. But if you want to apply the adaptive refinement and don’t want to refine faces, currently I think the following works:

geoModel.ngmesh.Elements2D().NumPy()['refine'] = 0

Thank you!

4

Yes, that seems like a typo. I’ve made a PR to fix it:

5

Thanks for fixing and replying so quickly!