One of the programs which I’ve been experimenting with, is called “FreeFem”. Actually, it exists both as a library, as well as a set of executables, the latter of which are meant to facilitate its use, in that scripts can be written in a loose but C++-like syntax, that will already produce interesting plots. But, plots can also be meaningless, unless the user understands how they were generated, as well as the underlying Math.

In an earlier posting, I already wrote that a basic weakness of the ~standard~ 2D integrals which FreeFem will solve for, was, that those do not correspond closely to “Double Integrals”, as I was taught those in Calculus 2.

Just to recap basic Calculus 2: If an integral is supposed to be plotted over two dimensions, then the standard way to do so is first, to declare a variable that defines one of the dimensions, and to declare a second variable that defines the second. Next, a 1-dimensional integral is defined over the inner variable, which results in a continuous function. Next, this continuous function is integrated a second time, over the outer variable, thereby resulting in a double integral.

Even though this could be extended to more dimensions than 2 or 3, there is usually little practical value in actually doing so, at least, in my limited life.

Further, the fact is clear to me that integrals exist, which go beyond this basic definition. For example, “Curl Integrals” also exist… However, for the moment I’m focusing on how to overcome some of the limitations which are imposed by FreeFem. And while working on this problem, I have found a way to force FreeFem to compute the sort of double integral which I was taught. The script below shows how I did this…

```
// For label definition.
int Bottom=1;
int Right=2;
int Top=3;
int Left=4;
// The triangulated domain Th is on the left side of its boundary.
mesh Th = square(10, 10, [(x*2)-1., (y*2)-1.]);
plot(Th, ps="ThRectX.eps", wait = true);
// Define a function f.
func f = x;
// The finite element space defined over Th is called Vh here.
fespace Vh(Th, P2);
Vh u, v; // Declare u and v as piecewise P2-continuous functions.
// Get the clock in seconds.
real cpu=clock();
// Define a simple PDE...
solve SimpleXInt(u, v, solver=LU)
= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) // Show me a kludge.
- int2d(Th)(
f*v
)
+ on(Bottom, u=0) // The Dirichlet boundary condition
+ on(Right, u=0)
+ on(Top, u=0)
+ on(Left, u=0);
// Plot the result...
plot(u, ps="RectX.eps", value=true, wait = true);
// Now, let's try to plot a double integral...
Vh A, Fh;
func real f2(real xx, real yy) {
// return xx * yy;
return xx;
}
func real Fx() {
real s = 0.;
for (int i = -5; (i * 0.2) < x; ++i) {
s += int1d(Th, Left) ( f2(i * 0.2, y) );
}
return s;
}
// Compute the x-integrals for reference purposes...
Fh = Fx();
solve Double(A, v)
= int2d(Th) (
dy(A) * dy(v)
)
- int2d(Th) (
Fx() * v
)
+ on(Bottom, A=0) // The Dirichlet boundary condition
+ on(Right, A=0)
+ on(Top, A=0)
+ on(Left, A=0);
plot(A, ps="DoubleIntX.eps", value=true, wait=true);
plot(Fh, value=true);
// Display the total computational time.
cout << "CPU time = " << (clock()-cpu) << endl;
```

This script actually outputs 4 plots and saves the first 3. But, for the sake of this posting, I’m going to assume that there is no reason to show the reader a plot, of a plain, rectangular mesh… The following plot shows the result, when ‘the Poisson Equation’ is simply given for a function of (x), such that:

f(x) = x

A simple linear function was given. It can be integrated once or twice, and the plot above shows how FreeFem usually does so, resulting in u(x,y).

One interesting fact about FreeFem is, that this program can be made to compute a 1-dimensional integral explicitly, i.e., without imposing any need, to solve an equation. The resulting sum (a real number) can then simply be returned with a variable, for every maximum value the function was summed to. So, what would be tempting to do next is, to integrate this variable a second time, which might result in a 2-dimensional array. However, there is a basic limitation in how FreeFem works, that would next stand in our way… (:1)

And so, in order for the second (outer) integral actually to be with respect to (y), I used the little trick in the script above. I defined the gradient of Finite Element function (A), with respect to (y) only, and stipulated that it must equal the 1-dimensional integral which was computed before, for all combinations of (x,y). Hence, this second step of my solution is similar to a degenerate Poissson equation, in that the addition of the (x) gradient has simply been dropped.

When this approach gets applied, the following plot results:

Now, there is still something which FreeFem does, which is supposed to be a feature, ~~and which I cannot switch off~~. That is, to impose that at the boundaries, a certain exact value needs to exist for (A), just as was needed for (u). Documentation states, that by merely defining an equality, the (2D) derivative of a function is supposed to have with some other function, one has not fully defined that function. A boundary value must also be given, and then the resulting “Partial Differential Equation” (‘PDE’) actually defines, either the FE function (u), or the FE function (A). (:2) Further, the way ‘`int1d()`

‘ is implemented, also applies such boundaries.

For that reason, I can also ~~not~~ just shut off, anything which I feel FreeFem might be doing wrong, every time it tries to solve a PDE. What I can do is, request that this boundary occupy the left-hand side of the rectangular plot, and additionally, make sure that I begin my summation with the same value, thus resulting in a definite and not an indefinite integral every time. Also, I had suspected somewhere, that in order to resolve this issue with 2D plots, it really only needs to be resolved along one inconspicuous axis. In this case, I resolved it explicitly for the X-axis, and the behaviour of the Y-axis fell in line.

Yet, what I can do is observe, that the two plots shown above don’t match. And as long as they don’t, what FreeFem computed in its first plot, was also not a double integral.

(Updated 7/25/2021, 17h30… )

**1:)**

There is a largely undocumented implementation detail to FreeFem, which slowed down my attempt to produce this double integral considerably. Apparently, the function ‘`int1d()`

‘ has as behaviour, just to evaluate the function it has been fed at the current values of (x,y), ~~but, to multiply that result by whatever the ~~. The reason the developers did this, was probably, the fact that doing so fits in well, with how FreeFem solves its 2D PDEs. But, for a user who wants to generate 1-dimensional integrals, what this does is, to force that user to do the work himself, of looping through the domain of x-values, that lead to the current x-value, and to perform the summation, ~of *width* of one element of the current mesh is*outputs of* ‘`int1d()`

‘~, to arrive at that 1-D integral, and to do what any Scientific Computing platform should automate.

Furthermore, this function has as evil, that if the convention is used, which works elsewhere within FreeFem, just to pass in a function reference, no function call will actually take place! The name of the function must be written, with a set of parentheses, in order actually to cause a function-call. (Which is also, how it usually works with pure programming languages. It’s just another inconsistency, that FreeFem will sometimes allow the user to pass in the name of the function by itself, causing function-calls.)

Eventually, one discovers these little details.

(Update 7/15/2021, 22h05: )

Astute readers will notice the fact that, even though I’ve complained about FreeFem developers ‘cheating’, * I also cheated*. Code which truly forces a double integral to be computed for a 10×10 mesh, requires that about 1,000,000 vertices be computed in total, and takes about 400 seconds of CPU time on my machine, to complete. Further, the results ‘look wrong’, because when:

```
f2(xx, yy) := xx ->
dy(f2) == 0.
```

This means that, if the ‘outer integration’ has been computed accurately, *vertical lines that correspond to individual values of (x)* will slope in a completely uniform way along the Y-axis. And the following is the *evil* code that does this:

```
// For label definition.
int Bottom=1;
int Right=2;
int Top=3;
int Left=4;
// The triangulated domain Th is on the left side of its boundary.
mesh Th = square(10, 10, [(x*2)-1., (y*2)-1.]);
plot(Th, ps="ThRectX.eps", wait = true);
// Define a function f.
func f = x;
// The finite element space defined over Th is called Vh here.
fespace Vh(Th, P2);
Vh u, v; // Declare u and v as piecewise P2-continuous functions.
// Get the clock in seconds.
real cpu=clock();
// Define a simple PDE...
solve SimpleXInt(u, v, solver=LU)
= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) // Show me a kludge.
- int2d(Th)(
f*v
)
+ on(Bottom, u=0) // The Dirichlet boundary condition
+ on(Right, u=0)
+ on(Top, u=0)
+ on(Left, u=0);
// Plot the result...
plot(u, ps="RectX.eps", value=true, wait = true);
// Now, let's try to plot a double integral...
Vh A, Fh;
func real f2(real xx, real yy) {
// return xx * yy;
return xx;
}
func real Fx(real xx, real yy) {
real s = 0.;
for (int i = -5; (i * 0.2) < xx; ++i) {
s += int1d(Th, Left) ( f2(i * 0.2, yy) );
}
return s;
}
// Compute the x-integrals for reference purposes...
Fh = Fx(x, y);
func real Fy(real xx, real yy) {
real s = 0.;
for (int i = -5; (i * 0.2) < yy; ++i) {
s += int1d(Th, Bottom) ( Fx(xx, i * 0.2) );
}
return s;
}
// This is where the double integral will be computed,
// with ugly results...
A = Fy(x, y);
// This takes about 400 seconds of CPU time on my machine.
plot(A, ps="DoubleIntX.eps", fill=true, value=true, wait=true);
// Notice how a plot, that has straight edges
// along the y-axis, is 'true', but also 'looks wrong'.
// Those straight edges really follow from the fact,
// that f2(x,y) = (x) has a slope of 0. with respect to (y).
plot(Fh, value=true);
// Display the total computational time.
cout << "CPU time = " << (clock()-cpu) << endl;
```

The resulting Y-axis slope will be strongest, where the X-axis *integral* was already, conspicuously negative and asymmetrical. (:3)

(Update 7/16/2021, 5h40: )

I suppose that a relevant question which readers might have, could be, whether it should always take 1,000,000 vertex-computations, to compute the true, double integral, on a 10×10 mesh. And the answer is, ‘Absolutely Not!’

One strategy which can and will be used, is, to save temporary summations, such as of one line while computing one x-axis integral, or of one array after having computed a full set of those, and before computing the y-axis integrals, so that the summation can be continued, from saved temporary summations. It’s a perfectly legitimate thing to do.

However, a basic limitation which FreeFem has is, that it stores its meshes in a way that’s easy to write to, but next to impossible to read information from. For example, *even in their rectangular meshes*, the vertices have *only one index*, while, to perform integrals, what’s really called for, is an x-index and a y-index. And, it’s because one cannot easily use these meshes as input, that I decided to solve the problem, by writing nested loops, that perform a ridiculously excessive number of summations.

Each of my vertices in this example, was being visited up to 10,000 times.

(Update 7/17/2021, 13h10: )

* Conclusion*:

After this amount of experimentation, I have decided that I will trust the output of FreeFem, with two caveats:

- The fact needs to be understood that, even though Integrals over 2 variables were always taught as second-order integrals in Calculus 2 – i.e., double integrals – FreeFem computes them as a sort of first-order integral each time, putting a first-order integral form on both the right-hand and the left-hand side of the equation to be solved for. By forcing FreeFem actually to have a second-order integral on the right-hand side, I’m creating an equation which, Mathematically, is less symmetrical than what I had before, but which FreeFem will be able to solve for anyway, And
- In order to solve the PDE, FreeFem ‘feels free’ to flip the signs of the bilinear terms, as is convenient. This may in fact be a necessity inherent in the problem, that FreeFem solves. And, this peculiarity may result passively, from the notion that the ‘
`int2d()`

‘ function itself ignores ‘`dx(v)`

‘ and ‘`dy(v)`

‘, as it’s computing its summation on both sides of the equation – including, in the summation of the bilinear terms. (:4)

Also, I did have some concerns, over the possibility that the ‘`int2d()`

‘ function computes more of an ad-hoc summation, than an integral. But, what I’m noticing is that FreeFem developers put enough attention into their ad-hoc summation, that it will truly pass for an integral, because of one consideration:

- The same type of ad-hoc summation is being performed on both sides of the equation, matching for each vertex.

(Update 7/17/2021, 14h25: )

**2:)**

Even though the answer to this question has already been stated elsewhere, using tough Mathematical expressions, the answer could be sought in layman’s terms, or in ‘horse sense’, to ‘Why is the Dirichlet problem valid, in that the PDE is under-defined, without a boundary condition being stated?’

The first basic notion needed, to answer that question, is, that if we are given an equation to solve, in which both (x) and (y) are the variables, but only given that:

x + y = A

Where, (A) is a given parameter, the problem is in fact underdefined, because there could be an infinite number of (x,y) combinations that result in (A). Nothing seems new there.

But, when a PDE is loosely stated as the stipulation, that ~The derivative of a function of 2 variables~ must equal some other function, of the same 2 variables, that function being referred to as ‘The Test Function’, what one is really saying is that, for any combination of (x) and (y), that Test Function generates *one real number*, and that the derivative must therefore also exist as one real number.

This poses a problem with most PDEs, such as for example, with the Poisson Equation, in which this derivative is simply given as:

```
(dx(u) * dx(v)) + (dy(u) * dy(v))
```

Because here, two gradients are being added, to arrive at one supposed derivative, it’s another form of the first expression I gave above, which was underdefined. Granted, given a field of ‘`u(x,y)`

‘, both the x-axis gradient and the y-axis gradient will follow. But, those remain 2 distinct gradients, which are to be summed, to arrive at one real number per (x,y).

The Dirichlet Problem also states itself to be solvable, as long as the boundary function is “sufficiently smooth”, which I take to mean, that ‘The derivative of the boundary function along the boundary must not exceed some small value.’ Well, FreeFem puts a narrower constraint on the problem, by stating that each boundary must have one exact value. This makes its derivative zero, so that it should qualify as ‘sufficiently smooth’. (:3)

But, what this observation also tells me is that, when I dropped the ‘`dx(u)`

‘ term from the Poisson Equation, to arrive at my degenerate equation above, the Dirichlet Boundary condition resulted in an overdefinition. Because in my second plot above, ~the derivative of the function ‘`u(x,y)`

‘~ was truly dependent on one variable, that being the y-axis gradient, this equation already defined that derivative fully. By adding *2* boundary conditions, I was really opening myself to the possibility, that the problem might become unsolvable.

It was really just, the fact that my explicit boundary conditions for ‘u(x,y)’ coincided at the ‘`Top`

‘ and ‘`Bottom`

‘ with the way I had defined the ‘outer integral’ of my problem, with respect to (y), at both extremes of (y), that no contradiction resulted, and that my overdetermined problem could still be solved.

I take the fact that ‘FreeFem’ *was able to generate an approximate solution*, as independent of the question, of whether *an exact* solution exists in fact.

(Update 7/17/2021, 19h35: )

**3:)**

I have made two discoveries, important to serious use of FreeFem:

- The ‘true’ second-order integral, which I showed a plot of above, the Y-axis lines of which just continuously slope away from the ‘
`Bottom`

‘ boundary, can be achieved efficiently from my first script, by just commenting out Line 64. And, - It’s not actually a fixed behaviour of FreeFem, to impose the boundaries as one value each. Documentation states, that The way the ‘
`int1d()`

‘ function is supposed to behave inside a ‘`solve() ...`

‘ clause is, to impose boundary conditions as a function. In other words, the Dirichlet, ‘`+ on(...)`

‘ statements can also be replaced in some plots, by putting code as follows…

```
- int1d(Th, Bottom) ( g(x, y) )
... ;
```

If ‘`g(real x, real y)`

‘ is a defined function, its returned real numbers would get mapped along the ‘`Bottom`

‘ boundary in the plot. However, ~~I have not tested this feature~~. It’s said to exist, because FreeFem doesn’t *officially* support 1-dimensional integration.

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

**4:)**

Just to make my thinking clear, I’d like to elaborate on how I imagine, FreeFem may perform a first-order integration, of a mesh, which has ‘a number of boundaries’, in general.

What the ‘`int2d()`

‘ function may do is, to set all the vertices that belong to a boundary, to their Dirichlet value, if one is set. This will initially assure, that there is a set of triangles, each of which have 2 vertices’ values out of 3 defined. The 3rd vertex of each triangle is then ~integrated~, by computing the mean of the 2 preceding vertices, and always adding the local value of the function to be integrated in the positive sense, meaning, that if the function happened to be negative, the 3rd vertex value will also result as more-negative.

This is essentially what allows complex meshes to be built with many boundaries.

When the user writes a script which computes a gradient, the onus is on him, to multiply the derivative of the function, which is accurately stated by ‘`dx(u)`

‘ for the derivative with respect to (x) of the Finite Element function (u), with the amount of change in coordinate (x), which is written as ‘`dx(v)`

‘ if the domain of integration was (v), for example. ‘`dx(v)`

‘ will then be negative, if the solver, by way of the ‘`int2d()`

‘ function, stepped from a more-positive element of (v) to a more-negative element of (v), according to the X-axis.

What this means is that, for example, if integration is starting from the ‘`Right`

‘ boundary, where (x) is considered to be positive, leftward, the way the script computes a delta in any function (u) that happens to be positive, is by multiplying a negative step revealed by ‘`dx(v)`

‘, with the eventual negative derivative ‘`dx(u)`

‘, to arrive at a positive value, that will ultimately balance with the positive value of the Test Function.

This solution (u) to the equation may only be reached, after the solver has made more than one pass over the expressions to be solved for. But (u) slopes to more-positive *right-to-left*, from the boundary ‘`Right`

‘, if integration was to start there, and, if the Test Function was positive there. Hence, this is opposite to how classical integration would integrate the definite integral of ‘`f(x,y) = x`

‘, if the indefinite integral at the boundary ‘`Right`

‘, was also the start of the interval of the definite integral.

This just happens to work well, if people want to build complex meshes, with numerous boundaries. Also, which boundaries integration starts from can be modified by stating the Dirichlet values for only some of them, instead of for all of them, or, by setting ‘regions’ to existing boundary labels. Those labels become the non-default starting points of the discrete integration.

Nowhere in my exercises have I worked with FreeFem Regions yet. But it seems, according to the documentation, that their existence is a minor detail, that happens to put labels on triangles (or, on tetrahedrea in the case of ‘`int3d()`

‘), so that regions of either can be integrated in their own subdomains. The fact that this needs to be done in explicit scripting, means that it does not kick in by itself. The ad-hoc integration can still proceed, as was just described.

When, instead of a gradient, ‘`int2d()`

‘ is integrating a linear value such as the Test Function, the script must multiply that value by the area of a quad, which is identified by just ‘`v`

‘ if that is the domain of the integration. And, because the use of cross-products can switch sign over something as trivial, as whether edges were counted clockwise or counter-clockwise, their absolute is used, and in this case, ‘`v`

‘ is always positive. Hence, if the value of the Test Function was negative, the delta will also be negative, regardless of whether progress is right-to-left, or left-to-right…

Hence, when (u) finally balances with the Test Function at the ‘`Left`

‘ boundary, if the Test Function was negative there, a negative derivative of (u) must be multiplying with a positive ‘`dx(v)`

‘.

(Update 7/25/2021, 17h30: )

I’d say that the explanation so far leaves an important observation unexplained: ‘If the sense of the ad-hoc integration is as stated, close to each boundary that has a Dirichlet value, why does the sense reverse near the centre of the plot, into the direction which seems more natural, such that the gradient of (u) agrees with the actual gradient of the Test Function? After all, most users did not program the script, to multiply those two variables.’

And, after much experimentation with the software, I’d say that the reason is the fact that something special happens, when two opposing boundaries are given Dirichlet values…

What I’ve concluded FreeFem does is,

- First, to interpolate its Linear Integral term, which does not depend on (u), for all boundaries that have a Dirichlet value set.
- To save the resulting, first-order integral in a hidden property belonging to (v). And then,
- Generally, to compute the terms ‘
`dx(v)`

‘ and ‘`dy(v)`

‘ as linear combinations of the (Δx) and (Δy) respectively, of the*two*triangle-edges, joining*with already-solved vertices*. And, - The concept should be inferred, that the ‘
`int2d()`

‘ function works with the same methodology, on the LHS of the expression,*recomputing*‘`dx(v)`

‘ and ‘`dy(v)`

‘*each time*. - Further, I have concluded that in order for (u) to have Dirichlet values
*maintained*, each integration which ‘`int2d()`

‘ performs, can more simply start from the value (0.), which will take place on both sides of the expression and lead to an error value of (0.)*at the boundaries themselves*. However, the values of (u) which those error-values accumulate in to, need to be initialized to their Dirichlet values, and then, doing so also needs to take place with an interpolation, such as the one I described below, in case there is more than one such boundary. - When the ‘
`int1d(...) (...)`

‘ term is setting ‘`v`

‘, it must do so with a simple interpolation of boundary values, the weight of which diminishes as the vertices get close to any Boundaries that have a Dirichlet value set.

(Update 7/19/2021, 15h25: )

I can extend my hypothesizing a little further, about what numerical methods the developers might have used (without actual proof), but in a way that will yield the observed results, preserve accuracy as much as possible, and remain as simple as the problem permits.

The form which was used to solve PDEs was ‘`solve PDEName(Vh u, Vh v) = expr;`

‘, where (u) and (v) are Finite Element functions, and (expr) was the expression, which was supposed to be made as close to zero as possible. Previously I had described this form such, that (u) simply receives output values, and that ~(v) was just along for the ride~, doing nothing specifically, other than to provide a vertex-by-vertex map, of how (u) was to be written to. Of course, this would have made (v) redundant.

One fact which certainly exists is that (v) can have properties, which the user has no access to, but which FreeFem uses internally. And, according to my recent ruminations, those properties should include:

(dx), (dy), (wv)

Where ‘`dx(v)`

‘ already states a possibly-interpolated X-coordinate delta, ‘`dy(v)`

‘ states a possibly-interpolated Y-coordinate delta, and according to me, something corresponding to (wv) states, ‘with how much weight the vertex has been written to’, as belonging both to (u) and to (v).

The way I visualize the ad-hoc integration is such, that each boundary replaces values into (v), but with ever-decreasing weight, as the vertices become a number of triangles more distant from the boundary. That weight could just be called (wb). One way to compute it, written in my pseudo-code, might be:

```
int d > 0
wb = 1. / d;
wv += wb;
new_value += wb * contributed_value;
```

Where (d) could be the number of triangles that the current vertex is distanced from the current boundary. This weight effectively reaches zero, ‘directly in front of’ the other boundaries.

Dirk

## One thought on “Implicitly restating a classical, double integral, using FreeFem.”