// file heatexample.edp
//
// solution of $-\nabla\cdot(\kappa(x,y)\nabla u) = 0$
// together with $u = u_0$ on $\Gamma_0$
//               $u = u_1$ on $\Gamma_1$
//
// $\kappa$ is piecewise constant
//


//
// NUMERICAL VALUES TO PROVIDE
//
int C0=100, C1=99, C2=98; // the labels
real u0=20., u1=100.;     // the values on the boundaries
real k0=1., k1=0.001;     // the values of $\kappa$

real k1mk0=k1-k0; // auxiliar

//
// DESCRIPTION OF THE BOUNDARY AND THE INTERFACE
//
// outer boundary ($\Gamma_0$):
border C00(t=0,2*pi){x=5*cos(t); y=5*sin(t); label=C0;}

// inner boundary ($\Gamma_1$):
border C11(t=0,1){ x=1+t;  y=3;      label=C1;}
border C12(t=0,1){ x=2;    y=3-6*t;  label=C1;}
border C13(t=0,1){ x=2-t;  y=-3;     label=C1;}
border C14(t=0,1){ x=1;    y=-3+6*t; label=C1;}

// interface ($\Gamma_2$):
border C21(t=0,1){ x=-2+t; y=3;      label=C2;}
border C22(t=0,1){ x=-1;   y=3-6*t;  label=C2;}
border C23(t=0,1){ x=-1-t; y=-3;     label=C2;}
border C24(t=0,1){ x=-2;   y=-3+6*t; label=C2;}

//
// FIGURES (I)
//
// boundary and interface:
plot (                C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20)
                    + C21(15)+C22(60)+C23(15)+C24(60),
                    wait=true, ps="heatexb.eps");

//
// MESH
//
mesh Th=buildmesh(    C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20)
                    + C21(-15)+C22(-60)+C23(-15)+C24(-60));

//
// FIGURES (II)
//
// mesh:
plot(Th, wait=true, ps="heatexTh.eps");

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//
fespace Vh(Th,P1); Vh u,v;
Vh kappa = k0 + k1mk0*(x<-1)*(x>-2)*(y<3)*(y>-3);

plot(kappa, wait=true, value=true, fill=true);

//
// THE PROBLEM (FORMULATION AND SOLUTION)
//
problem a(u,v) = int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
        +on(C0,u=u0) + on(C1,u=u1);
a;

//
// FIGURES (III)
//
// the solution:
plot(u, wait=true, value=true, fill=true, ps="heatex.eps");