// File Mesh3D.edp
//
// 3D ELLIPTIC EQUATION WITH DIRICHLET AND/OR NEUMANN CONDITIONS
//
// - \sum_{i,j=1}^3 a_{ij} \partial_i\partial_j u + c_0 u = f, \quad x \in \Omega
// u = g,  \quad x \in \Gamma_1
// \partial_n u = h,  \quad x \in \Gamma_2
//

load "msh3"
		
load "medit"
		
//
// NUMERICAL VALUES TO PROVIDE
//

verbosity = 10;

real x0 = 0., x1 = 1.;                          // $x_0$, $x_1$
real y0 = 0., y1 = 2.;                          // $y_0$, $y_1$
real z0 = 0., z1 = .5;                          // $_z0$, $z_1$

int NInt=100, NSigma=99, NBase=98, NTop=97;     // labels: interior, $\Sigma$, $z = 0$, $z = T$
int Nxtot=19, Nytot=39, Nztot=9;                // steps in directions $x$, $y$, $z$
real errmesh=.1;                                // initial error for mesh adaptation
		
real a11=10., a12=0., a13=0., a22=10., a23=0., a33=10.; // pde coefficients
real c0=5.;                                             // pde coefficients

real qx=.15*x0+.85*x1, qy=.15*y0+.85*y1, qz=.5*z0+.5*z1, sigma=.05; // data for rhs

real unsigma=1./sigma, errold;                  // auxiliar
int imesh;                                      // auxiliar

func fseg0=exp(-unsigma*((x-qx)*(x-qx)+(y-qy)*(y-qy)+(z-qz)*(z-qz))); // rhs
func f2d0=exp(-unsigma*((x-qx)*(x-qx)+(y-qy)*(y-qy)));                // rhs at z=0
func gseg0=0.;                                                        // Dirichlet boundary data
func hseg0=0.;                                                        // Neumann boundary data

//
// THE 2D TRIANGULATION
//

border C1(s=x0, x1){x=s; y=y0;     label=NSigma;}
border C2(s=y0, y1){x=x1; y=s;     label=NSigma;}
border C3(s=x1, x0){x=s; y=y1;     label=NSigma;}
border C4(s=y1, y0){x=x0; y=s;     label=NSigma;}
plot(C1(Nxtot)+C2(Nytot)+C3(Nxtot)+C4(Nytot), wait=true);

mesh Th2d=buildmesh(C1(Nxtot)+C2(Nytot)+C3(Nxtot)+C4(Nytot));
plot(Th2d, wait=true);
//savemesh(Th2d,"E-Mesh2D.msh");

errold=errmesh;
for (int i=0;i< 4;i++)
{
	imesh=i+1;
	Th2d=adaptmesh(Th2d,f2d0,err=errold);
	plot(Th2d, cmm="2D MESH NO. "+imesh, wait=true);
	errold = .5*errold;
};

//
// THE 3D TRIANGULATION
//

real px=.5*(x0+x1), py=.5*(y0+y1);
int refint = Th2d(px,py).region;
		
int[int] rup=[refint,NTop,NSigma,NSigma],              // upper face --> NTop
		 rdown=[refint,NBase,NSigma,NSigma],           // lower face --> NBase
		 rmid=[NSigma,NSigma];                         // lateral faces, NSigma--> NSigma

mesh3 Th3d = buildlayers(Th2d, Nztot, zbound=[z0,z1],
		 reffacemid=rmid, reffaceup=rup, reffacelow=rdown);
plot(Th3d,wait=true);
//plot(Th3d,wait=1,ps="mixedTh3d.eps");
//savemesh(Th3d,"E-Mesh3D.mesh");
//medit("3D Mesh", Th3d);

//
// 	THE SPACES
//

fespace  Vh(Th3d,P2);
Vh u,v,fseg=fseg0,gseg=gseg0,hseg=hseg0;

//
// THE PROBLEM
//

problem ellip3d(u,v) = int3d(Th3d)(a11*dx(u)*dx(v) +
		+ a12*(dx(u)*dy(v)+dy(u)*dx(v)) + a13*(dx(u)*dz(v)+dz(u)*dx(v))
		+ a22*dy(u)*dy(v)               + a23*(dy(u)*dz(v)+dz(u)*dy(v))
		+ a33*dx(u)*dx(v)               + c0*u*v)
		- int3d(Th3d)(fseg*v)
        - int2d(Th3d,NSigma)(hseg*v)
        + on(NBase,NTop, u=gseg)
//        + on(NBase,NSigma,NTop, u=gseg)
;

// SOLVING THE PROBLEM

ellip3d;

//
// DISPLAYING, WRITING AND PRINTING
//

// 3D VIEW

plot(Th3d, u, wait=true, value=true, fill=false, cmm="THE SOLUTION");

// ADDITIONAL (VERTICAL) 2D VIEWS

border CC1(s=x0, x1){x=s;  y=0.;    label=NBase;}
border CC2(s=z0, z1){x=x1; y=s;     label=NSigma;}
border CC3(s=x1, x0){x=s;  y=z1;    label=NTop;}
border CC4(s=z1, z0){x=x0; y=s;     label=NSigma;}
//plot(CC1(Nxtot)+CC2(Nztot)+CC3(Nxtot)+CC4(Nztot), wait=true);
plot(CC1(80)+CC2(40)+CC3(80)+CC4(40), wait=true);

//mesh Th2dv=buildmesh(CC1(Nxtot)+CC2(Nztot)+CC3(Nxtot)+CC4(Nztot));
mesh Th2dv=buildmesh(CC1(80)+CC2(40)+CC3(80)+CC4(40));
plot(Th2dv, wait=true, cmm="VERTICAL MESH");

fespace Vh2dv(Th2dv,P2); Vh2dv uxz;

real hy = .1*(y1-y0), yy;
for(int i=1;i<10;i++)
{
	yy = y0+i*hy;
	uxz = u(x,yy,y);
    plot(uxz, cmm="CUT AT y ="+yy, value=true, fill=true, wait=true);
}

// ADDITIONAL (HORIZONTAL) 2D VIEWS

border BB1(s=x0, x1){x=s; y=0.;   label=NSigma;}
border BB2(s=y0, y1){x=x1; y=s;   label=NSigma;}
border BB3(s=x1, x0){x=s; y=y1;   label=NSigma;}
border BB4(s=y1, y0){x=x0; y=s;   label=NSigma;}
//plot(BB1(Nxtot)+BB2(Nytot)+BB3(Nxtot)+BB4(Nytot), wait=true);
plot(BB1(40)+BB2(80)+BB3(40)+BB4(80), wait=true);

//mesh Th2dh=buildmesh(BB1(Nxtot)+BB2(Nytot)+BB3(Nxtot)+BB4(Nytot));
mesh Th2dh=buildmesh(BB1(40)+BB2(80)+BB3(40)+BB4(80));
plot(Th2dh, wait=true);

fespace Vh2dh(Th2dh,P2); Vh2dh uxy;

real hz = .1*(z1-z0), zz;
for(int i=1;i<10;i++)
{
	zz = z0+i*hz;
	uxy = u(x,y,zz);
    plot(uxy, cmm="CUT AT z ="+zz, value=true, fill=true, wait=true);
}

//{ // save solution
//ofstream file("Solution3D.txt");
//file << u[] << endl;
//} // close the file (end block)

// { // read
// ifstream file("f.txt");
// file >> g[] ;
// } // close reading file (end block)