Contouring unstructured triangular grids

Hello all,

I have some code to generate contour plots of unstructured triangular
grids that I have volunteered to include in matplotlib. Some
discussion of this has occurred in matplotlib-users, see

http://sourceforge.net/mailarchive/forum.php?thread_name=hnbpkq%24o03%241%40dough.gmane.org&forum_name=matplotlib-users

Before I go ahead, I want to check if others think this would be
useful functionality to have in mpl or not.

Thanks,
Ian Thomas

I am a little unclear on what this is and what it is used for. Is this to visualize the triangular grid for things like Finite Element Analysis (FEM) and Computational Fluid Dynamics (CFD)? What kind of format is the data in? Are there any standards for this type of thing? Do you have some example code or images?

Thanks,
-Ben

···

-----Original Message-----
From: Ian Thomas [mailto:ianthomas23@…564…]
Sent: Monday, March 15, 2010 7:01 AM
To: matplotlib-devel@lists.sourceforge.net
Subject: [matplotlib-devel] Contouring unstructured triangular grids

Hello all,

I have some code to generate contour plots of unstructured triangular grids that I have volunteered to include in matplotlib. Some discussion of this has occurred in matplotlib-users, see

http://sourceforge.net/mailarchive/forum.php?thread_name=hnbpkq%24o03%241%40dough.gmane.org&forum_name=matplotlib-users

Before I go ahead, I want to check if others think this would be useful functionality to have in mpl or not.

Thanks,
Ian Thomas

------------------------------------------------------------------------------
Download Intel® Parallel Studio Eval Try the new software tools for yourself. Speed compiling, find bugs proactively, and fine-tune applications for parallel performance.
See why Intel Parallel Studio got high marks during beta.
http://p.sf.net/sfu/intel-sw-dev
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel

I am a little unclear on what this is and what it is used for.

It is used to generate contour plots for data that is defined on
unstructured triangular grids. Currently mpl supports generating
contour plots on regular rectangular grids; if you have an
unstructured grid you have to interpolate it onto a regular grid
before contouring it. Contouring directly from the triangular grid
avoids the need for this interpolation.

Is this to visualize the triangular grid for things like Finite Element Analysis (FEM) and Computational Fluid Dynamics (CFD)?

FEM and CFD are indeed big application areas for unstructured triangular grids.

What kind of format is the data in?

At its simplest, assume you have n irregular spaced points defined by
x and y positions and you wish to contour some field z defined at
those points, where x, y and z are length n numpy arrays or lists.
You will use something like

  tricontour(x, y, z)

A default (Delaunay) triangulation will be constructed and contoured
for your pleasure. Alternatively, if you want to specify your own
triangulation instead you can do that using a indices array of shape
(3, ntri) where ntri is the number of triangles, and indices[:, i]
defines the indices of the three points that the i-th triangle is
composed of.

Are there any standards for this type of thing?

I have no idea. I'm not proposing any file converters for/from any
standard file formats as all the information required for such a grid
is contained in the rather simple x, y and indices arrays.

Do you have some example code or images?

My original posting of September last year includes source code,
examples and images. See

http://sourceforge.net/mailarchive/forum.php?thread_name=4AB3B95B.3090903%40noaa.gov&forum_name=matplotlib-users

This was packaged as a separate python module called 'mpl_tri' which
you can install and play with as much as you desire. I am now
proposing to integrate this into the core of matplotlib as I don't
wish to maintain it as a separate project. I've agreed to add a few
extras, as discussed in

http://sourceforge.net/mailarchive/forum.php?thread_name=hnbpkq%24o03%241%40dough.gmane.org&forum_name=matplotlib-users

Before going ahead with this I wished to ascertain how much interest
there was for this functionality as I don't want to spend time doing
something that isn't wanted or needed.

I hope that answers your questions,
Ian

···

Ben Axelrod <BAxelrod@...817...> wrote:

Ben Axelrod wrote:

I am a little unclear on what this is and what it is used for. Is
this to visualize the triangular grid for things like Finite Element Analysis
(FEM) and Computational Fluid Dynamics (CFD)?

In can be, yes -- though that's really only the part that renders the mesh. Contouring on a triangular mesh can be used to visualize data on such a mesh, as generated by FE or FV codes, but can also be used to contour and visualize data on any arbitrary points -- it's an alternative to interpolating to a regular grid before contouring.

Doing the contouring directly on a triangular mesh is more efficient and bit more robust, particularly near the bounds of the data.

What kind of format is the data in?

It would take simple x,y,z vectors for the contouring routines. There is some discussion of the API in that thread. I'd like to see a standard class that holds a triangular mesh while we're at it -- to make it easier to separate the generation of the mesh from the use of the mesh.

> Are there any standards for this type of thing?

I wish.

There is a netcdf standard for unstructured meshes being developed as we speak, primarily for oceanographic model, but that's all I know of -- and MPL doesn't (and shouldn't) support netcdf directly anyway.

There are common ways to express a triangular mesh -- ans we should use them:

ie. a NX3 array of indexes, where each row in the array is the three indexes into the point coordinate array.

some systems also store indexes to the neighbors of each triangle, etc. I think that should all be hidden in the mesh class.

> Do you have some example code or images?

There were some links in that thread. Ian, maybe you should post them again.

-Chris

···

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

Ian Thomas wrote:

Alternatively, if you want to specify your own
triangulation instead you can do that using a indices array of shape
(3, ntri) where ntri is the number of triangles, and indices[:, i]
defines the indices of the three points that the i-th triangle is
composed of.

That would be better as (ntri, 3), to be compatible with the usual C ordering of numpy arrays.

just a nit,

- Chris

···

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

Well then the meteorologist in me is an overwhelmingly +1 on having
such functionality available. Are you looking at making it possible to
construct a triangulation from the delaunay triangulation that is used
by griddata? (Sorry, I didn't follow the thread that closely.)

Ryan

···

On Mon, Mar 15, 2010 at 12:37 PM, Ian Thomas <ianthomas23@...564...> wrote:

Ben Axelrod <BAxelrod@...817...> wrote:

I am a little unclear on what this is and what it is used for.

It is used to generate contour plots for data that is defined on
unstructured triangular grids. Currently mpl supports generating
contour plots on regular rectangular grids; if you have an
unstructured grid you have to interpolate it onto a regular grid
before contouring it. Contouring directly from the triangular grid
avoids the need for this interpolation.

--
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma

He's using matplotlib.delaunay.

···

On 2010-03-15 13:30 PM, Ryan May wrote:

On Mon, Mar 15, 2010 at 12:37 PM, Ian Thomas<ianthomas23@...564...> wrote:

Ben Axelrod<BAxelrod@...817...> wrote:

I am a little unclear on what this is and what it is used for.

It is used to generate contour plots for data that is defined on
unstructured triangular grids. Currently mpl supports generating
contour plots on regular rectangular grids; if you have an
unstructured grid you have to interpolate it onto a regular grid
before contouring it. Contouring directly from the triangular grid
avoids the need for this interpolation.

Well then the meteorologist in me is an overwhelmingly +1 on having
such functionality available. Are you looking at making it possible to
construct a triangulation from the delaunay triangulation that is used
by griddata? (Sorry, I didn't follow the thread that closely.)

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco

I'm definitely interested, but I defer to Eric as the ultimate arbiter
of all things contour related. I'm not concerned about the C++, we
depend on plenty of that and std C++ won't increase the compile or
distribution burden. I am a little concerned about the "hand wrapped
python" part because of the possibilities of introducing memory leaks
and difficulties it may pose in the matplotlib 2->3 transition. We
use CXX mostly to expose our C++ (eg see src/ft2font.cpp or
src/_backend_agg.cpp) and I'm hopeful that this will ease our
transition to python3 since CXX6 supports python3. It would probably
be beneficial, but by no means required, to use CXX to expose the C++
part to python so take a look if you are inclined.

I'm also a bit concerned by the maintenance aspect, since this is a
fairly sophisticated piece of code. We'd be looking for a commitment
to support it if we include it in mpl.

···

On Mon, Mar 15, 2010 at 12:37 PM, Ian Thomas <ianthomas23@...564...> wrote:

Before going ahead with this I wished to ascertain how much interest
there was for this functionality as I don't want to spend time doing
something that isn't wanted or needed.

Robert Kern wrote:

···

On 2010-03-15 13:30 PM, Ryan May wrote:
  Are you looking at making it possible to

construct a triangulation from the delaunay triangulation that is used
by griddata? (Sorry, I didn't follow the thread that closely.)

He's using matplotlib.delaunay.

right, and the goal is to make it pretty easy to plug in another triangulation routine, or a pre-defined triangulation, if you have one.

-Chris

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

John Hunter wrote:

It would probably
be beneficial, but by no means required, to use CXX to expose the C++
part to python so take a look if you are inclined.

What about Cython -- is any one using Cython in MPL yet?

-Chris

···

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

Christopher Barker wrote:

John Hunter wrote:

It would probably
be beneficial, but by no means required, to use CXX to expose the C++
part to python so take a look if you are inclined.

What about Cython -- is any one using Cython in MPL yet?

No, not yet, but I think we should be looking at moving in that direction. C++ support in cython is being introduced right now, if I understand correctly. I don't know what its limitations are, or how its ease of use and maintenance will compare with CXX.

Eric

···

-Chris

Replying to everyone at the same time: (spot the person in the UK!)

Chris Barker wrote:

That would be better as (ntri, 3), to be compatible with the usual C
ordering of numpy arrays.

Of course.

John Hunter wrote:

I am a little concerned about the "hand wrapped python" part
because of the possibilities of introducing memory leaks and
difficulties it may pose in the matplotlib 2->3 transition. <snip>
It would probably be beneficial, but by no means required, to use
CXX to expose the C++ part to python so take a look if you are
inclined.

I've taken a look at CXX and it looks like a much more pleasant
solution, so I'm happy to use that instead.

I'm also a bit concerned by the maintenance aspect, since this
is a fairly sophisticated piece of code. We'd be looking for a
commitment to support it if we include it in mpl.

Understood. I'm uneasy about making a major commitment at this stage
as I'm new to open source contributions, so how about you ask me again
when I've submitted the code and people have had the chance to
review/correct/admire/destroy it?

Eric Firing wrote:
One possibility would be to develop the tri* functionality at least
initially as a module in lib/mpl_toolkits, like axes_grid and mplot3d;
or in a module in lib/matplotlib. This could still take advantage of
refactoring in contour.py. An advantage is that it would consolidate
the triangle functionality so it would be easier to find, track, and
document.

I like the idea of putting it in its own directory under lib/matplotlib.

What sort of timeline do you have in mind?

I am currently "between jobs", which leaves me plenty of spare time
for interesting projects. I reckon I should have something worthy of
consideration in 2 weeks time.

Ian

Hello all,

Attached is a patch file against svn head to add triangular grid
plotting and contouring. It will need some serious checking/reviewing
before it can be added to MPL. I've tested it on 32 and 64-bit Linux,
but I don't use other operating systems very often.

Most of the new code is in a new directory lib/matplotlib/tri. There
is a Triangulation class to store an unstructured triangular grid for
reuse; the user can either specify the triangles or allow
matplotlib.delaunay to create a Delaunay triangulation. There are
three plotting functions: triplot to plot grid lines and points,
tripcolor to draw a pseudoplot, and tricontour/tricontourf to
calculate and draw contour lines and filled contours. You can either
pass in a Triangulation object to these functions, or pass in the x,y
points, etc and have MPL create a temporary Triangulation for you (for
convenience).

The underlying contouring code is C++, and I've used CXX to access it.

There are some changes to axes.py and pyplot.py to expose the new
functionality, and changes to the build scripts to build the new C++
module. I've also changed contour.py, splitting the previous
ContourSet class into ContourSet and a derived QuadContourSet so there
is separate responsibility for creating contours (QuadContourSet) and
storage/display (ContourSet). ContourSet allows you to specify your
own polygons to draw, which should allow easier extension to other
contouring algorithms. QuadContourSet is now used by contour/contourf
to create contour lines/polygons for a quad grid and the base class
stores and draws them. The benefit of this approach is that
QuadContourSet is a relatively small class, allowing the equivalent
triangular grid contouring class TriContourSet to be similarly small,
and there is little code duplication.

I've added examples to demonstrate the new functionality, including
manually creating your own contours, and also a comparison between
griddata and tricontourf to contour unstructured grids.

Both the C++ and python are documented, but as I'm not familiar with
the relationship between pydoc strings and sphinx, you may need to
guide me further to increase/reduce links between and duplication of
pydocs.

All questions and comments gratefully received.
Ian

patch.diff (137 KB)

Ian Thomas wrote:

Hello all,

Attached is a patch file against svn head to add triangular grid
plotting and contouring. It will need some serious checking/reviewing
before it can be added to MPL. I've tested it on 32 and 64-bit Linux,
but I don't use other operating systems very often.

Most of the new code is in a new directory lib/matplotlib/tri. There
is a Triangulation class to store an unstructured triangular grid for
reuse; the user can either specify the triangles or allow
matplotlib.delaunay to create a Delaunay triangulation. There are
three plotting functions: triplot to plot grid lines and points,
tripcolor to draw a pseudoplot, and tricontour/tricontourf to
calculate and draw contour lines and filled contours. You can either
pass in a Triangulation object to these functions, or pass in the x,y
points, etc and have MPL create a temporary Triangulation for you (for
convenience).

The underlying contouring code is C++, and I've used CXX to access it.

There are some changes to axes.py and pyplot.py to expose the new
functionality, and changes to the build scripts to build the new C++
module. I've also changed contour.py, splitting the previous
ContourSet class into ContourSet and a derived QuadContourSet so there
is separate responsibility for creating contours (QuadContourSet) and
storage/display (ContourSet). ContourSet allows you to specify your
own polygons to draw, which should allow easier extension to other
contouring algorithms. QuadContourSet is now used by contour/contourf
to create contour lines/polygons for a quad grid and the base class
stores and draws them. The benefit of this approach is that
QuadContourSet is a relatively small class, allowing the equivalent
triangular grid contouring class TriContourSet to be similarly small,
and there is little code duplication.

I've added examples to demonstrate the new functionality, including
manually creating your own contours, and also a comparison between
griddata and tricontourf to contour unstructured grids.

Both the C++ and python are documented, but as I'm not familiar with
the relationship between pydoc strings and sphinx, you may need to
guide me further to increase/reduce links between and duplication of
pydocs.

All questions and comments gratefully received.
Ian

Ian,

I've only glanced so far, but one thing that caught my eye is your documentation changes and code regarding the need for simply-connected paths. This is obsolete--mpl now handles multiply-connected paths. If you look in cntr.c, you will find a "reorder()" function which converts the simply-connected paths to multiply-connected paths, which render better, because they don't have the cuts. (Maybe the core routines in cntr.c could be modified so that the multiply-connected paths would be generated directly, so that reorder() would not be needed; but I find those routines very difficult to follow, so writing reorder() was easier.)

A second initial suggestion is that you divide the work into two patches, one of which provides the refactoring of existing code (presumably only in contour.py), and a second which adds the new functionality. This would make reviewing and debugging easier. (Your contour_manual.py could also go in the first patch.)

Thanks for all the good work!

Eric

Eric Firing wrote:

I've only glanced so far, but one thing that caught my eye is your
documentation changes and code regarding the need for simply-connected
paths. This is obsolete--mpl now handles multiply-connected paths.

Thanks for clarifying this, I now understand the 'point kinds' in
cntr.c which I originally thought were for debug purposes. It will
make my code simpler in the end.

A second initial suggestion is that you divide the work into two patches,
one of which provides the refactoring of existing code (presumably only in
contour.py), and a second which adds the new functionality.

Good idea. I've attached the first patch which changes axes.py and
contour.py, plus adds the new example contour_manual.py. My approach
to the refactoring is probably more of a C++ than pythonic way of
thinking with derived classes that override base class methods, but it
avoids changing much of the low level code and avoids code duplication
when tricontour is added in.

Ian

contour_refactor.diff (28 KB)

Ian Thomas wrote:

Eric Firing wrote:

I've only glanced so far, but one thing that caught my eye is your
documentation changes and code regarding the need for simply-connected
paths. This is obsolete--mpl now handles multiply-connected paths.

Thanks for clarifying this, I now understand the 'point kinds' in
cntr.c which I originally thought were for debug purposes. It will
make my code simpler in the end.

A second initial suggestion is that you divide the work into two patches,
one of which provides the refactoring of existing code (presumably only in
contour.py), and a second which adds the new functionality.

Good idea. I've attached the first patch which changes axes.py and
contour.py, plus adds the new example contour_manual.py. My approach
to the refactoring is probably more of a C++ than pythonic way of
thinking with derived classes that override base class methods, but it
avoids changing much of the low level code and avoids code duplication
when tricontour is added in.

Ian

Ian,

I went ahead and committed your refactoring. It seems to work fine. Contour.py was not pretty before the refactoring, and it is still not pretty, but I don't think it is much worse, if at all. I have a vague sense that ideally it would be factored somewhat differently, but I don't seem to have the time or vision to pursue it right now, and maybe if I did, I would still run into a dead end. So, let's proceed with what you have. If inspiration strikes later, I suspect it won't be too hard to modify both contour and tricontour accordingly.

Eric

Eric Firing wrote:

I went ahead and committed your refactoring. It seems to work fine.
Contour.py was not pretty before the refactoring, and it is still not
pretty, but I don't think it is much worse, if at all. I have a vague sense
that ideally it would be factored somewhat differently, but I don't seem to
have the time or vision to pursue it right now, and maybe if I did, I would
still run into a dead end. So, let's proceed with what you have. If
inspiration strikes later, I suspect it won't be too hard to modify both
contour and tricontour accordingly.

Understood.

I've attached the patch for the new triangular grid functionality.
I've altered it to use multiply-connected paths for filled contours
which has made the code a bit simpler overall.

Ian

tri_patch.diff (104 KB)

As nobody has commented on the patch I submitted to add triangular
grid functions, I can only assume that nobody has looked at it. This
is a final friendly reminder that some people considered that it would
be useful and I'd rather not throw away useful work.

Ian

···

On 15 April, Ian Thomas wrote:

I've attached the patch for the new triangular grid functionality.

I've attached the patch for the new triangular grid functionality.

As nobody has commented on the patch I submitted to add triangular
grid functions, I can only assume that nobody has looked at it. This
is a final friendly reminder that some people considered that it would
be useful and I'd rather not throw away useful work.

Ian,

I'm sorry, this has been on my mind, but I have put it off day by day. In addition, I'm not a C++ writer.

I am somewhat inclined to simply commit it, or to give you commit rights, in view not only of the triangular grid work but of the great work you did in fixing problems with cntr.c. Do you have a sourceforge login?

Eric

···

On 05/10/2010 10:09 PM, Ian Thomas wrote:

On 15 April, Ian Thomas wrote:

Ian

------------------------------------------------------------------------------

_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
matplotlib-devel List Signup and Options

Ian Thomas wrote:

···

On 15 April, Ian Thomas wrote:

I've attached the patch for the new triangular grid functionality.

As nobody has commented on the patch I submitted to add triangular
grid functions, I can only assume that nobody has looked at it.

I have NO time to look at it, but I think it's great think to add.

-Chris

--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...