using System;
using System.Collections.Generic;

namespace Science.Mathematics.PartialDifferentialEquation
{
	/// <summary>
	/// LinearEllipticEquationTwoDimension
    /// C_xxU_xx+C_xU_x+C_yyU_yy+C_yU_y+C_zzU_zz+C_zU_z
    /// +C_xyU_xy+C_xzU_xz+C_yzU_yz+C_propU=C_source
	/// </summary>
	public class LinearEllipticEquation3D
	{
		public delegate double Coefficient(double x, double y, double z);
		
		public LinearEllipticEquation3D()
		{
		}
	
		private double delta, rj;
		public double MeshSize
		{
			get{return delta;}
			set{delta=value;}
		}
		public double JacobiRadius
		{
			get{return rj;}
			set{rj=value;}
		}

		private Domain3D dtd;
		public Domain3D Domain
		{
    		set{dtd=value;}
		}

		private Coefficient cxx, cx, cyy, cy, czz, cz, cxy, cxz, cyz, cprop, csource;
		public Coefficient CoefficientOfSecondDerivativeOfX
		{
			set{cxx=value;}
		}
		public Coefficient CoefficientOfFirstDerivativeOfX
		{
			set{cx=value;}
		}
		public Coefficient CoefficientOfSecondDerivativeOfY
		{
			set{cyy=value;}
		}
		public Coefficient CoefficientOfFirstDerivativeOfY
		{
			set{cy=value;}
		}
        public Coefficient CoefficientOfSecondDerivativeOfZ
        {
            set { czz = value; }
        }
        public Coefficient CoefficientOfFirstDerivativeOfZ
        {
            set { cz = value; }
        }
		public Coefficient CoefficientOfDerivativeOfXandY
		{
			set{cxy=value;}
		}
        public Coefficient CoefficientOfDerivativeOfXandZ
        {
            set { cxz = value; }
        }
        public Coefficient CoefficientOfDerivativeOfYandZ
        {
            set { cyz = value; }
        }
		public Coefficient CoefficientOfProportional
		{
			set{cprop=value;}
		}
		public Coefficient CoefficientOfSourceTerm
		{
			set{csource=value;}
		}

        private List<double> al = new List<double>();
        public void SetBoundaryValue(double xFrom, double yFrom, double zFrom,
            double xTo, double yTo, double zTo, double boundaryValue)
        {
            al.Add(xFrom);
            al.Add(yFrom);
            al.Add(zFrom);
            al.Add(xTo);
            al.Add(yTo);
            al.Add(zTo);
            al.Add(boundaryValue);
        }


		private bool[,,] isItInDomain;
		private double[,,] bv;
		public void Solve()
		{
            double x, y, z;
            int nbv = al.Count / 7;
            double xf = dtd.LowerBoundOfX;
            double xt = dtd.UpperBoundOfX;
            double yf = dtd.LowerBoundOfY;
            double yt = dtd.UpperBoundOfY;
            double zf = dtd.LowerBoundOfZ;
            double zt = dtd.UpperBoundOfZ;
            int jx = (int)((xt - xf) / delta) + 1;
            int jy = (int)((yt - yf) / delta) + 1;
            int jz = (int)((zt - zf) / delta) + 1;

            double[, ,] axp = new double[jx, jy, jz];
            double[, ,] axm = new double[jx, jy, jz];
            double[, ,] ayp = new double[jx, jy, jz];
            double[, ,] aym = new double[jx, jy, jz];
            double[, ,] azp = new double[jx, jy, jz];
            double[, ,] azm = new double[jx, jy, jz];
            double[, ,] axy = new double[jx, jy, jz];
            double[, ,] axz = new double[jx, jy, jz];
            double[, ,] ayz = new double[jx, jy, jz];
            double[, ,] aprop = new double[jx, jy, jz];
            double[, ,] asource = new double[jx, jy, jz];
			
			isItInDomain = new bool[jx,jy,jz];
			bv = new double[jx,jy,jz];
			solution = new double[jx,jy,jz];

			for(int k = 0; k < jx; k++)
			{
				x = xf + k*delta;
				for(int m = 0; m < jy; m++)
				{
					y = yf + m*delta;
                    for (int n = 0; n < jz; n++)
                    {
                        z = zf + n*delta;

                        axp[k, m, n] = cxx(x, y, z) + cx(x, y, z) * delta / 2.0;
                        axm[k, m, n] = cxx(x, y, z) - cx(x, y, z) * delta / 2.0;
                        ayp[k, m, n] = cyy(x, y, z) + cy(x, y, z) * delta / 2.0;
                        aym[k, m, n] = cyy(x, y, z) - cy(x, y, z) * delta / 2.0;
                        azp[k, m, n] = czz(x, y, z) + cz(x, y, z) * delta / 2.0;
                        azm[k, m, n] = czz(x, y, z) - cz(x, y, z) * delta / 2.0;

                        aprop[k, m, n] = - 2.0 * cxx(x, y, z) - 2.0 * cyy(x, y, z) - 2.0 * czz(x, y, z) + cprop(x, y, z) * delta * delta;
                        asource[k, m, n] = csource(x, y, z) * delta * delta;

                        axy[k, m, n] = cxy(x, y, z) / 4.0;
                        axz[k, m, n] = cxz(x, y, z) / 4.0;
                        ayz[k, m, n] = cyz(x, y, z) / 4.0;

                        isItInDomain[k, m, n] = dtd.Check(x, y, z);
                        if (!isItInDomain[k, m, n])
                        {
                            for (int kk = 0; kk < nbv; kk++)
                            {
                                if (Convert.ToDouble(al[7 * kk]) < x && x < Convert.ToDouble(al[7 * kk + 3])
                                    && Convert.ToDouble(al[7 * kk + 1]) < y && y < Convert.ToDouble(al[7 * kk + 4])
                                    && Convert.ToDouble(al[7 * kk + 2]) < z && z < Convert.ToDouble(al[7 * kk + 5]))
                                    bv[k, m, n] = Convert.ToDouble(al[7 * kk + 6]);
                            }
                        }                       
                    }
				}
			}
            for (int k = 0; k < jx; k++)
            {
                for (int m = 0; m < jy; m++)
                {
                    for (int n = 0; n < jz; n++)
                    {
                        solution[k, m, n] = bv[k, m, n];
                    }
                }
            }
            sor(axp,axm,ayp,aym,azp,azm,aprop,asource,axy,axz,ayz,solution, rj);
		}

        private double[, ,] solution;
        public double[, ,] Solution
        {
            get { return solution; }
        }
	

		// numerical recipes

		private int MAXITS = 1000;
		private double EPS = 1.0e-5;

		private void sor(double[,,]axp,double[,,]axm,double[,,]ayp,double[,,]aym,
            double[,,]azp,double[,,]azm,double[,,]aprop,double[,,]asource,
            double[, ,] axy, double[, ,] axz, double[, ,] ayz, double[, ,] u, double rjac)
		{
			int jmax = u.GetLength(0);
			int ipass,j,jsw,l,lsw,n,m;
			double anorm = 0.0, omega=1.0, resid;

            for (n = 1; n <= MAXITS; n++)
            {
                anorm = 0.0;
                jsw = 1;
                for (ipass = 1; ipass <= 2; ipass++)
                {
                    lsw = jsw;
                    for (j = 1; j < jmax - 1; j++)
                    {
                        for (m = 1; m < jmax - 1; m++)
                        {
                            for (l = lsw; l < jmax - 1; l += 2)
                            {
                                if (isItInDomain[j, m, l])
                                {
                                    resid = axp[j, m, l] * u[j + 1, m, l]
                                          + axm[j, m, l] * u[j - 1, m, l]
                                          + ayp[j, m, l] * u[j, m + 1, l]
                                          + aym[j, m, l] * u[j, m - 1, l]
                                          + azp[j, m, l] * u[j, m, l + 1]
                                          + azm[j, m, l] * u[j, m, l - 1]
                                          + aprop[j, m, l] * u[j, m, l]
                                          - asource[j, m, l]
                                          + axy[j, m, l] * (u[j + 1, m + 1, l] + u[j - 1, m - 1, l] - u[j + 1, m - 1, l] - u[j - 1, m + 1, l])
                                          + axz[j, m, l] * (u[j + 1, m, l + 1] + u[j - 1, m, l - 1] - u[j + 1, m, l - 1] - u[j - 1, m, l + 1])
                                          + ayz[j, m, l] * (u[j, m + 1, l + 1] + u[j, m - 1, l - 1] - u[j, m + 1, l - 1] - u[j, m - 1, l + 1]);
                                    anorm += Math.Abs(resid);
                                    u[j, m, l] += omega * resid / 6.0;
                                }
                                else u[j, m, l] = bv[j, m, l];
                            }
                            lsw = 3 - lsw;
                        }
                    }
                    jsw=3-jsw;
					omega=(n == 1 && ipass == 1 ? 1.0/(1.0-0.5*rjac*rjac) :
						1.0/(1.0-0.25*rjac*rjac*omega));
				}
				if (anorm < EPS) return; 
			}
            try { throw new Science.Error(); }
            catch (Science.Error er)
            {
                er.Write("MAXITS exceeded");
            }			
		}
	}
}