6.2.2 Data Output

For the direct output of graphics images, the PNG file format has been chosen, because it is widely supported (also by all newer WWW clients), generates very small, highly compressed files, it is not encumbered by any patents (like the GIF file format), and the reference library libpng [66] is freely available.

However, first the magnetization, which is defined on the nodes of the unstructured finite element mesh, has to be interpolated on a regular mesh in the desired slice plane. The slice plane is defined in the configuration file in the form

$$ax+by+cz=d \quad, \quad \boldsymbol{n} = \left( \begin{array}{c} a \\ b \\ c \\ \end{array} \right)$$ (6.1)

with the normal vector $\boldsymbol{n}$. Then all elements, which are cut by the slice plane have to be found. In order to simplify this task the whole model (simply the Cartesian coordinates of the vertices) is translated and rotated in such a way, that the slice plane coincides with the $x$-$y$-plane. Then all elements which are cut by the $x$-$y$-plane - this is easily determined using the $z$-coordinates of the vertices - are tagged, and the maximum dimensions (the ``bounding rectangle'') of the finite element model in the slice plane are determined. For a given image size we can now calculate the Cartesian coordinates $\boldsymbol{P}$ of each pixel $P$, which then has to be assigned to one of the tagged elements, within which the pixel is located. This ``point location problem'', which is a standard problem in geographic information systems and computer-aided design and engineering [74], is the most time consuming part. However, the number of tagged elements (which are cut by the slice plane) is usually far smaller than the total number of elements, and so we can use a brute force yet clever approach, which takes advantage of some peculiarities of our specific problem.

In principle, for each pixel we have to search all tagged elements and determine if the pixel is inside or outside. But we can speed up the search by skipping those elements, for which the distance of the pixel from the first node is larger than the maximum edge length of any tetrahedron.

Figure 6.2: Barycentric coordinates in a regular tetrahedron.

If a suitable element is found we calculate the barycentric coordinates $\boldsymbol{B}_P=(\alpha, \beta, \gamma, \delta)$ of the pixel. The barycentric coordinates are defined as the relative position of $P$ inside the tetrahedron $ABCD$. In a ``regular'' tetrahedron (which has three edges coinciding with the coordinate axes) these correspond to the relative coordinates as indicated in Fig. 6.2.

$$\boldsymbol{P} = \alpha \boldsymbol{A} + \beta \boldsymbol{B... ...\delta \boldsymbol{D} \quad,\quad \delta=1-\alpha-\beta-\gamma$$ (6.2)

If we solve for $(\alpha, \beta, \gamma)$ we get a simple linear system of equations

$$\left( \begin{array}{c} \alpha \\ \beta \\ \gamma ... ...bol{D} \right)^{-1} \cdot (\boldsymbol{P} - \boldsymbol{D})$$ (6.3)

The point $P$ is inside the tetrahedron if and only if

$$\alpha, \beta, \gamma, \delta \ge 0 \quad,$$ (6.4)

which is equivalent to

$$\alpha, \beta, \gamma \ge 0 \quad \mathrm{and} \quad \alpha + \beta + \gamma \le 1 \quad.$$ (6.5)

The barycentric coordinates $(\alpha, \beta, \gamma, \delta)$ give the weight, by which the data values at the corresponding vertices of the tetrahedron contribute to the linear interpolation at $P$. Thus, for some data $x_i$ defined on the vertices of the tetrahedron, we find the linear interpolation at $P$ with

$$x_P=\alpha x_A + \beta x_B + \gamma x_C + \delta x_D \quad.$$ (6.6)

As a result we can assemble an interpolation matrix, which calculates the magnetization for all pixels of the image. Let $\boldsymbol{M}_\mathrm{FE}$ denote the PETSc vector of one Cartesian component of the magnetization of all nodes of the finite element mesh and $\boldsymbol{M}_\mathrm{ip}$ the PETSc vector of interpolated values, then we have

$$\boldsymbol{M}_\mathrm{ip} = A_\mathrm{ip} \boldsymbol{M}_\mathrm{FE}$$ (6.7)

with the $(m \times n)$ interpolation matrix ($m=a\times b$ rows for images with $a \times b$ pixels and $n$ columns for a finite element mesh with $n$ nodes). The matrix elements $A_{ij,\mathrm{ip}}$ are given by the barycentric coordinate of pixel $P_i$ corresponding to node with the global index $j$ of the finite element, in which the pixel is located. This matrix is very sparse (cf. Fig. 6.3 for the sparsity pattern of the nanodot model), as it has only four entries per row (the four barycentric coordinates of a given pixel).

Figure 6.3: Sparsity pattern of the interpolation matrix for graphics output.

Since the interpolation matrix depends only on the triangulation, we can calculate it during the initialization phase of the program. Whenever a snapshot of one magnetization component should be stored to disk, a simple matrix-vector multiplication gives the interpolated values at the pixel positions. Then the data are color coded using the scheme shown in Fig. 6.4. Finally, the PNG library [66] is called to convert these data into a PNG graphics file.

Figure 6.4: RGB color intensities for data encoding and resulting color map.

In addition to snapshots of the magnetization it is often desirable to measure some data along a probing line through the model. This feature has been implemented based on the algorithm described above. For the probing line we translate and rotate the model to make the probing line coincide with the $x$-axis. Then we tag all elements, which are cut by the $x$-axis and determine the barycentric coordinates of every point on the sampling line. We finally end up with another interpolation matrix, which calculates the data values along the probing line. These can then be stored in a data file on disk.

Finally, the complete set of data, which includes the magnetization, the magnetic scalar potential, and the local fields are stored in files in UCD format [75]. This file format can be read by many advanced visualization packages like AVS [76] and MicroAVS [77]. However, the (human-readable) ASCII UCD format generates huge files, because the whole finite element mesh is stored in each file. Thus, two measures have been taken to tackle this problem.

First, the finite element mesh is stored in a single separate file in UCD format during the initialization phase. The data files, which are generated during the simulation are then stored without the mesh definition, which saves about 50% of disk space. Secondly, the data files are compressed using the zlib-library [65] in gzip-compatible format, which shrinks the data to about 30% of their initial size. Thus, we end up with data files, which require only 15% of the uncompressed storage space. To restore the data in proper UCD format, a simple shell script uncompresses the data and merges them with the mesh data. Alternatively, there is also a binary UCD format, which gives smaller files than the standard ASCII format, but it is not as widely supported and may lead to compatiblity problems, when the files are transferred to machines with different hardware architecture and operating system.