Plotting and much more

We describe here some aspects of the graphics library available within the castor project. It is based on the well-known VTK library which must be installed by the user. Depending on the operating system, binary files may be directly available. However, one should pay attention to the fact that it is only compatible with version 9.x.

After installing VTK (see Installation), the user only needs to include the castor/graphics.hpp header file.

The path is very close to the one chosen by Matlab or Matplotlib (for the Python language). First, the user creates a figure which is a container which will hold the plots. Then,he adds one or multiple plots to the figure. Finally, he asks the system to display the figures (one window per figure will be created) with the drawnow(...) function. Each call to drawnow displays all the figure objects which have been defined since the last call.

Warning

Each call to drawnow is blocking, meaning that it will suspend the execution of the program until all the opened figure are closed. This behavior cannot be modified as it is due to VTK.

On Windows and linux, when the opened figure are closed, the program stops. A workaround consists in to call the drawnow` function at the end of the program.

The graphics library features 2D/3D customizable plotting and basic triangular/tetrahedral mesh generation. We detail the use of some of these features below. The demo files may be found in the demo/demo_graphics/ subfolder.

It is possible to interact with the figures in multiple ways:

  • “right-click and hold” : rotate the content in the direction of the mouse-pointer for 3D figures. For 2D plots, the user can use it to drag the plots in the direction of choice.

  • “middle-click and hold”: translate the content in the direction of the mouse-pointer.

  • “left-click and hold”: zoom-in if the pointer is in the upper window, otherwise zoom-out.

  • “key-pressed e or q”: close window. Warning: in some systems, it closes all windows.

  • “key-pressed r”: resets the view.

  • “key-pressed w”: switches to wireframe view when applicable.

  • “key-pressed s”: switches to surface view when applicable.

We advise the user to experiment with the different features mentioned above. The summary of all the available functions is available at Basic plot, Graphical input/output, Graphical tools, Mesh management and Mesh plot.

Basic 2D plotting

In this example, we will plot a sine and a square root function with different styles but on the same display. First, we initialize some data and we create a figure object. The full code is available at demo/demo_graphics/basic.cpp.

matrix<> X = linspace(0,10,100);
matrix<> Y = cos(X);
figure fig;

Then, we add our first plot as a red solid line with diamond markers with the correct legend.

plot(fig,X,Y,{"r-d"},{"sin(x)"});

Now we create a second set of data, reusing the already-defined variables. This is possible because the inputs are copied. This time, we plot it as a simple dotted blue line.

X = linspace(-2,2,30);
Y = sqrt(X);
plot(fig,X,Y,{"b"},{"sqrt(x)"});
// display
drawnow(fig);

The figure should look like this:

_images/basic.png

Remark: Currently, it is not possible to save the output to a graphic file.

Surface plot

Surface plotting consists in plotting a surface defined by the equation Z = f(X,Y). First we create the grid (X,Y). The full code is available at demo/demo_graphics/surface_plot.cpp.

matrix<> X,Y;
std::tie(X,Y) = meshgrid(linspace(-M_PI,M_PI,100));

Then, we create the surface which we want to plot, create a figure, adjust the color axis and display everything.

matrix<> Z = 2*sin(X)/X * sin(Y)/Y;
// create the figure
figure fig;
caxis(fig,{-1,1}); // scaled color in the range [-1,1]
mesh(fig,X,Y,Z);
// display
drawnow(fig);

The result is a 3-dimensional plot which should look like this :

_images/surface_plot_wireframe.png

TIP: It is possible to switch to a full surface plot by pressing the s key and switch back to a wireframe display by pressing the w key.

_images/surface_plot_surface.png

Displaying nodal values

This feature is particularly useful if, for example, one needs to display the result of a finite element computation where the data is known at the vertices. In the following example, we create a planar mesh with triangular elements. Then we define a linearly increasing nodal data along the x-axis. The full code is available at demo/demo_graphics/nodal_values.cpp.

// geometric data
matrix<> X,Y;
std::tie(X,Y) = meshgrid(linspace(-1,1,10),linspace(-1,1,5));
X = vertcat(X,X);
Y = vertcat(Y,Y);
matrix<> Z = zeros(size(X));

// create mesh
matrix<> vtx;
matrix<std::size_t> elt;
std::tie(elt,vtx) = tridelaunay(X,Y,Z);

// display
figure fig;
trimesh(fig,elt,vtx,eval(vtx(row(vtx),0)));
drawnow(fig);

The result should look like this:

_images/nodal_values.png

From mesh generation to file output

In this example, we create a spherical mesh and compute the normals to the triangles. We plot both on the same figure. Finally, we save the mesh to a .ply file. We also illustrate the use of the quiver function. The full code is available at demo/demo_graphics/advanced_mesh.cpp.

First, we create the mesh using the sphere2 function which creates a Fibonacci sphere.

std::size_t nvtx=100;
// 1. Create the mesh
matrix<> vtx;
matrix<std::size_t> elt;

matrix<> X,Y,Z;
std::tie(X,Y,Z) = sphere2(nvtx); // Fibonacci sphere
// Delaunay tetrahedrisation
std::tie(elt,vtx) = tetdelaunay(X,Y,Z);
// Boundary extraction
std::tie(elt,vtx) = tetboundary(elt,vtx);

Then, we compute the center ctr of the triangles and the normal vector nrm. We recall that the coordinates of ctr are simply the averaged sum of the coordinates of the vertices of the triangles. The coordinates of nrm may be determined by computing the cross-product between the tangent of, two consecutive edges. In this example we normalize the length to 0.25 to get an equilibrated display.

// 2. Compute the normal vectors and the centers of the triangles
std::size_t nelt = size(elt,1);
matrix<> nrm = zeros(nelt,3);
matrix<> ctr = zeros(nelt,3);
for(std::size_t ie=0; ie<nelt; ++ie)
{
    // center
    for(std::size_t i=0; i<3; ++i)
    {
        ctr(ie,i) = (vtx(elt(ie,0),i)+vtx(elt(ie,1),i)+vtx(elt(ie,2),i))/3.;
    }
    // normal vector to triangle {A,B,C}
    matrix<> AB=zeros(1,3), BC=zeros(1,3), nv = zeros(1,3);
    for(std::size_t i=0; i<3; ++i)
    {
        // tangent to first edge
        AB(i) = vtx(elt(ie,1),i) - vtx(elt(ie,0),i);
        tangent to second edge
        BC(i) = vtx(elt(ie,2),i) - vtx(elt(ie,1),i);
    }
    // for single vectors, this code is faster than
    // a call to 'cross' (for the cross-product) or 'norm'.
    nv(0) = AB(1)*BC(2) - AB(2)*BC(1);
    nv(1) = AB(2)*BC(0) - AB(0)*BC(2);
    nv(2) = AB(0)*BC(1) - AB(1)*BC(0);
    double l = std::sqrt(nv(0)*nv(0)+nv(1)*nv(1)+nv(2)*nv(2));
    // normalization
    for(std::size_t i=0; i<3; ++i)
    {
        nrm(ie,i) = nv(i)/(2*l); // arrows of length 0.5
    }
}

Now, we plot the result.

// 3. Plot everything
figure fig;
trimesh(fig,elt,vtx); // plot the mesh
quiver(fig,ctr,nrm);  // plot the normal vectors at the centers
drawnow(fig);

Finally, save the mesh into the current directory in the .ply format.

// 4. save to .ply
std::string path="./", name="testfile.ply";
triwrite(path,name,elt,vtx);

The result should look like an urchin, see below. Please note that the normal vectors may not appear when the window appears. In that case, simply clicking inside the window should do the trick.

_images/advanced_mesh.png