Thank you, Chris. What I would like to build is essentially a step function that I can easily change how many steps it has with one parameter. For example, how can I build something like:

f = 0, x < 0 and x > J
i, i - 1< x < i, i = 1, …, J.

I’d like to build it in a loop where I can easily change how many steps there are with a simple change of the parameter J. Especially since sometimes I want many steps.

Yes, I have, and I’ve tried many things. I can’t seem to get it to run properly though. For example, I can get a running code doing the following, but f doesn’t plot like g.

int J = 6;
func f = 0;
for(int i = 0; i < J; i+=2){
cout << "i = " << i << endl;
func h = f + (x > i && x < i+1 && y > i && y < i+1);
func f = h;
}
func g = (x > 0 && x < 1 && y > 0 && y < 1) + (x > 2 && x < 3 && y > 2 && y < 3) + (x > 4 && x < 5 && y > 4 && y < 5);
border B1(t = 0, 10){x = t; y = 0; label = 1;};
border B2(t = 0, 10){x = 10; y = t; label = 2;};
border B3(t = 0, 10){x = 10-t; y = 10; label = 3;};
border B4(t = 0, 10){x = 0; y = 10-t; label = 4;};
mesh Th = buildmesh(B1(50) + B2(50) + B3(50) + B4(50));
fespace Vh(Th, P1);
Vh gproj = g, fproj = f;
plot(fproj, wait = true, cmm = "fproj");
plot(gproj, Th, wait = true, cmm = "gproj");

I want to be able to build the function in the loop because I want to be able to make J as large or small as I want without having to manually add or subtract terms.

Thank you, I think I see your problem now. Here’s a version of your example that works. Take a look and let me know if you have any questions.

int J = 6;
/* This section of code does not make sense
func f = 0;
for(int i = 0; i < J; i+=2){
cout << "i = " << i << endl;
func h = f + (x > i && x < i+1 && y > i && y < i+1);
func f = h;
}
*/
// Maybe this is what you're looking for?
func real f(int J){
real h;
for (int i = 0; i < J; i+=2){
//cout << "i = " << i << endl; This will print a lot if kept
h += (x > i && x < i+1 && y > i && y < i+1);
}
return h;
}
func g = (x > 0 && x < 1 && y > 0 && y < 1) + (x > 2 && x < 3 && y > 2 && y < 3) + (x > 4 && x < 5 && y > 4 && y < 5);
//border B1(t = 0, 10){x = t; y = 0; label = 1;};
//border B2(t = 0, 10){x = 10; y = t; label = 2;};
//border B3(t = 0, 10){x = 10-t; y = 10; label = 3;};
//border B4(t = 0, 10){x = 0; y = 10-t; label = 4;};
//mesh Th = buildmesh(B1(50) + B2(50) + B3(50) + B4(50));
// Here's an easier option to mesh the 10x10 square:
mesh Th = square(50,50,[10*x,10*y]);
fespace Vh(Th, P1);
Vh gproj = g, fproj = f(J);
//plot(fproj, wait = true, cmm = "fproj", fill = 1);
//plot(gproj, Th, wait = true, cmm = "gproj", fill = 1);
// Let's compare the difference instead
Vh diff = fproj-gproj;
plot(diff, wait = true, cmm = "fproj-gproj", fill = 1);