loading netgen mesh into ngsolve


I seem to have a problem importing a Netgen-generated mesh into Ngsolve for a PDE solution (Laplace eq. for capacitance calculation). I attached the files I used to generate the mesh: export.geo and export.py as well as the mesh file generated by those two files respectively: mesh.vol.gz and mesh_py.vol.gz

When I import mesh.vol.gz into the attached NGSpy script calculate.py (attached), none of the boundaries are imported with it and print(mesh.GetBoundaries()) produces (‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’). This causes a problem for the calculation since I use Dirichlet boundary conditions on selected boundaries.

When I import mesh_py.vol.gz into the same calculation script calculate.py, everything shows up correctly and I get the following list of boundaries (‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘default’, ‘2’, ‘2’, ‘2’, ‘2’, ‘2’, ‘2’, ‘1’, ‘1’, ‘1’, ‘1’, ‘1’, ‘1’)

The reason I can’t generate the mesh in NGSpy is that I need extruded polygons in the structure and the Python interface cannot extrude polygonal curves, only splines, and that doesn’t work for my structures.

If someone can please point me to a solution for the above problem I would very much appreciate the help. I suspect there is a simple mistake somewhere but cannot find it.

Thanks very much,



Attachment: calculate.py

Attachment: export.geo

Attachment: export.py

Attachment: mesh.vol.gz

Attachment: mesh_py.vol.gz

Hi Valery,

you are mixing up integers and strings.
In the export.py file the bcmod sets a “bcname”, which is a string, to “1”.
In the export.geo file the “boundarycondition” keyword set the boundary condition number to the integer 1. Thus the approaches give you different results.

Solution 1:
You could use the keyword “boundaryconditionname” and set the third argument to you bcname. Here you have to use some kind of string, just 1 or “1” is not recognized. These bcnames are properly stored in the mesh and you can use them in your calculate.py

Solution 2:
Stick to the keyword “boundarycondition” and set the bcnames in python after generating the Netgen mesh using “ngmesh.SetBCName(…)”. Be aware that the first argument is 0-based and Netgen boundary condtions are 1-based, thus 0 sets the first bc number.

geo = CSGeometry("export.geo")
ngmesh = geo.GenerateMesh(maxh=0.5)
mesh = Mesh(ngmesh)

Could you give a simple example of such an extrusion which is not working in with NGSpy and works with geo-files? Would you like to extrude a 2d geometry along a polygon or is the 2d geometry you extrude a polygon?


Hi Christoph,

Thanks very much for the prompt and helpful reply! I tried your solution #1 and it works very well like this:

boundaryconditionname p_M1_0_1 p_M1_0_1 e1;

In general I think I would prefer to keep everything in python, since I use python in calculate.py anyway and as far as I understood this is you preferred way of doing things, please correct me if I’m wrong. So if this all can be made to work in python, it does seem like a somewhat cleaner solution.

To answer your question about the structure I am trying to set up and simulate I attached a more complete set of files. The application is on-chip interconnect (capacitance) simulation for small cells such as memory (SRAM), image sensors, etc. The conductors (aluminum, copper or highly doped polysilicon) are defined as polygons on a number of layers. Each layer is planar with a defined thickness, there are a number of polygons on each layer. Some polygons are simple rectangles, some are more complex like p_M1_0_11 in the attached export.geo file. This works well in .geo but I have not been able to implement this type of extruded polygons in NGSpy.

In reality some of the polygons may need to be optionally represented by a much larger number of nodes. In this case the layout patterns (polygons) are obtained from a lithography simulation of the layout, I attached an example image (two layers of an SRAM cell). In this case a layout polygon would need typically 60-100 nodes for an adequate representation. So far this level of complexity seemed to stress Netgen a bit far, I attached an example to go with the above screenshot in export_large.geo. If you think it should be still doable in Netgen, perhaps with some simplification of the polygon, removing some nodes or similar, I would appreciate your advice on this.

I wonder if you may be able to help with another related question. The goal here is calculation of the capacitance matrix, i.e. the coupling capacitances between all electrodes in the structure. I attached an example matrix calculated by calculate.py. Capacitance is calculated as C = dQ/dV after applying a voltage (Dirichlet BC) to one electrode and 0V to all others. To calculate dQ I need a surface integral of -grad(u)*normal_vector on a closed surface around the electrode. However I cannot simple do:

     charge = (Integrate(gfE*bfv2), mesh, definedon=mesh.Boundaries("e"+str(i+1))))

This is because it turns out that bfv2 has some normal vectors pointing into and some out of the electrode surface:

func_domain2 = -CoefficientFunction (specialcf.normal(3))
bfv2 = BoundaryFromVolumeCF(func_domain2)

So as a workaround I use the absolute value of the scalar product instead:

     charge = (Integrate(sqrt((gfE*bfv2)**2), mesh, definedon=mesh.Boundaries("e"+str(i+1))))

This works because all normal field on the electrode surface must have the same sign (you can show this using the mean value theorem). Is there a better (cleaner) way to do this?

Thanks again for your help!




https://ngsolve.org/media/kunena/attachments/1538/Screenshot_20210428_113256.png [attachment=undefined]Screenshot_20210428_115219.png[/attachment] [attachment=undefined]Screenshot_20210428_115219.png[/attachment]

Hi Valery,

defining everything in python would be the way to go.
The attachments still do not show up. Could you try to upload them again or maybe upload all of them as zip archive?

I do not see a cleaner way for your integral. The normal vector depends on how you define the geometry and it might be painful (if possible) to get the orientation “correct” for complex geometries.
Just a comment, you should be able to you “specialcf.normal(3)” directly. I do not see why you should use the BoundaryFromVolumeCF.


Hi Christoph,

Sorry about the attachment problem, I used a zip archive this time. My current set of files is attached: export.geo defines the structure, calculate.py loads export.geo, meshes it and runs simulations. You can see that 4 of the geometry features on the Poly layer are trapezoidal, so in the .geo file I use an extrusion like this:

solid Poly_1 = ob and bot_Poly and top_Poly and extrusion(pathcurve2z;curve_Poly_1;0,1,0);

In the documentation at Constructive Solid Geometry CSG — NGS-Py 6.2.2302 documentation I only found half-space, sphere, cylinder and brick. After some searching I did find a Python extrusion function like this Extrusion(path,spline,Vec(0,1,0)), however only for path = SplineCurve3d() and spline = SplineCurve2d(). However, splines don’t work well for my polygonal structures. You can see my experiments with this in the attached file export_extr_ref.py, trying to create a small box using the extrusion function produced a curved object instead. My problem is with the spline created as SplineCurve2d(), if there is a way to use simple line segments like curve2d does in .geo, that would solve my problem.

Thanks also for the recommendation to simplify the calculation of surface normals, it worked well and is included in the attached calculate.py

Thanks and regards,



Attachment: ngsolve_forum_files2_2021-05-04-2.zip

Hi Valery,

attached is a python file generating the trapezoidal “Poly_1” extrusion object.

These spline objects can be composed of lines (2 points) and splines (3 points). The first argument of “Extrusion” (the extrusion path) has to be smooth, thus a single line is totally fine. The second argument can describe any 2d polygon you like.

Concerning you code:

 spline.AddSegment (0,1) 

adds a line segment from point 0 to point 1 to the spline and

 spline.AddSegment (0,1,2) 

adds a spline3 segment from point 0 to point 2 using point 1 as control point for the curvature.

The extrusion in “export_extr_ref.py” should be defined as follows.

spline = SplineCurve2d() 
spline.AddPoint (-1.1,-1.0)
spline.AddPoint (-1.1,-1.1)
spline.AddPoint (-1.2,-1.1)
spline.AddPoint (-1.2,-1.0)
spline.AddSegment (0,1)
spline.AddSegment (1,2)
spline.AddSegment (2,3)
spline.AddSegment (3,0)
extr = Plane (Pnt(0,0,0.2), Vec(0,0,-1) )*Plane (Pnt(0,0,0.3), Vec( 0,0,1) )*Extrusion(path,spline,Vec(0,1,0)).bc("1")


Attachment: extrude.py

Thanks very much! Exactly the answer I was looking for

Best regards,