-
Notifications
You must be signed in to change notification settings - Fork 4
PyBCPyTaps
= Python/C Extension: PyTAPS =
== What is PyTAPS? ==
PyTAPS is a set of Python bindings to [http://www.tstt-scidac.org/ ITAPS] (Interoperable Tools for Advanced Petascale Simulations), a set of standardized interfaces for use in simulation code written in C, C++, or FORTRAN. Specifically, ITAPS includes:
- ''iBase'': Common definitions used by all other ITAPS interfaces
- ''iMesh'': Creation, querying, and manipulation of discrete polygonal/polyhedral meshes
- ''iGeom'': Creation, querying, and manipulation of continuous geometry (e.g. CAD)
- ''iRel'': Relations between pairs of other ITAPS interface instances
ITAPS strives to remain implementation-independent, allowing people to write generic simulation code that can be used on any system that implements an ITAPS interface. There are several software packages that implement ITAPS interfaces, including [http://trac.mcs.anl.gov/projects/ITAPS/wiki/MOAB MOAB], [http://tetra.mech.ubc.ca/GRUMMP/ GRUMMP], [http://www.scorec.rpi.edu/FMDB FMDB], [http://www.emsl.pnl.gov:2080/nwgrid/ NWGrid], [http://frontier.ams.sunysb.edu/news/news.php Frontier], [http://trac.mcs.anl.gov/projects/ITAPS/wiki/CGM CGM], and [http://trac.mcs.anl.gov/projects/ITAPS/wiki/Lasso Lasso]. PyTAPS wraps all of this up and makes it available to Python.
{{{ #!div style="border: 1px solid #d7d7d7; margin: 1em 1.75em; padding: .25em; overflow: auto;" '''Aside: Hands-on Examples''[[BR]] Unlike many of the other talks, this one does not include any hands-on examples. To use all of the features in PyTAPS, users would need to install implementations for all the ITAPS interfaces in addition to PyTAPS itself. On my system, this translates to a total of five different packages: PyTAPS, MOAB, CGM, CUBIT (a geometry backend for CGM), and Lasso.
For those with access to these libraries, PyTAPS is available via {{{easy_install}}} or on [http://pypi.python.org/pypi/PyTAPS PyPI]. }}}
== Why Does it Matter? ==
Well, for the same reason that we'd want to use Python for any kind of scientific computing: it's easier! 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. We'll start by writing this in C, using iMesh.
{{{ #!CodeExample #!c #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: %sn", 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, %fn", x, y, z);
} printf("n");
}
free(adj_ents); free(offsets); free(faces);
iMesh_dtor(mesh, &err); if(err != 0)
handleError(mesh);return 0;
Now, let's try the same thing in Python, using PyTAPS:
{{{ #!CodeExample #!python from itaps import iBase, iMesh
mesh = iMesh.Mesh() mesh.load("mesh.vtk")
faces = mesh.getEntities(iBase.Type.face, iMesh.Topology.all) 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
}}}
It turns out that the source code for the C version of our program is over '''5 times''' the size (in bytes) of that of our Python program! Because of the sheer volume 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!
== Using iMesh ==
As mentioned above, iMesh is an interface for dealing with discrete polygonal/polyhedral meshes. iMesh has a notion of three basic types of data: ''entities'' (individual elements in the mesh, such as vertices or tetrahedra), ''entity sets'' (collections of entities), and ''tags'' (labels for storing data on entities or entity sets).
In the examples below, we'll assume that we've already typed the following: {{{ #!CodeExample #!python from itaps import iBase, iMesh mesh = iMesh.Mesh() }}}
=== Creating Entities ===
Once we have our mesh object, we need to get some data into it before we can do anything. We can either load the data from a file, as we did above, or create entities on our own. The simplest entity to create is a vertex, specified by a 3-vector. We can also create several vertices at once by providing a list of points:
{{{ #!CodeExample #!python vertex = mesh.createVtx([1, 2, 3]) vertices = mesh.createVtx(1, 1, 1], [2, 2, 2], [3, 3, 3,
iBase.StorageOrder.interleaved)
}}}
With vertices, we can create higher-order elements, such as quadrilaterals, triangles, tetrahedrons, and so forth. Let's create four new vertices and use them to create a quadrilateral:
{{{ #!CodeExample #!python vertices = mesh.createVtx(0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1,
iBase.StorageOrder.interleaved)
quad = mesh.createEnt(iMesh.Topology.quadrilateral, vertices) }}}
=== Getting/Querying Entities ===
As we saw above, we can easily retrieve entities of a particular ''Type'' (dimension) and/or ''Topology'' (shape). Likewise, we can retrieve the coordinates of any vertex: {{{ #!CodeExample #!python quads = mesh.getEntities(iBase.Type.all, iMesh.Topology.quadrilateral) verts = mesh.getEntities(iBase.Type.vertex, iMesh.Topology.all) coords = mesh.getVtxCoords(verts, iBase.StorageOrder.interleaved) }}}
There are several other basic things we can do with these entities. For instances, with vertices, we can change their coordinates:
{{{ #!CodeExample #!python verts = mesh.getEntities(iBase.Type.vertex, iMesh.Topology.all) mesh.setVtxCoords(verts, 0, 0, 0*len(verts), iBase.StorageOrder.interleaved) }}}
=== Sets ===
Entity sets are a way of organizing entities in ITAPS. As you may expect, they are a collection of entities in the mesh. Like mathematical sets (and Python sets), we can perform the usual set operations on them: union, intersection, and difference. The entire mesh is contained in the "root set", referenced by {{{mesh.rootSet}}}. This is so common that in PyTAPS, the "{{{.rootSet}}}" part is optional. We can see this above when using {{{getEntities}}}. This method actually operates on a set, but we can make it operate on the root set by typing {{{mesh.getEntities(...)}}}.
In addition, we can create and manipulate our own sets: {{{ #!CodeExample #!python set1 = mesh.createEntSet(ordered=True) set2 = mesh.createEntSet(ordered=True)
union = set1 | set2 intersection = set1 & set2 difference = set1 - set2
mesh.destroyEntSet(set1) }}}
=== Tags ===
Tags allow you to apply data to entities or entity sets. This data could be virtually anything, from floating point values describing temperature at each vertex, to strings providing IDs for important entities. When creating a tag, we must specify what type of data the tag will hold, and how many units of that data will be present for each entity with that tag. We can then set data to (and get data from) our entities and sets:
{{{ #!CodeExample #!python temp = mesh.createTag('temperature', size=1, type='d') # 'd' = C double (or a Python float)
temp.setData(verts, [40.0, 41.2, 45.4, 50.8, 44.5]) temp_values = temp.getData(verts)
mesh.destroyTag(temp) }}}
== Using iGeom ==
iGeom is very similar to iMesh in many ways. Though it operates on continuous (CAD-like) geometry instead of discrete meshes, most of the same concepts are present. iGeom instances consist of entities, entity sets, and tags, and the interfaces for entity sets and tags are almost identical. As such, we'll only cover the important differences in the iGeom interface.
In the examples below, we'll assume that we've already typed the following: {{{ #!CodeExample #!python from itaps import iBase, iGeom geom = iGeom.Geom() }}}
=== Creating Entities ===
Rather than creating individual vertices as we would in iMesh, creation of entities in iGeom is performed with basic three-dimensional shapes, such as prisms or cones. This means that, unlike in iMesh, where we explicitly create each entity, iGeom may implicitly create many entities. For instance, if we create a cube, iGeom will create the cube itself, 6 faces, 12 edges, and 8 vertices:
{{{ #!python >>> brick = geom.createBrick(1, 1, 1) >>> print geom.getNumOfType(iBase.Type.all) 27 >>> print geom.getNumOfType(iBase.Type.edge) 12 }}}
{{{ #!CodeExample #!python torus = geom.createTorus(major_rad=4, minor_rad=1) }}}
=== Transforming Entities ===
Additionally, iGeom supports a variety of transformation functions to manipulate our entities. These are especially necessary since all the creation functions create their entities at the origin.
{{{ #!CodeExample #!python
geom.reflectEnt(brick, [0, 1, 1]) geom.moveEnt(brick, [5, 0, 0]) geom.intersectEnts(brick, torus) }}}
== The Future ==
PyTAPS is still under active development and is finalizing its implementation of the iRel interface, which will allow users to create a relation between an iGeom instance and an iMesh instance, useful when the mesh was generated from that geometry.
PyTAPS is also looking towards interactive mesh manipulation. Since it uses Python, we are able to query and manipulate meshes in the Python interactive interpreter. However, this has the unfortunate downside that we can't ''see'' the changes we're making. Thus, we have been pushing to get a visualization application that works directly with the ITAPS interfaces and allows Python scripting. The vis tool VisIT works with ITAPS (but only via files) and scriptable via Python, so we hope to be able to take advantage of this to create an interactive mesh exploration tool using PyTAPS.