C.1 The simple algorithm

I am afraid that my routine for contouring an irregular area is a bit long to present in this book but the principles are very simple and in this instance are easier to describe than to demonstrate - in contrast to the normal situation in life when it’s easier to show someone than to describe what they need to do.

Contouring may be done with line contours or coloured bands. In both cases an essential prerequisite is to divide the area to be contoured into triangles within which the variation of the quantity to be contoured varies linearly or approximately so. Some algorithms in analysis operate using triangles in the first instance and others with shapes that are very easily divided into triangles.

You then decide what the contour interval is to be, and given the total range of field values in the problem, how they are distributed. If you have both positive and negative values, then the contours are distributed on both sides of zero. If you only have positive (or negative) values then the question is whether or not to divide the problem into a finite number of intervals with irrational contour levels or whether to again take a base of zero and then accept, for a given contour interval, however many contours there should be.

The first step with each triangle is to examine the field values at each corner or node. If they are all the same, or all three lie within one contour band then with line contours nothing needs to be drawn, whereas with area fill the colour selected needs to be appropriate to the contour band, and the triangle is drawn filled with that colour.

The second case is a triangle where two nodes have the same value, but the third is different either higher or lower. Initially, fill the triangle with the colour appropriate to the values that are the same, and then proceed to draw successive triangles using the remaining apex node and two points established by linear interpolation along the sides between that remaining apex node and the side where the two points are the same. This procedure may require working to successively higher-level contours, or lower. Instead of requiring area fill a line contour is simply a matter of drawing the line.

I have assumed that you are working with the ClearWin+ functions and a Windows printer so that the painter’s algorithm applies, that is you only see the last colour and the later graphics areas obscure the earlier. It slightly more complicated if the painter’s algorithm does not apply and that earlier colouring shows through.

The third case is where all three nodes have different values. The simplest way of treating this is to divide the triangle into two further triangles with a point on the side between the lowest and highest values picked to match the node value that is intermediate. Each of the triangles then becomes case 2.

A version of this algorithm (more or less, and for line contours) was published in Byte magazine, but I see in the development notes for a graphics program of my own that I had developed it independently somewhat earlier, although from experience of the time it takes to go from a paper draft in its submission to a journal to publishing then I am not comparing like with like. However, I certainly used it to generate the results published in a paper (Pugh & Bromhead, 1985) although they were there redrawn by a draughtsman from a computer plot both for an initial technical report and for the paper, and one of the figures appears in a book of mine published in 1985.

Figure C1. The sequence of triangles to show different contour bands. Left: the individual triangles, Right: triangles drawn over each other. The apex with the value higher or lower than the other two is on the right in each case.

C.2 Using a plotter

Having used various pen plotters when they were available attached to mainframes or personal computers I found several things that may still be relevant or may not. One of the things was that much of the time taken to produce a plot was taken up by pen changes and that plots could be speeded considerably by doing multiple passes through the contouring, finishing one colour first before another. This is particularly the case where, for example, the zero level contour is in a different colour, positive and negative contours are coloured differently, or periodically, contours are given different emphasis for example when contouring topography where the contour interval might be say 10 m and every 50 or 100 m the contour is drawn in a different colour or weight. Plotters also slowed down dramatically when dotted or dashed lines needed to be drawn. The point is that even the slowest computer is dramatically faster than even a very fast pen plotter. Pen plotters also are very poor at drawing filled areas and because they use wet ink they have all the defects in causing paper stretch that are found with the use of highly saturated colours with inkjet printers. These things are not ClearWin+ issues but are things that you (might) need to worry about within your Fortran code.

C.3 Labelling the contour values

I mention this because it is an issue that arises very quickly with any contour graphic.

If you can be certain that the contour lines always finish at the boundary of the graphic, then the annotation can be applied there. I have a particular application where this is so and the contour lines always finish at the sides and/or bottom of the plot area and I have chosen to label them along the bottom and on the right-hand side, with the right-hand side selected because that there is a known position to start the text.

With coloured bands, labelling on the body of the graphic is not an effective way to do things and a key or legend is better. Assuming that one is operating with ClearWin+ graphics, then to ensure that the key does not overwrite any part of the graphic, it is better to assume that one part of the available drawing surface is dedicated to the key and that the remainder is dedicated to the plot itself.

When using line contours, labelling is best done on specific contours. One method is to determine when the individual segments of a particular contour are long enough to be omitted and replaced with the contour label, and another method is to simply overwrite the contour, perhaps using a different colour for the label, or blanking out the line by overwriting it with the background colour before writing the label. This task is trivial for a skilled draughtsman, but is rather difficult to program, especially as the text values for specific contour levels may vary in length. A method that I have found simplifies matters is to imagine that there are one or more virtual (i.e. invisible and not plotted) lines through the graphic, perhaps a vertical line and a horizontal line, and to only label contours where they pass through the virtual lines. Labels can be a single character if a key is provided, and labelling with letters enables you to (for example) associate uppercase with positive values and lower case with negative. Short contour labels are in any case very desirable on screen as the resolution is likely to be limited.

C.4 A slight warning (or two)

Triangulation contouring can be somewhat misleading in the vicinity of isolated high (or low) points, as the effect of the high point is distributed over a larger area than may actually be the case. It is a problem that occurs when plotting topographic surveys where the elevations are measured as a discrete number of points. It can also be a problem when the triangulation is automatic and triangles with a very long and narrow aspect ratio are thereby generated. Misleading contours are also generated if they are smoothed over a discontinuity as for example in plotting strata contours on a geological map without taking into account a fault or faults, especially when the existence of such a fault is unknown.

Pen plotters were notorious for running out of ink at critical points in a generation of a big and complicated plot and therefore requiring the whole lot to be plotted a second time. In the past I found it useful to be able to rerun the plot from a particular point where the pen had run out of ink rather than to run a complete second drawing.

C.5 Finite elements and smoothing

I have experiences in this field that date back a long way. Some software takes the approach that the problem is divided into simple triangles anyway, with field values determined at nodes. A variant of the triangular element approach was to divide the problem into quadrilaterals, each of which was formed by four triangles. In both of these cases, the contours can only be smooth(er) if a finer mesh of elements is used. Where quantities are determined inside the elements, then contours that are continuous from element to element can be obtained but with the loss of some accuracy by averaging at the nodes, including, if necessary, any appropriate extrapolation from the inside of an element to its nodes.

My own particular preference in the finite element field is to use the isoparametric elements that are basically quadrilaterals using 4 (linear), 8 (quadratic) or 12 (cubic) nodes. In much of my work I have used Serendipity shape functions, (see Zienkiewicz & Taylor), but there may be some merits in using Lagrangian shape functions at the cost of more nodes, some of which are internal and not simply distributed around the periphery of each element.

Describing only serendipity elements, eight and 12-noded variants may need to be plotted with linear segments around the periphery and again - provided enough elements are specified - this barely shows. The four- and eight-node elements produce reasonable contours if the triangles are obtained from a central point and segments of the periphery. However, as the field value defined at nodes has respectively a quadratic or cubic distribution, it is possible to determine field values at internal points. I have found it extremely useful in a way that works with quadratic or cubic variation to divide each element into a number of smaller quadrilaterals with virtual nodes that are interpolated, and a 6 by 6 grid guarantees that using either element the actual nodes will coincide with at least some of the virtual (sub-element) nodes. The figures below illustrate this subdivision.

Figure C.2 Top row: Elements of the Serendipity family, which have simple polynomial shape functions, Bottom row: the equivalent Lagrangian elements, these having internal as well as peripheral nodes. Left to right: first, second and third order elements.

Figure C.3 A 6x6 subdivision of each element including a subsequent 4-triangle subdivision of each sub-element works with Serendipity or Lagrangian elements of 1^st, 2^nd and 3^rd order, with the appropriate parent element nodes fitting exactly onto the sub-element nodes. Field values obtained at the nodes can be interpolated to the sub-element nodes, and whereas the central node for division of a sub-element into triangles could be obtained by straightforward averaging from the four corners, it too could be obtained via shape function interpolation from the parent element.

FORTRAN and the ART of Windows Programming, Copyright © Eddie Bromhead, 2023.