This page was generated from unit-1.8-meshtopology/meshtopology.ipynb.

1.8 Exploring the mesh topology

In this tutorial, we learn how to iterate over mesh entities and obtain information about how one mesh entity is connected to others. We will be able to answer questions like the following after this tutorial.

• What edges are connected to a vertex?

• What are the edges of an element?

• What are the points of an element?

• What elements are adjacent to a face?

In some algorithms, such as custom block smoothers (that we shall see in 2.1.2), we will need to iterate over mesh entities using the concepts in this tutorial.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg import unit_cube

mesh = Mesh(unit_cube.GenerateMesh(maxh=1))
Draw(mesh)

 Start Findpoints
 main-solids: 1
 Found points 8
 Analyze spec points
 Find edges
 12 edges found
 Check intersecting edges
 CalcLocalH: 8 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 12 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 12 edges found
 Check intersecting edges
 Surface 1 / 6
 load internal triangle rules
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 Surface 2 / 6
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 Surface 3 / 6
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 Surface 4 / 6
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 Surface 5 / 6
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 Surface 6 / 6
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 2 elements, 8 points
 SplitSeparateFaces
 Check subdomain 1 / 1

 Meshing subdomain 1 of 1
 Use internal rules
 Delaunay meshing
 number of points: 8
 blockfill local h
 number of points: 8
 Points: 8
 Elements: 33
 Tree data entries per element: 3.0303
 Tree nodes per element: 0.030303
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 8 points, 0 elements
 start tetmeshing
 Use internal rules
 Success !
 10 points, 16 elements
 Remove Illegal Elements
 Volume Optimization
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 Update mesh topology
 Update clusters

BaseWebGuiScene

Iterating over mesh objects

Iterating over vertices:

for v in mesh.vertices:
    print (v, v.point)

V0 (0.0, 0.0, 0.0)
V1 (0.0, 0.0, 1.0)
V2 (1.0, 0.0, 0.0)
V3 (0.0, 1.0, 0.0)
V4 (1.0, 0.0, 1.0)
V5 (0.0, 1.0, 1.0)
V6 (1.0, 1.0, 0.0)
V7 (1.0, 1.0, 1.0)
V8 (0.5000000008190509, 0.5000000010144133, 0.5000000008188038)

What is the type of v in mesh.vertices?

for v in mesh.vertices:
    print(type(v))

<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>
<class 'ngsolve.comp.MeshNode'>

We’ll shortly return to this MeshNode class.

Iterating over elements:

for el in mesh.Elements(VOL):
    print(type(el))
    print ("vertices: ", el.vertices)   # get vertices of an element
    print ("edges: ", el.edges)         # get edges of an element

<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V3, V5, V7, V8)
edges:  (E16, E22, E25, E13, E15, E21)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V0, V1, V8, V2)
edges:  (E1, E6, E12, E0, E5, E9)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V1, V2, V4, V8)
edges:  (E9, E12, E20, E6, E7, E10)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V0, V3, V8, V5)
edges:  (E3, E13, E22, E2, E5, E16)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V1, V4, V5, V8)
edges:  (E9, E20, E22, E7, E8, E17)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V0, V2, V8, V6)
edges:  (E4, E11, E24, E1, E5, E12)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V0, V3, V6, V8)
edges:  (E5, E16, E24, E2, E4, E14)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V0, V1, V5, V8)
edges:  (E5, E9, E22, E0, E3, E8)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V2, V4, V8, V6)
edges:  (E11, E18, E24, E10, E12, E20)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V3, V6, V8, V7)
edges:  (E15, E23, E25, E14, E16, E24)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V4, V6, V7, V8)
edges:  (E20, E24, E25, E18, E19, E23)
<class 'ngsolve.comp.Ngs_Element'>
vertices:  (V4, V5, V8, V7)
edges:  (E19, E21, E25, E17, E20, E22)

NodeId and MeshNodes

An object of type NodeId in NGSolve is just a number together with a type of the mesh entity it describes.

v = NodeId(VERTEX,0)   # a standalone NodeId object
type(v)
print ("type = ", v.type, "v.nr =", v.nr)

type =  NODE_TYPE.VERTEX v.nr = 0

NODE_TYPE can be one of the following:

• VERTEX: dimension 0

• EDGE: dimension 1

• FACE: dimension 2

• CELL: dimension 3

• ELEMENT: codimension 0

• FACET: codimension 1

E.g., in \(n\) space dimensions, facets are mesh objects of dimension \(n-1\). When \(n=2\) there are no CELL entities.

Nodes can be associated to an existing mesh. Consider the above-defined node

v = NodeId(VERTEX,0)

When it is associated to mesh, it becomes an object of type MeshNode which has coordinate information.

meshv = mesh[v]
print ("type = ", type(meshv))
print ("point = ", meshv.point)

type =  <class 'ngsolve.comp.MeshNode'>
point =  (0.0, 0.0, 0.0)

type(v), type(meshv)   # note the differnt types

(ngsolve.comp.NodeId, ngsolve.comp.MeshNode)

MeshNode objects like meshv can be queried for topology information.

E.g., what are edges connected to the mesh vertex meshv?

meshv.edges

(E0, E1, E2, E3, E4, E5, E6, E8, E9, E11, E12, E13, E14, E16, E22, E24)

c = NodeId(CELL, 1)
meshc = mesh[c]
meshc.faces          # faces of a cell

(F11, F4, F0, F2)

f = NodeId(FACE, 2)
meshf = mesh[f]
meshf.edges, meshf.elements

((E0, E9, E5),
 (<ngsolve.comp.ElementId at 0x7ff5a46a05f0>,
  <ngsolve.comp.ElementId at 0x7ff5a46a05b0>))

ElementId and Ngs_Element

An ElementId is made using a number together with an object like BND or VOL that knows the codimension.

ei = ElementId(BND,0)
type(ei)

ngsolve.comp.ElementId

As with NodeId, we may associate ElementId with a mesh to get an object of type Ngs_Element, which can be queried for topology information.

meshel = mesh[ei]
type(meshel)
print ("type of meshel = \n ", meshel)
print ("vertices =", meshel.vertices)

type of meshel =
  <ngsolve.comp.Ngs_Element object at 0x7ff5a46a0870>
vertices = (V0, V1, V5)

Note that meshel has only three vertices because it is a boundary element. Volume elements have the same type as boundary elements:

elt0 = mesh[ElementId(VOL, 0)]
type(elt0)

ngsolve.comp.Ngs_Element

elt0.vertices, elt0.edges, elt0.facets

((V3, V5, V7, V8), (E16, E22, E25, E13, E15, E21), (F28, F22, F19, F18))

You can also access the same volume element using NodeId:

el0 = mesh[NodeId(ELEMENT, 0)]
type(el0)

ngsolve.comp.MeshNode

el0.vertices, el0.edges

((V3, V5, V7, V8), (E16, E22, E25, E13, E15, E21))

Dofs

Dofs are numbers enumerating the global degrees of freedom of a finite element space. Dofs are associated to mesh entities of the previously described types.

E.g., all dofs of the Lagrange finite element space associated to edges of the mesh can be obtained as follows.

fes = H1(mesh, order=4)
for edge in mesh.edges:
    print ("type = ", type(edge))
    print ("dofs = ", fes.GetDofNrs(edge))

type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (9, 10, 11)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (12, 13, 14)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (15, 16, 17)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (18, 19, 20)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (21, 22, 23)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (24, 25, 26)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (27, 28, 29)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (30, 31, 32)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (33, 34, 35)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (36, 37, 38)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (39, 40, 41)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (42, 43, 44)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (45, 46, 47)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (48, 49, 50)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (51, 52, 53)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (54, 55, 56)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (57, 58, 59)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (60, 61, 62)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (63, 64, 65)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (66, 67, 68)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (69, 70, 71)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (72, 73, 74)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (75, 76, 77)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (78, 79, 80)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (81, 82, 83)
type =  <class 'ngsolve.comp.MeshNode'>
dofs =  (84, 85, 86)

From the output, it is evident that these edge dofs are associated to MeshNode objects.

What about dofs associated to elements?

for el in fes.Elements(VOL):
    print(type(el))
    print (el.dofs)

<class 'ngsolve.comp.FESpaceElement'>
[3, 5, 7, 8, 57, 58, 59, 75, 76, 77, 84, 85, 86, 48, 49, 50, 54, 55, 56, 72, 73, 74, 171, 172, 173, 153, 154, 155, 144, 145, 146, 141, 142, 143, 177]
<class 'ngsolve.comp.FESpaceElement'>
[0, 1, 8, 2, 12, 13, 14, 27, 28, 29, 45, 46, 47, 9, 10, 11, 24, 25, 26, 36, 37, 38, 120, 121, 122, 99, 100, 101, 87, 88, 89, 93, 94, 95, 178]
<class 'ngsolve.comp.FESpaceElement'>
[1, 2, 4, 8, 36, 37, 38, 45, 46, 47, 69, 70, 71, 27, 28, 29, 30, 31, 32, 39, 40, 41, 135, 136, 137, 126, 127, 128, 120, 121, 122, 117, 118, 119, 179]
<class 'ngsolve.comp.FESpaceElement'>
[0, 3, 8, 5, 18, 19, 20, 48, 49, 50, 75, 76, 77, 15, 16, 17, 24, 25, 26, 57, 58, 59, 144, 145, 146, 111, 112, 113, 102, 103, 104, 108, 109, 110, 180]
<class 'ngsolve.comp.FESpaceElement'>
[1, 4, 5, 8, 36, 37, 38, 69, 70, 71, 75, 76, 77, 30, 31, 32, 33, 34, 35, 60, 61, 62, 159, 160, 161, 129, 130, 131, 126, 127, 128, 123, 124, 125, 181]
<class 'ngsolve.comp.FESpaceElement'>
[0, 2, 8, 6, 21, 22, 23, 42, 43, 44, 81, 82, 83, 12, 13, 14, 24, 25, 26, 45, 46, 47, 138, 139, 140, 114, 115, 116, 96, 97, 98, 99, 100, 101, 182]
<class 'ngsolve.comp.FESpaceElement'>
[0, 3, 6, 8, 24, 25, 26, 57, 58, 59, 81, 82, 83, 15, 16, 17, 21, 22, 23, 51, 52, 53, 150, 151, 152, 114, 115, 116, 108, 109, 110, 105, 106, 107, 183]
<class 'ngsolve.comp.FESpaceElement'>
[0, 1, 5, 8, 24, 25, 26, 36, 37, 38, 75, 76, 77, 9, 10, 11, 18, 19, 20, 33, 34, 35, 129, 130, 131, 111, 112, 113, 93, 94, 95, 90, 91, 92, 184]
<class 'ngsolve.comp.FESpaceElement'>
[2, 4, 8, 6, 42, 43, 44, 63, 64, 65, 81, 82, 83, 39, 40, 41, 45, 46, 47, 69, 70, 71, 165, 166, 167, 138, 139, 140, 132, 133, 134, 135, 136, 137, 185]
<class 'ngsolve.comp.FESpaceElement'>
[3, 6, 8, 7, 54, 55, 56, 78, 79, 80, 84, 85, 86, 51, 52, 53, 57, 58, 59, 81, 82, 83, 174, 175, 176, 153, 154, 155, 147, 148, 149, 150, 151, 152, 186]
<class 'ngsolve.comp.FESpaceElement'>
[4, 6, 7, 8, 69, 70, 71, 81, 82, 83, 84, 85, 86, 63, 64, 65, 66, 67, 68, 78, 79, 80, 174, 175, 176, 168, 169, 170, 165, 166, 167, 162, 163, 164, 187]
<class 'ngsolve.comp.FESpaceElement'>
[4, 5, 8, 7, 66, 67, 68, 72, 73, 74, 84, 85, 86, 60, 61, 62, 69, 70, 71, 75, 76, 77, 171, 172, 173, 168, 169, 170, 156, 157, 158, 159, 160, 161, 188]

The type ngsolve.comp.FESpaceElement appearing in the output is derived from Ngs_Element as can be seen from the documentation:

from ngsolve.comp import FESpaceElement
help(FESpaceElement)

Help on class FESpaceElement in module ngsolve.comp:

class FESpaceElement(Ngs_Element)
 |  Method resolution order:
 |      FESpaceElement
 |      Ngs_Element
 |      pybind11_builtins.pybind11_object
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  GetFE(...)
 |      GetFE(self: ngsolve.comp.FESpaceElement) -> ngsolve.fem.FiniteElement
 |
 |      the finite element containing shape functions
 |
 |  GetTrafo(...)
 |      GetTrafo(self: ngsolve.comp.FESpaceElement) -> ngsolve.fem.ElementTransformation
 |
 |      the transformation from reference element to physical element
 |
 |  __init__(self, /, *args, **kwargs)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |
 |  dofs
 |      degrees of freedom of element
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from Ngs_Element:
 |
 |  VB(...)
 |      VB(self: ngsolve.comp.Ngs_Element) -> ngsolve.comp.VorB
 |
 |      VorB of element
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from Ngs_Element:
 |
 |  edges
 |      tuple of global edge numbers
 |
 |  elementnode
 |      inner node, i.e. cell, face or edge node for 3D/2D/1D
 |
 |  faces
 |      tuple of global face numbers
 |
 |  facets
 |      tuple of global face, edge or vertex numbers
 |
 |  index
 |      material or boundary condition index
 |
 |  mat
 |      material or boundary condition label
 |
 |  nr
 |      the element number
 |
 |  type
 |      geometric shape of element
 |
 |  valid
 |      is element valid
 |
 |  vertices
 |      tuple of global vertex numbers
 |
 |  ----------------------------------------------------------------------
 |  Static methods inherited from pybind11_builtins.pybind11_object:
 |
 |  __new__(*args, **kwargs) from pybind11_builtins.pybind11_type
 |      Create and return a new object.  See help(type) for accurate signature.