Getting FreeFem++ to display impressive visuals under Linux.

One of my Computing habits is, to acquire many frameworks, for performing Scientific or Analytical Computing, even though, in all honesty, I have little practical use for them, most of the time. They are usually of some Academic curiosity to me.

Some of the examples familiar to me are, ‘wxMaxima‘ (which can also be installed under Linux, directly from the package manager), ‘Euler Math Toolbox‘ (which, under Linux, is best run using Wine), and ‘SageMath‘ (which IMHO, is best installed under Linux, as a lengthy collection of packages, from the standard repositories, using the package manager, that include certain ‘Jupyter’ packages). In addition to that, I’d say that ‘Python‘ can also be a bit of a numerical toolbox, beyond what most programming languages can be, such as C++, yet, a programming language primarily, which under Linux, is best installed as a lengthy collection of packages through the package manager. And a single important reason is the fact that a Python script can perform arbitrary-precision integer arithmetic natively, and, with a small package named ‘python3-gmpy2′, can also perform arbitrary-precision floating-point arithmetic easily. If a Linux user wanted to do the same, using C, he or she would need to learn the library ‘GMP’ first, and that’s not an easy library to use. Also, there exists IPython, although I don’t know how to use that well. AFAICT, this consists mainly of an alternative shell, for interacting with Python, which makes it available through the Web-interface called “Jupyter”. Under Debian Linux, it is best installed as the packages ‘ipython3′, ‘python3-ipython-genutils’, ‘python3-ipykernel’, ‘python3-nbconvert’, and ‘python3-notebook’, although simply installing those packages, does not provide a truly complete installation… Just as one would want a ‘regular’ Python installation to have many additional packages, one would want ‘IPython’ to benefit from many additional packages as well.

But then, that previous paragraph also touches on an important issue. Each Scientific Computing platform I learn, represents yet-another scripting language I’d need to learn, and if I had to learn 50 scripting languages, ultimately, my brain capacity would become diluted, so that I’d master none of them. So, too much of a good thing can actually become a bad thing.

As a counter-balance to that, it can attract me to a given Scientific Computing platform, if it can be made to output good graphics. And, another Math platform which can, is called “FreeFem“. What is it? It’s a platform for solving Partial Differential Equations. Those equations need to be distinguished from simple derivatives, in that they are generally equations, in which a derivative of a variable is being stated on one side (the “left, bilinear side”), but in which a non-derivative function of the same variable is being stated on the other (the “right side”). What this does, is to make the equation a kind of recursive problem, the complexity of which really exceeds that of simple integrals. (:2)  Partial Differential Equations, or ‘PDE’s, are to multi-variable Calculus, as Ordinary Differential Equations, or ‘ODE’s, are to single-variable Calculus. Their being “partial” derives from their derivatives being “partial derivatives”.

In truth, Calculus at any level should first be studied at a University, before computers should be used as a simplified way of solving its equations.

FreeFem is a computing package, that solves PDEs using the “Finite Element Method”. This is a fancy way of saying, that the software foregoes finding an exact analytical solution-set, instead providing an approximation, in return for which, it will guarantee some sort of solution, in situations, where an exact, analytical solution-set could not even be found. There are several good applications. (:1)

But I just found myself following a false idea tonight. In search of getting FreeFem to output its results graphically, instead of just running in text mode, I next wasted a lot of my time, custom-compiling FreeFem, with linkage to my many libraries. In truth, such custom-compilation is only useful under Linux, if the results are also going to be installed to the root file-system, where the libraries of the custom-compile are also going to be linked to at run-time. Otherwise, a lot of similar custom-compiled software simply won’t run.

What needs to be understood about FreeFem++ – the executable and not the libraries – is, that it’s a compiler. It’s not an application with a GUI, from which features could be explored and evoked. And this means that a script, which FreeFem can execute, is written much like a C++ program, except that it has no ‘main()‘ function, and isn’t entirely procedural in its semantics.

And, all that a FreeFem++ script needs, to produce a good 2D plot, is the use of the ‘plot()‘ function! The example below shows what I mean:

 

Screenshot_20210708_232521

 

I was able to use an IDE, which I’d normally use to write my C++ programs, and which is named “Geany”, to produce this – admittedly, plagiarized – visual. The only thing I needed to change in my GUI was, the command that should be used, to execute the program, without compiling it first. I simply changed that command to ‘FreeFem++ "./%f"‘.

Of course, if the reader wants in-depth documentation on how to use this – additional – scripting language, then This would be a good link to find that at, provided by the developers of FreeFem themselves. Such in-depth information will be needed, before FreeFem will solve any PDEs which may come up within the course of the reader’s life.

But, what is not really needed would be, to compile FreeFem with support for many back-ends, or to display itself as a GUI-based application. In fact, the standard Debian version was compiled by its package maintainers, to have as few dependencies as possible (‘X11′), and thus, only to offer a minimal back-end.

(Updated 7/14/2021, 21h45… )

(As of 7/09/2021, 0h10: )

 

Screenshot_20210709_003013

 

Surprisingly, I discovered that even the bare-bones Debian 9 / Stretch version of FreeFem++ has ‘ffmedit’ and ‘VTK 2′ (file-writer) support. However, the way the ffmedit viewer displays a mesh is disappointing, because it doesn’t colour-code the mesh. This makes the display unworthy to be shown in public. (:3) However, how well I can view VTK Files, really only depends, on how much I want to play with the settings, of the stock version of ‘Paraview’…

 

Screenshot_20210709_011313

 


 

(Update 7/10/2021, 20h40: )

1:)

An approach that can be used for polynomials, which will find a scalar of roots sequentially, each time zeroing in on a root that ‘works by itself’, cannot be used for PDEs. Because PDEs state that a Partial Derivative, a function of perhaps 3 variables, needs to equal another value, also computed from the same variables, it follows that the solution of such an example is inherently multi-dimensional, perhaps consisting of a membrane occupying 3 dimensions, with an inequality instead, for all parts of a volume not exactly on that membrane. Or, with similar probability, the solution could exist as one unique combination of the 3 variables, none of which can be found, without also finding the other 2.

I have to admit, however, that the way PDEs are “solved” in practice, differs, from how I had imagined they would be solved. My previous idea was, that each variable would be subdivided into increments within its domain, but that a fixed function would be fed values of these variables, to determine whether those values satisfy the PDE.

What happens instead is, that an approximation of a function is computed – essentially as an array of output-values – in such a way as to get as close as possible to satisfying the PDE. This array is not an analytically defined function, but then acts as a replacement for one. (:4)


 

The Debian version of the packages, also offers a solver named ‘ff3d‘, which accepts a script as well as a POV-Ray File as input. AFAIK, that solver simply takes the defined mesh of the POV-Ray File as its domain – which is referred to as “a fictitious domain”.



 

2:)

A question which the uninitiated may not know the answer to could be, ‘How do differential equations imply integral equations, or vice-versa?’ And the answer which I would offer is, that if a person has been given ‘an integral equation’ to solve, that person can just differentiate both sides of it. What was the integral side of the equation will become non-integral, and what was the non-integral side, will become the derivative side.

It’s just that, because derivatives are generally slightly easier to solve, than integrals, transforming integral equations into differential equations in this way, may ease solving either form.

Also, an integral equation which states fluid flow, can also be differentiated.

 

Another piece of insight which I can offer, about the subject of ‘Only differentiating an implicit function on both sides of the equation, and then solving’, is that in my own experience, when one does this, the solution one will find, is the derivative of the original, implicit function, not, the solution of the original, implicit function.

The reason for which the PDE is formally written as an equality-or-inequality between integrals, of derivatives, would seem to be the fact that in order for the function’s derivative to satisfy a constraint, the function itself must be the (negative) integral of that constraint.


 

3:)

Above, I had written, that the way ‘ffmedit‘ works, where that is supposed to be the preferred way to view ‘FreeFem++’ meshes, was disappointing, because it was not colour-coding the meshes.

Further, something which I had written about ‘ff3d‘ left the question unanswered, whether that tool, when outputting data for a POV-Ray Mesh, would be able to generate a visual.

To my great satisfaction, I’ve discovered that both goals can, after all, be achieved! :D

The actual executable ‘ff3d‘ (under Debian 9) has no dependency on ‘X11′, which means, that it will not display anything directly. However, it can be scripted to output 1 or more .MESH Files, which in turn can be viewed with ‘ffmedit‘.

When doing this, it may be helpful to type the name of the mesh to be viewed on the ‘ffmedit’ command-line, without the ‘.mesh’ file-name extension, just so that the viewer will load both the ‘.mesh’ and the ‘.bb’ Files with the same base-name.

This viewer, in turn, has a context-menu which I was not aware of before, which one obtains by right-clicking on the displayed mesh, and that menu notably has a ‘Data’ sub-menu, which allows colour-coding to be set as desired. So, it should be possible to obtain everything from these two executables, that one might desire, from 3D PDEs specifically, which are being computed along a POV-Ray Mesh:

 

Screenshot_20210709_141658

 

Screenshot_20210711_214444

 

Tada!


 

4:)

I have just been studying the FreeFem documentation, to dig deeper into how it can be used. And one fact which I learned was, that the Finite Element Method I had believed in earlier, falls short, of what the actual software package can do. Looking at the tutorial, what I found was, that the following code:

 


fespace Vh(Th, P1);
Vh u, v;// Define u and v as piecewise-P1 continuous functions

 

Does not truly “define” the functions (u) and (v), but rather, ‘declares’ them to be (continuous) functions, which are to receive one value for each of  the elements belonging to ‘Th’. What the scripting language offers next, goes further than just, to compute what the derivative of defined functions (u) and (v) would be. The following code:

 


solve Poisson(u, v, solver=LU)
    = int2d(Th)(    // The bilinear part
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(Th)(    // The right hand side
          f*v
    )
    + on(C, u=0);   // The Dirichlet boundary condition

 

Actually populates the elements of (u), with values that satisfy the condition, that the expression to be solved for will be close to zero in value. This same block of code does not populate (v) with elements. But it is only due to this second block of code ‘solving’, that (u) can actually be plotted.


 

(Update 7/11/2021, 15h55: )

There is something specific, which FreeFem does, which can make it harder for people like me, who have zero experience with FreeFem, to learn its basic usage. I tend to be a person, who needs to have some idea of what, approximately, a software package is doing procedurally, before I’d be able to set up my own problems for it. And in the case of public FreeFem documentation, there is no mention of that.

As can be seen in the code block above, the user is allowed – and required – to make use of the terms ‘dx(v)‘ as well as ‘dy(v)‘, even though (v) has not been initialized. And, after the solver exits, (v) still has no values. (v) is simply mapped over the same domain of (x) and (y), which (u) is mapped over, except that (u) ends up receiving the values, that will form the solution.

What, exactly, the solver does, could be primitive or complex. If the solver is primitive, what it can do is, to keep cycling through the elements of (u) and (v) concurrently, obtain non-zero values from the expression to be solved for, and subtract those values from the elements of (u), until the final pass over the domain yields results from the expression that are desirably close to zero…

At the very least, this would require that the solver be able to obtain a numerical result from the expression to be solved for.

It seems to me, that when the solver encounters the term ‘dx(v)‘, for an element of (v) that has no parameter, it simply evaluates that to (1.), without setting the parameter. Similarly, when the solver encounters ‘f*v‘, it simply treats (v) as being equal to (1.).

However, attempting to plot (v) afterwards simply generates an error message, over the same issue.

The user can do many things with FreeFem. What the user may not do is, to leave out the (v) term, or alternatively, the ‘dx(v)‘ and ‘dy(v)‘ terms, from the integrals to be solved for. And the reason seems to be the fact that, in a fixed way, the solver looks for those terms, to define the interval over which the integrals are to be computed (numerically).


 

(Update 7/11/2021, 22h50: )

What I have also determined by experimentation is that, if the user initializes (v) in the example above, before calling ‘solve ...‘, what FreeFem will do is, just ignore whatever those elements were initialized to…

Hence, the following code is fully legal, although it gives a completely different result, from what the Poisson code gave:

 


// Load Cubic polynomial interpolator
load "Element_P3";

// Caution:
// Cubic interpolations are often unneccessary, and may
// overshoot their endpoints.
// They're usually only needed, for second-degree
// Derivatives, which lead to second-order Differential Equations.

// A P3-derived fespace will nevertheless be used here.

// For label definition.
int Dirichlet=1;

// Define mesh boundary.
border C1(t=0, 2*pi){x=cos(t); y=sin(t); label=Dirichlet;}
border C2(t=0, 2*pi){x=cos(t) * 0.2 + 0.3; y=sin(t) * 0.2 + 0.3;
						label=Dirichlet;}

// The triangulated domain Th is on the left side of its boundary.
mesh Th = buildmesh(C1(100) + C2(-40));
plot(Th, ps="ThWithHole.eps", wait = true);

// Define a function f.
func f = x * y;

// The finite element space defined over Th is called Vh here.
fespace Vh(Th, P3);
Vh u, v = f;	// Declare u and define v as piecewise P3-continuous functions.

// Define a Vh-fespace function equivalent to f.
Vh fh = f;

// Get the clock in seconds.
real cpu=clock();

// Define a simple PDE...
solve Poisson(u, v, solver=LU)
	= int2d(Th)(    // The bilinear part
		  (dx(u) * dx(v))
		+ (dy(u) * dy(v))
	)
	- int2d(Th)(
		v
	)
	+ on(Dirichlet, u=0);   // The Dirichlet boundary condition

// Plot the result...
plot(u, ps="PoissonForUnity.eps", wait = true);
plot(v, ps="VDoesNothing.eps", wait = true);


u[] = 0.;

// Define a twisted PDE...
solve TwistedPDE(u, v, solver=LU)
	= int2d(Th)(    // The bilinear part
		  (dx(fh) * dx(u) * dx(v))
		+ (dy(fh) * dy(u) * dy(v))
	)
	- int2d(Th)(
		f*v
	)
	+ on(Dirichlet, u=0);   // The Dirichlet boundary condition

// What I had hoped, the first plot would show...
plot(u, ps="TwistedPDE.eps", fill=true, nbiso=64);

// Display the total computational time.
cout << "CPU time = " << (clock()-cpu) << endl;

 

Also, just to practice, my hypothetical code put a hole in the mesh, which it is solving for. It generates the following two plots…

 

Screenshot_20210711_173216

 

Screenshot_20210711_182645

 

So far, my code has only achieved that the second of these two plots, which is of (v), no longer generates an error message. Yet, inspection of the first plot reveals, that it will behave as though ‘v‘, which was supposed to be equivalent tof‘, had been defined as (1.). Hence, to have initialized (v), had no effect whatsoever, on the first plot.


 

Now, the following plot shows what I would have expected the first plot to show, which, translated into English, would mean, ‘Find the fespace function (u) such that its gradients with respect to (x) and (y), multiplied with the corresponding gradients of (f), and the results added, will yield the values of (f).’ …

 

TwistedPDE

 

And there I have that twisted PDE, which has no practical value, but solved non-trivially.


 

(Update 7/12/2021, 8h20: )

I suppose that I have yet another observation to add. The way the Finite-Element Space function (u) is being colour-coded, seems somewhat lacklustre, in the bottom-left and top-right quadrants of the plot…

I find it to be revealing, that the plot of (u) is in fact negative, where (f) is negative. The gradient of (f) itself is positive – at least, going from left to right, along the top half of the plot… But, the gradient of (f) is negative, going from bottom to top, along the left half of the plot. Along the bottom half of the plot, the gradient of (f) should then be negative again, going from left to right… I guess that the PDE’s plot has proved, that its Y-axis is facing upward, and not, facing downward.

What this results in is the fact, that the X- and Y-gradients of (f) oppose each other in the top-left and bottom-right quadrants, while agreeing with each other, in the top-right and bottom-left quadrants. And this, in turn, seems to require less amplitude from (u), than is required in the top-left and bottom-right quadrants.


 

(Update 7/12/2021, 18h00: )

There is an admission which I must make about this posting, and what it fails to explain about the integration which FreeFem can perform, on a disk. In Calculus 2, what Students are generally taught is twofold:

  • There exist definite and indefinite integrals, where, the definite variety is derived from the indefinite variety, by computing the indefinite integral at a point in its domain, where the definite integral is predetermined to equal zero, and then to subtract this value from all the other indefinite integrals, to arrive at the definite integrals, And
  • There exists integration over 2 dimensions, in which each dimension is categorized as an orthogonal interval, defined by one coordinate-variable.

The problem which this posting requests that FreeFem solve, poses two issues with respect to those teachings:

  1. The mesh is a disk-shaped mesh, not a rectangle, And
  2. The value of (u), in this case, is not just supposed to equal zero at one point, but is actually supposed to equal zero, everywhere along the boundaries of the disk, in addition to satisfying its other constraints.

In truth, I did not study Calculus beyond Calculus 2, for which reason I can only give a partial answer, to how these two issues are probably resolved.

Issue 1, which is also referred to as “The Dirichlet Problem“, has as added challenge, that the vertices of the mesh generated by this script are poorly ordered; however, each vertex has an (x) and a (y) coordinate. And I think that it gets resolved, in the way the mesh was built, which was, as a set of concentric rings, the angle of which is supposedly orthogonal to the radius. Hence, the problem can be broken down into a classical, 2-variable problem, due to how the mesh was created. The natural order of the elements, which only have one official index in this case, completes one ring after the other. Each ring passively receives a different weight, according to the circumference it represents.

This issue would have gotten sidestepped completely, had I instructed FreeFem to generate a rectangular mesh instead.

The reason for which FreeFem is able to compute ‘dx(u)‘ and ‘dy(u)‘ correctly, which translate into ‘the derivative of (u) with respect to (x)’ and ‘the derivative of (u) with respect to (y)’, is the fact that a mesh has in fact been built, so that vertices additionally belong to triangles.

The following plot displays the mesh, which FreeFem created for me:

 

ThWithHole

 


 

(Update 7/12/2021, 18h30: )

According to my latest ruminations, the terms ‘dx(v)‘ and ‘dy(v)‘ should really represent numerically, By how much the x- and y-coordinate have changed, from the previous vertices’.

What’s confusing about that is, the fact that it’s inconsistent, with the way ‘dx(u)‘ and ‘dy(u)‘ are parsed.


 

(Content removed 7/13/2021, 1h45 because inaccurate. )


 

(Update 7/14/2021, 14h30: )

There is a hypothetical computation which I can imagine FreeFem to be performing, which would act as a surrogate for the 2-dimensional integrals which I learned in Calculus 2, and which could act as a point of reference, when using the program. Given that the disk’s vertices only have 1 official index, but, that the vertices all belong to a triangular mesh, if the definite integral of 2 points of each triangle is already known, a sort of definite integral of the 3rd point can also be computed.

But, in order for such an approach really to work, what would first need to be done is, that the ‘dx(v)‘ and ‘dy(v)‘ of each edge must form the basis of a linear equation, so that two resulting linear combinations will yield:

 


ax * edge1 + bx * edge2 ->
dx(v) == 1., dy(v) == 0.

ay * edge1 + by * edge2 ->
dx(v) == 0., dy(v) == 1.

 

In other words, the solution would generate error messages, as soon as any two edges are parallel, or, as soon as the length of any edge is zero. (:5)

 

In order to arrive at ‘f*v‘, the function that can be computed at the unsolved vertex (x0,y0) would be:

 


2DCross( (s1,t1), (s2,t2) ) :=
    (s2 * t1) - (s1 * t2)

a_sign = sign( sign(dx(v)) + (0.5 * sign(dy(v))) )

AreaQuad( (x0,y0), (x1,y1), (x2,y2) ) :=
    |2DCross( ((x0-x1),(y0-y1)), ((x0-x2),(y0-y2)) )|

u0 = (u1 + u2) / 2 +
    ( a_sign *
    f(x0,y0) * AreaQuad( (x0,y0), (x1,y1), (x2,y2) ) )

 

(Revised 7/14/2021, 19h10. )

The result could simply be added to the mean of the already-solved vertices’ solutions, and the sum offered as the solution for the 3rd vertex.

 

In order to arrive at either partial gradient by itself, only 1 out of the 2 available linear combinations would get used.

 


dx(v) = x0 - ((x1 + x2) / 2)
dy(v) = y0 - ((y1 + y2) / 2)

( ax * (f(x0,y0) - f(x1,y1)) + bx * (f(x0,y0) - f(x2,y2)) ) * dx(v)  OR
( ay * (f(x0,y0) - f(x1,y1)) + by * (f(x0,y0) - f(x2,y2)) ) * dy(v)

 

This ‘sounds nice’, especially since the boundary value of the disk’s outer boundary is set, and constitutes 2 vertices of a triangle each time (see plot above). A ring of vertices could then be derived, ‘1 step inside’ each ring of vertices that has already been solved.

Such a kludge would additionally remain consistent with the way in which I was taught integrals, in that the values computed in this way for ‘f*v‘ would end up having amplitudes, proportional both to those of ‘f(x,y)‘ and to the area of the plot.


 

(Update 7/13/2021, 18h20: )

5:)

Breaking down this thought process into baby-steps, the methodology of computing the multipliers (ax), (bx), (ay) and (by), which I mentioned above, involves a simple application of Linear Algebra, but is included in the following work-sheet, for readers who do not wish to perform the work, to see some sort of results:

http://dirkmittler.homeip.net/FreeFem_1.html

On most browsers, the reader will need to enable JavaScript from my site, as well as from ‘mathjax.org’, to be able to view the sheet.


 

(Update 7/14/2021, 21h45: )

FreeFem cheats, at least, when using a disk-shaped 2D plot !

There was a basic question going around my head, about how the developers of FreeFem solved a problem, that would be inherent, in outputting a numerically computed integral to a 2D Mesh, in the form of a disk with unordered vertices, that have a single index number. And the answer which I have found was, that they did not, in fact, solve the following problem…

Let’s say that the function to be integrated, f(x), is just (x). This is the simplest linear function, I’d say. The indefinite integral is:

0.5 * x2 + C

Where (C) is the arbitrary constant.

It’s important to understand the significance of this arbitrary constant, that also characterizes all indefinite integrals as distinct from definite integrals. These integrals would be a parabola, concave upwards.

If it was our intent to make the endpoints of this parabola equal a specific boundary value, such as zero, the classical way to do that would be, to set the arbitrary constant to (-0.5). The following plot shows what I mean:

 

Parabola_C_1

 

‘The problem’ to be considered, if this function needed to be computed as a discrete summation from both ends, is, that the summation from the right-hand side of the plot, where (x= +1), towards the origin, would need to be done with inverted sign, because the integral is becoming increasingly negative, even though the function being summed, is positive. OTOH, the summation coming from the left-hand-side, where (x = -1), also needs to lead to symmetrically negative values, because the function being summed was itself negative.

What this means is that, where a summation from left to right took place positively, the same summation going from right to left needs to be negated, in order to preserve the sense of a true integral.

I’ve known this from the start of this posting. And, one of the questions I’ve been breaking my head over was, ‘Where does FreeFem derive the negation?’

The answer is, ‘It does NOT!’ Given a disk to be summed from the outside in, FreeFem applies each summation with the same, positive sign. The following plot shows the result:

 

UnityIntegrated

 

What’s more, if an attempt were made, as I found out, to introduce negations to the summations in any systematic way, doing so would lead to artifacts in the way the plot ‘looks’. It would just never seem to come out right.

For this reason, when plotting a disk, even though by naming, FreeFem appears to be doing something ~in 2D~, it’s really just a 1-dimensional integral, radially inwards.

For that reason, its amplitudes will actually reach those of the original function squared, but never, those of the original function cubed.

 

Dirk

 

Print Friendly, PDF & Email

2 thoughts on “Getting FreeFem++ to display impressive visuals under Linux.”

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>