//The mesh is uniform with 240  60  30 points.
double N = 15;
vector n = (8*N, 2*N, N);
vector a = (0,0,0);
//vector b = (25000,695*L,6000);
vector b = (25000,695,6000);
scene S=pov("dechet2.pov");
mesh Mesh = structured(n,a,b);
domain D= domain(S,outside(<1,0,0>));
function dogger = one(y<=199.99);
function clay =one(y <= (295.01+ 55.01*x/25000))*one(y>199.99);
function
limestone=one(y<=594.99)*one(y>(295.01+55.01*x/25000));
function marl = one(y>594.99);
function K = (3.1536e-5 * marl + 6.3072 * limestone
                        + 3.1536e-6*clay + 25.2288*dogger);

solve(H) in D by Mesh
        
cg(maxiter=900,epsilon=1E-10),krylov(precond=ichol)

{
 pde(H) 
        -div(K*grad(H))=0;
        dnu(H)=0 on <1,0,0>;
        H =289* dogger + 310* limestone
                + clay* (289 + (y-200)*21.0/150)
                + marl * (310 + (y - 595)*0.3) on Mesh xmax;
        H =286*dogger + 200*limestone
                + clay * (286 - (y-200)*86/95)
                + marl * (200 - (y - 595)*0.2) on Mesh xmin;
        H =180 + 160 * x/25000 on Mesh ymax;
};

save(medit,"H",H,Mesh,dos);
save(medit,"H",Mesh,dos);