Defining by parts properties

General Discussion
1

Hello fellows!

I need to define different properties values for three or more regions in my domain.
Following the FreeFem documentation, I was able to build the mesh in a way that this regions are known.

Searching in this community I found how to do that. But, only the code for two differents properties works. When I try the code for three, the result is odd numbers (not necessary related to the ones I gave).

Let’s to the code

Building the mesh

int n=20;
 real x0 = 0.0, x1 = 0.5, y0 = 0.0, y1 = 0.5; 
 real r = 0.15, cx = x1/2, cy = y1/2; 
 real rd = r/5; 

 border Q01(t=0, 1){x=x1*t; y=y0; label=1;}
 border Q03(t=0, 1){x=x1-x1*t; y=y1; label=3;}
 border Q04(t=0, 1){x=x0; y=y1-y1*t; label=4;}
 border Q02(t=0, 1){x=x1; y=y1*t; label=2;}

border C01(t=0, 2*pi){x=cx+r*cos(t); y=cy+r*sin(t); label=5;}
border C02(t=0, 2*pi){x=cx+(rd)*cos(t); y=cy+(rd)*sin(t); label=6;}

Th0 = buildmesh (Q01(n) + Q02(n) + Q03(n) + Q04(n) + C01(-n));
Th1 = buildmesh (C01(n) + C02(-n)); 
Th2 = buildmesh (C02(n)); 
Th = buildmesh (Q01(n) + Q02(n) + Q03(n) + Q04(n) + C01(n) + C02(n));

Defining by parts property (this code works just fine, but it is only for two regions)

real reg1 = Th(0,0).region; 
real reg2 = Th(0.35,0.35).region; 
real reg3 = Th(0.25,0.25).region; 
vk[0] = 1.0e-4; vk[1] = 0.12; vk[2] = 0.5;
fespace Ph(Th, P0);
Ph k = (region == reg1 ? vk[0] : vk[1]);

The code for three regions is the below (which is giving some strange result)

real reg1 = Th(0,0).region; 
real reg2 = Th(0.35,0.35).region; 
real reg3 = Th(0.25,0.25).region; 
vk[0] = 1.0e-4; vk[1] = 0.12; vk[2] = 0.5;
fespace Ph(Th, P0);
Ph k = vk[0]*(region==reg1) + vk[1]*(region==reg2) + vk[2]*(region=reg3);

Can someone give me a help?

Thanks in advance!

2

For reg3 you use assignment (=) instead of comparison (==).

Although it doesn’t seem to cause a problem I would expect region numbers to be int rather than real.

1 Like
3

Good evening!
There is only a little mistake within your script. Indeed in the last line you shared you wrote:
… (region=reg3), but of course you are required to use a proper logical operator which is in this case the == sign! So your last line should become something like:

Ph k = vk[0]*(region == reg1) + vk[1]*(region == reg2) + vk[2]*(region == reg3);

In such a way if you access your property in the region you defined you find the correct definition:

cout<<"k(reg1) = "<<k(0,0)<<endl;
cout<<"k(reg2) = "<<k(0.35,0.35)<<endl;
cout<<"k(reg3) = "<<k(0.25,0.25)<<endl;
4

Thank you so much @gero and @marcodeba !

I’m so embarrassed rn.
Since yesterday I haven’t be able to find the error in this code.
And it was dancing in front of me.

Thanks for the nice response.