Grids in Bempp

This tutorial demonstrates basic features of dealing with grids in Bempp. Simple grids can be easily created using built-in commands. More complicated grids can be imported in the Gmsh format.

In order for the commands in this tutorial to work, Gmsh must be installed and the command gmsh must be available in the path.

Creation of basic grid objects

Let us create our first grid, a simple regular sphere.

import bempp.api
import numpy
as np
grid = bempp.api.shapes.regular_sphere(3)

The command regular_sphere creates a sphere by refining a base octahedron. The number of elements in the sphere is given by * 4**n, where n is the refinement level. We can plot the grid with the following command.

grid.plot()


This uses Gmsh to plot the grid externally. Another way to create a sphere is by specifying the width of the elements. The command

grid = bempp.api.shapes.sphere(h=0.1)

will create an unstructured spherical grid with a grid size of roughly 0.1.

The shapes module contains functions for a number of shapes, including spheres, ellipsoids, cubes and the Nasa almond shape. Sometimes, it is desired to create a regular structured 2D grid (such as a screen). For this Bempp offers the function bempp.api.structured_grid.

Creating grids from connectivity data

If the shape you require is not provided, a grid may be created using connectivity data: arrays containing the nodes (or vertices) and the element definitions.

The vertices may be defined as follows.

vertices = np.array([[0,1,1,0],
                     [0,0,1,1],
                     [0,0,0,0]])

The array vertices contains the (x,y,z) coordinates of the four vertices of the unit square in the xy plane. Note that each column is the coordinate of one point.

We now define two elements by specifying how the vertices are connected.

elements = np.array([[0,1],
                     [1,2],
                     [3,3]])

The first element connects the vertices 0, 1 and 3. The second element connects the vertices 1, 2 and 3. As before, each column is an element.

To create a grid from these two elements we simply call the following command.

grid = bempp.api.grid_from_element_data(vertices, elements)

Bempp assumes that each element is defined such that the normal direction obtained with the right-hand rule is outward pointing. Elements with inward pointing normals will cause errors in computations. Normal directions can be visually checked for example by loading a mesh in Gmsh and displaying the normals. Also, it is not guaranteed that elements are stored in the grid object using the same numbering as during insertion. Later, we will explain this in detail using the methods grid.vertex_insertion_index and grid.element_insertion_index.

Importing grids from files

Grids can be imported from external files. Bempp natively supports the Gmsh v2.2 msh format (only ASCII and not binary). The Gmsh documentation gives details of this format. Gmsh can easily convert from various formats into .msh. The following command imports a grid into Bempp from an external file.

grid = bempp.api.import_grid('my_grid.msh')

Please note that it is important to choose the correct file ending .msh. Bempp uses the extension to recognize the file format.

More details of this interface are given in the Gmsh interface tutorial.

Iterating through grids

For many Bempp usage scenarios, the internals of the grid object are not important. However, sometimes it may be useful to iterate through a grid and retrieve information from individual elements. Internally, Bempp uses Dune-Grid to represent a grid. The basic objects in Dune are entities of a given codimension. For surface meshes in Bempp this translates as follows:

  • Codim-0 entities: Elements (triangles) of the mesh
  • Codim-1 entities: Edges of the mesh
  • Codim-2 entities: Vertices of the mesh

For example, in order to show the number of elements in the sphere mesh above the following command can be used.

# Create the grid
import bempp.api
grid = bempp.api.shapes.regular_sphere(3)

# Print out the number of elements
number_of_elements = grid.leaf_view.entity_count(0)
print("The grid has {0} elements.".format(number_of_elements))

The grid has 512 elements.

In order to obtain all vertices and elements of a mesh, the following commands can be used.

vertices = grid.leaf_view.vertices
elements = grid.leaf_view.elements

We can also iterate through entities and obtain geometric information about them. The following command stores references to all elements in a Python array.

elements = list(grid.leaf_view.entity_iterator(0))

Alternatively, the iterator may be directly used in a for loop.

for element in grid.leaf_view.entity_iterator(0):
    pass

Let us now print out the corners of the first element in the list above and the area of the corresponding triangle.

corners = elements[0].geometry.corners
area = elements[0].geometry.volume
print("Corners: {0}".format(corners))
print("Area: {0}".format(area))
Corners: [[ 0. 0. 0.19509032]
[ 1. 0.98078528 0.98078528]
[ 0. 0.19509032 0. ]]
Area: 0.0192138328039

Corners are always interpreted column-wise. Hence, the first corner of the element is [0, 1, 0]. The volume attribute depends on the entity. For elements it gives the area and for edges it gives the length of an edge.

Indexing in Dune is slightly more complicated. The above order from the iterator is not guaranteed to agree with the internal indices of the elements. To find out the index of an element or other entity type one can query an IndexSet object.

index_set = grid.leaf_view.index_set()
index = index_set.entity_index(elements[0])
print("The element index is {0}.".format(index))
The element index is 0.

In this case the index of the first element returned by the iterator is indeed 0. However, this is dependent on the implementation of the underlying grid manager and not guaranteed.

Furthermore, the indices from the IndexSet do not need to agree with the order in which elements and vertices were entered into the grid. However, this information is often needed to associate given physical data with mesh entities. For this case the functions grid.element_insertion_index and grid.vertex_insertion_index are provided.

The insertion index of the first element in the elements list is given as follows.

insertion_index = grid.element_insertion_index(elements[0])
print("Insertion index of first element: {0}.".format(insertion_index))
Insertion index of first element: 0

Again, in this case it agrees with the ordering returned by the iterator, but this behavior is not guaranteed and indices should always be computed using an IndexSet or the insertion_index methods.