# Unexpected output when summing two fespace variables

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?

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.

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]`).