This page was generated from unit-4.3-manualmesh/manualmeshing.ipynb.

4.3 Working with meshes

One-dimensional meshes

Meshes in one-dimension can be constructed using the netgen.meshing module. We just have to add segments (Element1D) and the boundaries (Element0D).

import netgen.gui
%gui tk
from netgen.meshing import *

First we create a new Mesh and set the spatial dimension to 1.

m = Mesh()
m.dim = 1

Then we define and add the MeshPoint’s to the mesh. The function m.Add returns PointId’s which we store in an array to be able to construct the segments in the next step.

N = 10
pnums = []
for i in range(0, N+1):
    pnums.append (m.Add (MeshPoint (Pnt(i/N, 0, 0))))

type(pnums[0])

netgen.libngpy._meshing.PointId

Now we can loop over the array with the PointId’s and add one-dimensional elements to the mesh. Further we can set the material for our domain.

for i in range(0,N):
    m.Add (Element1D ([pnums[i],pnums[i+1]], index=1))

m.SetMaterial(1,'material')

Finally we have to add the boundary elements and set boundary conditions.

m.Add (Element0D (pnums[0], index=1))
m.Add (Element0D (pnums[N], index=2))

To be able to visualize one-dimensional meshes and solution activate Show edges in the menu View > Viewing options > Mesh.

import ngsolve
mesh = ngsolve.Mesh(m)

Two-dimensional meshes

As example we mesh a unit square [0,1]x[0,1] using quadrilaterals.

from netgen.geom2d import unit_square

We create an empty mesh and initialize the geometry and the dimension

ngmesh = Mesh()
ngmesh.SetGeometry(unit_square)
ngmesh.dim = 2

and add all the MeshPoint’s we will need for the final mesh. Similar to the one-dimensional mesh we store the PointId’s in the pnums array.

N=5
pnums = []
for i in range(N + 1):
    for j in range(N + 1):
        pnums.append(ngmesh.Add(MeshPoint(Pnt(i / N, j / N, 0))))

Next, we add the quadrilaterals to the mesh. Before that we have so add a FaceDescriptor to the mesh.

ngmesh.Add (FaceDescriptor(surfnr=1,domin=1,bc=1))
ngmesh.SetMaterial(1, "mat")
for j in range(N):
    for i in range(N):
        ngmesh.Add(Element2D(1, [pnums[i + j * (N + 1)],
                               pnums[i + (j + 1) * (N + 1)],
                               pnums[i + 1 + (j + 1) * (N + 1)],
                               pnums[i + 1 + j * (N + 1)]]))

Finally we have to add boundary elements and set boundary conditions.

# horizontal boundaries
for i in range(N):
   ngmesh.Add(Element1D([pnums[N + i * (N + 1)],
                       pnums[N + (i + 1) * (N + 1)]], index=1))
   ngmesh.Add(Element1D([pnums[0 + i * (N + 1)], pnums[0 + (i + 1) * (N + 1)]], index=1))

# vertical boundaries
for i in range(N):
   ngmesh.Add(Element1D([pnums[i], pnums[i + 1]], index=2))
   ngmesh.Add(Element1D([pnums[i + N * (N + 1)], pnums[i + 1 + N * (N + 1)]], index=2))

from ngsolve.solve import Draw
mesh = ngsolve.Mesh(ngmesh)
Draw(mesh)

Merge three-dimensional meshes

In the following example we will merge two surface meshes and generate a unified volume mesh.

from netgen.meshing import *
from netgen.csg import *

from ngsolve import ngsglobals
ngsglobals.msg_level = 2

As starting point we create two geometries and mesh them.

# generate brick and mesh it
geo1 = CSGeometry()
geo1.Add (OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) ))
m1 = geo1.GenerateMesh (maxh=0.1)
# m1.Refine()

# generate sphere and mesh it
geo2 = CSGeometry()
geo2.Add (Sphere (Pnt(0.5,0.5,0.5), 0.1))
m2 = geo2.GenerateMesh (maxh=0.05)
m2.Refine()
# m2.Refine()

Now we start the merging process. Therefore we create an empty mesh and add a FaceDescriptor for each of the surfaces.

# create an empty mesh
ngmesh = Mesh()

# a face-descriptor stores properties associated with a set of surface elements
# bc .. boundary condition marker,
# domin/domout .. domain-number in front/back of surface elements (0 = void),
# surfnr .. number of the surface described by the face-descriptor

fd_outside = ngmesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1))
fd_inside = ngmesh.Add (FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))

Since the surface elements stay the same in the merged mesh, we copy the points on the surface and the surface elements to the new mesh.

# copy all boundary points from first mesh to new mesh.
# pmap1 maps point-numbers from old to new mesh
pmap1 = { }
for e in m1.Elements2D():
    for v in e.vertices:
        if (v not in pmap1):
            pmap1[v] = ngmesh.Add (m1[v])


# copy surface elements from first mesh to new mesh
# we have to map point-numbers:

for e in m1.Elements2D():
    ngmesh.Add (Element2D (fd_outside, [pmap1[v] for v in e.vertices]))

# same for the second mesh:
pmap2 = { }
for e in m2.Elements2D():
    for v in e.vertices:
        if (v not in pmap2):
            pmap2[v] = ngmesh.Add (m2[v])

for e in m2.Elements2D():
    ngmesh.Add (Element2D (fd_inside, [pmap2[v] for v in e.vertices]))

Finally we have to generate the new volume mesh.

ngmesh.GenerateVolumeMesh()
import ngsolve
mesh = ngsolve.Mesh(ngmesh)
Draw(mesh)