real a = 40, b = 40, c = 0.5; border a00(t=0, 1) {x=a*t; y=0; label=1;} border a10(t=0, 1) {x=a; y=b*t; label=1;} border a20(t=1, 0) {x=a*t; y=b; label=1;} border a30(t=1, 0) {x=0; y=b*t; label=1;} border a01(t=0, 1) {x=c+(a-c*2)*t; y=c; label=1;} border a11(t=0, 1) {x=a-c; y=c+(b-c*2)*t; label=1;} border a21(t=1, 0) {x=c+(a-c*2)*t; y=b-c; label=1;} border a31(t=1, 0) {x=c; y=c+(b-c*2)*t; label=1;} real p = 5, q = 20, d = 34, e = 1; border b00(t=0, 1) {x=p+d*t; y=q; label=3;} border b10(t=0, 1) {x=p+d; y=q+e*t; label=3;} border b20(t=1, 0) {x=p+d*t; y=q+e; label=3;} border b30(t=1, 0) {x=p; y=q+e*t; label=3;} real r = 30, s =1 , j = 1, u = 15; border c00(t=0, 1) {x=r+j*t; y=s; label=3;} border c10(t=0, 1) {x=r+j; y=s+u*t; label=3;} border c20(t=1, 0) {x=r+j*t; y=s+u; label=3;} border c30(t=1, 0) {x=r; y=s+u*t; label=3;} int n=13; mesh Sh = buildmesh(a00(10*n) + a10(10*n) + a20(10*n) + a30(10*n) + a01(10*n) + a11(10*n) + a21(10*n) + a31(10*n) + b00(5*n) + b10(5*n) + b20(5*n) + b30(5*n) + c00(5*n) + c10(5*n) + c20(5*n) + c30(5*n)); plot(Sh, wait=1); int bx; for (bx = 1; bx <= 7; bx++){ border C(t=0, 2*pi){x=2+cos(t); y=bx*5+sin(t); label=2;} mesh Th = buildmesh(a00(10*n) + a10(10*n) + a20(10*n) + a30(10*n) + a01(10*n) + a11(10*n) + a21(10*n) + a31(10*n) + C(10) + b00(5*n) + b10(5*n) + b20(5*n) + b30(5*n) + c00(5*n) + c10(5*n) + c20(5*n) + c30(5*n)); fespace Vh(Th, P1); func real wall() { if (Th(x,y).region == Th(0.5,0.5).region || Th(x,y).region == Th(7,20.5).region || Th(x,y).region == Th(30.5,2).region) { return 1; } else { return 0; } } Vh v,w; randinit(900); Vh wallreflexion = randreal1(); Vh wallabsorption = randreal1()*0.5i; Vh k = 6; cout << "Reflexion of walls min/max: " << wallreflexion[].min << " " << wallreflexion[].max << "\n"; cout << "Absorption of walls : " << wallabsorption[].min << " "<< wallabsorption[].max << "\n"; problem muwave(v,w) = int2d(Th)( (v*w*k^2)/(1+(wallreflexion+wallabsorption)*wall())^2 - (dx(v)*dx(w)+dy(v)*dy(w)) ) + on(2, v=1) ; muwave; Vh vm = log(real(v)^2 + imag(v)^2); plot(vm, wait=1, fill=true, value=0, nbiso=65); }