// file prob-n-doors.edp

// Probability of reaching one, two, three or four doors
//
// solution of $-\Delta p = 0$
// together with $p = 0$ on the walls ($B_1$)
//               $p = 1$ on the door(s) ($B_2$)
//


//
// NUMERICAL VALUES TO PROVIDE
//
int B1=99, B2=98; // the labels of the walls and the doors

//
// ONE DOOR
//

border C1(t=0,4.75){x=t;     y=0;    label=B1;}
border C2(t=4.75,5.25){x=t;  y=0;    label=B2;}
border C3(t=5.25,10){x=t;    y=0;    label=B1;}
border C4(t=0,10){x=10;      y=t;    label=B1;}
border C5(t=0,10){x=10-t;    y=10;   label=B1;}
border C6(t=0,10){x=0;       y=10-t; label=B1;}

plot (                C1(19)+C2(3)+C3(19)
                    + C4(41)+C5(41)+C6(41),
                    wait=true, ps="proba1b.eps");

mesh Th=buildmesh(    C1(16)+C2(3)+C3(16)
                    + C4(41)+C5(41)+C6(41));
plot(Th,wait=true,ps="proba1Th.eps");

fespace Vh(Th,P2); Vh u,up;

solve a(u,up)= int2d(Th)(dx(u)*dx(up)+dy(u)*dy(up))
                +on(B1,u=0)+on(B2,u=1);
// plot(u,wait=true, value=true, ps="proba1.eps");
plot(u, wait=true, value=true, fill=true, ps="proba1.eps");

//
// TWO DOORS
//

border K1(t=0,4.75){x=t;     y=0;    label=B1;}
border K2(t=4.75,5.25){x=t;  y=0;    label=B2;}
border K3(t=5.25,10){x=t;    y=0;    label=B1;}
border K4(t=0,4.75){x=10;    y=t;    label=B1;}
border K5(t=4.75,5.25){x=10; y=t;    label=B2;}
border K6(t=5.25,10){x=10;   y=t;    label=B1;}
border K7(t=0,10){x=10-t;    y=10;   label=B1;}
border K8(t=0,10){x=0;       y=10-t; label=B1;}

plot (                K1(19)+K2(3)+K3(19)
                    + K4(19)+K5(3)+K6(19)
                    + K7(41)+K8(41),
                    wait=true, ps="proba2b.eps");

mesh Sh=buildmesh(    K1(19)+K2(3)+K3(19)
                    + K4(19)+K5(3)+K6(19)
                    + K7(41)+K8(41));
plot(Sh,wait=true,ps="proba2Th.eps");

fespace Wh(Sh,P2); Wh w,wp;

solve b(w,wp)= int2d(Sh)(dx(w)*dx(wp)+dy(w)*dy(wp))
                +on(B1,w=0)+on(B2,w=1);
// plot(w,wait=true, value=true, ps="proba2.eps");
plot(w, wait=true, value=true, fill=true, ps="proba2.eps");

//
// THREE DOORS
//

border L1(t=0,4.75){x=t;     y=0;    label=B1;}
border L2(t=4.75,5.25){x=t;  y=0;    label=B2;}
border L3(t=5.25,10){x=t;    y=0;    label=B1;}
border L4(t=0,4.75){x=10;    y=t;    label=B1;}
border L5(t=4.75,5.25){x=10; y=t;    label=B2;}
border L6(t=5.25,10){x=10;   y=t;    label=B1;}
border L7(t=0,4.75){x=10-t;  y=10;   label=B1;}
border L8(t=4.75,5.25){x=10-t;y=10;  label=B2;}
border L9(t=5.25,10){x=10-t; y=10;   label=B1;}
border L10(t=0,10){x=0;      y=10-t; label=B1;}

plot (                L1(19)+L2(3)+L3(19)
                    + L4(19)+L5(3)+L6(19)
                    + L7(19)+L8(3)+L9(19)
                    + L10(41),
                    wait=true, ps="proba3b.eps");

mesh Rh=buildmesh(    L1(19)+L2(3)+L3(19)
                    + L4(19)+L5(3)+L6(19)
                    + L7(19)+L8(3)+L9(19)
                    + L10(41));
plot(Rh,wait=true,ps="proba3Th.eps");

fespace Zh(Rh,P2); Zh p,q;

solve c(p,q)= int2d(Rh)(dx(p)*dx(q)+dy(p)*dy(q))
                +on(B1,p=0)+on(B2,p=1);
// plot(p, wait=true, value=true, ps="proba3.eps");
plot(p, wait=true, value=true, fill=true, ps="proba3.eps");

//
// FOUR DOORS
//

border M1(t=0,4.75){x=t;     y=0;    label=B1;}
border M2(t=4.75,5.25){x=t;  y=0;    label=B2;}
border M3(t=5.25,10){x=t;    y=0;    label=B1;}
border M4(t=0,4.75){x=10;    y=t;    label=B1;}
border M5(t=4.75,5.25){x=10; y=t;    label=B2;}
border M6(t=5.25,10){x=10;   y=t;    label=B1;}
border M7(t=0,4.75){x=10-t;  y=10;   label=B1;}
border M8(t=4.75,5.25){x=10-t;y=10;  label=B2;}
border M9(t=5.25,10){x=10-t; y=10;   label=B1;}
border M10(t=0,4.75){x=0;    y=10-t; label=B1;}
border M11(t=4.75,5.25){x=0; y=10-t; label=B2;}
border M12(t=5.25,10){x=0;   y=10-t; label=B1;}

plot (                M1(19)+M2(3)+M3(19)
                    + M4(19)+M5(3)+M6(19)
                    + M7(19)+M8(3)+M9(19)
                    + M10(19)+M11(3)+M12(19),
                    wait=true, ps="proba4b.eps");

mesh Qh=buildmesh(    M1(19)+M2(3)+M3(19)
                    + M4(19)+M5(3)+M6(19)
                    + M7(19)+M8(3)+M9(19)
                    + M10(19)+M11(3)+M12(19));
plot(Qh,wait=true,ps="proba4Th.eps");

fespace Yh(Qh,P2); Yh v, vp;

solve d(v,vp)= int2d(Qh)(dx(v)*dx(vp)+dy(v)*dy(vp))
                +on(B1,v=0)+on(B2,v=1);
//plot(v,wait=true, value=true, ps="proba4.eps");
plot(v, wait=true, value=true, fill=true, ps="proba4.eps");