Why Use PyTAPS

Programmers already familiar with using the ITAPS interfaces in C or Fortran may be wondering why they should use PyTAPS instead. Well, for the same reason that one would want to use Python for any kind of scientific computing: it’s easier!

An Example

For a motivating example, let’s consider an extremely simple program. We want to load in a VTK mesh containing several quadrilateral faces, and then print out the coordinates of each point on each quad, like so:

0.0, 0.0, 0.0
0.0, 1.0, 0.0
1.0, 1.0, 0.0
1.0, 0.0, 0.0

0.0, 0.0, 1.0
0.0, 1.0, 1.0
1.0, 1.0, 1.0
1.0, 0.0, 1.0

0.0, 0.0, 2.0
0.0, 1.0, 2.0
1.0, 1.0, 2.0
1.0, 0.0, 2.0

C Version

We’ll start by writing this in C, using iMesh.

#include <iBase.h>
#include <iMesh.h>
#include <stdio.h>
#include <stdlib.h>

void handleError(iMesh_Instance mesh)
{
    int err;
    char error[120];
    iMesh_getDescription(mesh, error, &err, sizeof(error));
    printf("Error: %s\n", error);
    exit(1);
}

int main()
{
    iMesh_Instance mesh;
    int err;

    iMesh_newMesh("", &mesh, &err, 0);
    if(err != 0)
        handleError(mesh);

    iMesh_load(mesh, NULL, "mesh.vtk", "", &err, 8, 0);
    if(err != 0)
        handleError(mesh);

    iBase_EntityHandle *faces = NULL;
    int faces_allocated = 0;
    int faces_size;

    iMesh_getEntities(mesh, NULL, iBase_FACE, iMesh_ALL_TOPOLOGIES,
                      &faces, &faces_allocated, &faces_size, &err);
    if(err != 0)
        handleError(mesh);

    iBase_EntityHandle *adj_ents = NULL;
    int adj_allocated = 0;
    int adj_size;
    int *offsets = NULL;
    int offsets_allocated = 0;
    int offsets_size;

    iMesh_getEntArrAdj(mesh, faces, faces_size, iBase_VERTEX,
                       &adj_ents, &adj_allocated, &adj_size,
                       &offsets, &offsets_allocated, &offsets_size,
                       &err);
    if(err != 0)
        handleError(mesh);

    int i,j;
    for(i=0; i<faces_size; i++)
    {
        for(j=offsets[i]; j<offsets[i+1]; j++)
        {
            double x, y, z;
            iMesh_getVtxCoord(mesh, adj_ents[j], &x, &y, &z, &err);
            if(err != 0)
                handleError(mesh);
            printf("%f, %f, %f\n", x, y, z);
        }
        printf("\n");
    }

    free(adj_ents);
    free(offsets);
    free(faces);

    iMesh_dtor(mesh, &err);
    if(err != 0)
        handleError(mesh);

    return 0;
}

Python Version

Now, let’s try the same thing in Python, using PyTAPS.

from itaps import iBase, iMesh

mesh = iMesh.Mesh()
mesh.load("mesh.vtk")

faces = mesh.getEntities(iBase.Type.face)
adj = mesh.getEntAdj(faces, iBase.Type.vertex)

for i in adj:
    for j in i:
        x, y, z = mesh.getVtxCoords(j)
        print "%f, %f, %f" % (x, y, z)
    print

Conclusions

It turns out that the source code for the C version of this program is over 5 times the size (in bytes) of that of the Python version! Because of the sheer amount of code necessary to work with ITAPS in C, we end up spending much more time typing and much less time getting to the interesting bits. Besides that, having more code just leaves you with more places for bugs to hide! In fact, the C version already has a bug waiting to bite us: if an error occurs midway through the program, we don’t remember to free() the already-allocated arrays. In this case, we’re lucky because the operating system will reclaim that memory when we call exit(1), but imagine if this were part of a program that loaded many large meshes and tried to continue if one failed!