Unexpected output when summing two fespace variables

General Discussion
1

Here is my sample code

real x0 = 0, x1 = 1.;
real y0 = 0, y1 = 0.5;
int n = 2, m = 1;
mesh Th = square(n, m, [x0 + (x1 - x0) * x, y0 + (y1 - y0) * y]);
mesh uTh = Th;

macro res [rx, ry] // EOM // Residuum.
macro resEl [rxEl, ryEl] // EOM // Residuum.

// *********************
// Finite element space.
// *********************
fespace Vh(Th, [P1, P1]);
// FEM variables.
Vh res, resEl;

rx[] = 1;
ry[] = 2;
rxEl[] = 1;
ryEl[] = 2;

for (int i = 0; i < rx.n; i++){
    rx[](i) = rx[](i) - rx[](i);
    ry[](i) = ry[](i) - ryEl[](i);
}
cout << rx[] << endl;

The final print should be all zeros. Instead it is all “-2”.

12
         -2      -2      -2      -2      -2
         -2      -2      -2      -2      -2
         -2      -2

Anybody understands why the values in rx were not updated?

2

Concerning the output, the initial 12 looks like the number of elements in the array.

With both of the following syntaxes, I get an array of zeros as output.

rx[] = 1;
for [i, rxi : rx[]]{
    rxi = rxi - rxi;
}
cout << rx[] << endl;

rx[] = 1;
for (int i = 0; i < rx[].n; i++) {
    rx[][i] = rx[][i] - rx[][i];
}
cout << rx[] << endl;

Edit: the actual problem was the fact the arrays rx[] and ry[] are the same, and not the element access syntax. See my post below.

3

Concerning why your initial example produces -2 as output.

You are using a vector fespace. If you look the the output of:

cout << "Th.nv = " + Th.nv << endl;
cout << "Vh.ndof = " + Vh.ndof << endl;

Th.nv = 6
Vh.ndof = 12

With P1 elements, the number of degrees of freedom is equal to the number of vertexes. But since you have a vector function with two components it is multipled by two. Now I suspect that the arrays rx[] and ry[] actually point to the same array of size 12 because the idea is their dof’s are coupled.

Your lines were actually doing the following:

rx[] = 1;
ry[] = 2;
rxEl[] = 1;
ryEl[] = 2;
for (int i = 0; i < rx.n; i++){
    rx[](i) = rx[](i) - rx[](i);     // 2 - 2  -->  0
    ry[](i) = ry[](i) - ryEl[](i);   // note ry[][i] points to the same address as as rx[][i],
                                     // 0 - 2  -->  -2
}

Edit: indeed, when you have a vector-valued FE-function, you need to assign both components simultaneously:

[rx, ry]  = [1, 2];
for (int i = 0; i < Vh.ndof; ++i) {
	cout << i << " " << rx[][i] << " " << ry[][i] << endl;
}

This gives as output:

0 1 1
1 2 2
2 1 1
3 2 2
4 1 1
5 2 2
6 1 1
7 2 2
8 1 1
9 2 2
10 1 1
11 2 2

This would seem to indicate that the x and y components of a particular vertex, are stored together in memory.

If you would like to use the array syntax, you could do something like:

int ndof = Vh.ndof;
rx[](0:ndof:2) = 1;
rx[](1:ndof:2) = 2; // actually the y-component

Note, that the indexing appears to exclude ndof (kind of like indexing in Python, i.e. a[0:3] are the elements a[0],a[1],a[2]).