// file simpleexample.edp
//
// solution of $-D\Delta u + cu = f$
// together with appropriate boundary conditions

//
// NUMERICAL VALUES TO PROVIDE
//

int n=20, m=20;                         // steps
//int ig1=101, ig2=102, ig3=103, ig4=104; // the boundary labels
real x0=0., x1=5., y0=-2.5, y1=2.5;     // the corners
real D=.5, c=1., alp=1., f0=10.;     // coefficients and rhs data

//
// THE RIGHT HAND SIDE (f IS f0 NEAR (3.5,.5), ZERO OTHERWISE)
//

func fseg=f0*(abs(x-3.5) < 1)*(abs(y-.5) < 1);

//
// DESCRIPTION OF THE DOMAIN (A SQUARE) AND MESH:
//

mesh Th=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

//
// FIGURES
//

//int[int] labs=[ig1,ig2,ig3,ig4];
for (int i=0;i<5;++i)
{
	mesh Sh=square(n,m,flags=i);
	plot(Sh,wait=1,cmm=" square flags = "+i );
}

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//

fespace Vh(Th,P1); Vh u, v, f=fseg;
//fespace Vh(Th,P2); Vh u, v, f=fseg;

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

//
// THE PROBLEMS (FORMULATIONS, SOLUTIONS AND FIGURES)
//
// DIRICHLET
//

problem aDir(u,v)=int2d(Th)(c*u*v+D*(dx(u)*dx(v)+dy(u)*dy(v)))
        -int2d(Th)(f*v)
        +on(1,u=0.)+on(2,u=0.)+on(3,u=0.)+on(4,u=0.);
aDir;
plot(u, wait=true, value=true, fill=true, cmm="DIRICHLET");

//
// NEUMANN
//

problem aNeu(u,v) = int2d(Th)(c*u*v+D*(dx(u)*dx(v)+dy(u)*dy(v)))
        -int2d(Th)(f*v);
aNeu;
plot(u, wait=true, value=true, fill=true, cmm="NEUMANN");

//
// DIRICHLET-NEUMANN
//

problem aDNeu(u,v) = int2d(Th)(c*u*v+D*(dx(u)*dx(v)+dy(u)*dy(v)))
        -int2d(Th)(f*v)
        +on(1,u=0.)+on(3,u=0.);
aDNeu;
plot(u, wait=true, value=true, fill=true, cmm="DIRICHLET-NEUMANN");

//
// DIRICHLET-FOURIER
//

problem aDFou(u,v) = int2d(Th)(c*u*v+D*(dx(u)*dx(v)+dy(u)*dy(v)))
        +int1d(Th,2)(alp*u*v)+int1d(Th,4)(alp*u*v)
        -int2d(Th)(f*v)
        +on(1,u=0.)+on(3,u=0.);
aDFou;
plot(u, wait=true, value=true, fill=true, cmm="DIRICHLET-FOURIER");