using System;
using Science.Mathematics;
using L=Science.Mathematics.PartialDifferentialEquation;

namespace ScienceTest.MathematicsTest.PartialDifferentialEquationTest
{
	/// <summary>
	/// LaplaceEquationThreeDimensionTest
	/// </summary>
	public class LaplaceEquation3DTest
	{
		public LaplaceEquation3DTest()
		{
		}
		private string result;
		public string Result
		{
			get{return result;}
		}
		public void Compute()
		{
			Function.DoubleDoubleDoubleToDouble[] cf
				= new Function.DoubleDoubleDoubleToDouble[2];
            cf[0] = new Function.DoubleDoubleDoubleToDouble(InnerSphere);
            cf[1] = new Function.DoubleDoubleDoubleToDouble(OuterSphere);
		    L.Domain3D xd = new L.Domain3D();
            xd.ConditionFunctionsLessThanZero = cf;
			xd.LowerBoundOfX = -1.0;
			xd.LowerBoundOfY = -1.0;
            xd.LowerBoundOfZ = -1.0;
            xd.UpperBoundOfX = 1.0;
			xd.UpperBoundOfY = 1.0;
			xd.UpperBoundOfZ = 1.0;

			L.LaplaceEquation3D eq = new L.LaplaceEquation3D();

			eq.MeshSize = 1.0/32.0;
			eq.JacobiRadius = Math.Cos(Math.PI/65.0);

			eq.Domain = xd;
            eq.SetBoundaryValue(-0.11, -0.11, -0.11, 0.11, 0.11, 0.11, 0.0);
            eq.SetBoundaryValue(-1.12, -1.12, -1.12, 1.12, 1.12, - 0.2, 10.0);
            eq.SetBoundaryValue(-1.12, -1.12, -1.12, 1.12, -0.2, 1.12, 10.0);
            eq.SetBoundaryValue(-1.12, -1.12, -1.12, -0.2, 1.12, 1.12, 10.0);
            eq.SetBoundaryValue(0.2, -1.12, -1.12, 1.12, 1.12, 1.12, 10.0);
            eq.SetBoundaryValue(-1.12, 0.2, -1.12, 1.12, 1.12, 1.12, 10.0);
            eq.SetBoundaryValue(-1.12, -1.12, 0.2, 1.12, 1.12, 1.12, 10.0);
			
	    	eq.Solve();

  /*        for (int i = 0; i < eq.Solution.GetLength(0); i++)
            {
                for (int j = 0; j < eq.Solution.GetLength(2) - 1; j++)
                    result += Convert.ToString(eq.Solution[i, 32, j]) + "\t";
                result += Convert.ToString(eq.Solution[i, 32, eq.Solution.GetLength(1) - 1]) + "\r\n";
            } 	*/

            L.EqualValueSurface surface = new L.EqualValueSurface();
            surface.Solution3D = eq.Solution;
            surface.EqualValue = 5.0;
            surface.FindSurface();

			for(int i = 0; i < surface.SurfacePointX.Count; i++)
			{
                result += Convert.ToString(xd.LowerBoundOfX+eq.MeshSize*Convert.ToDouble(surface.SurfacePointX[i])) + "\t";
                result += Convert.ToString(xd.LowerBoundOfY+eq.MeshSize*Convert.ToDouble(surface.SurfacePointY[i])) + "\t";
                result += Convert.ToString(xd.LowerBoundOfX+eq.MeshSize*Convert.ToDouble(surface.SurfacePointZ[i])) + "\r\n";
			} 	
		}
				
		private double InnerSphere(double x, double y, double z)
		{
			return -x*x-y*y-z*z+0.1*0.1;
		}
		private double OuterSphere(double x, double y, double z)
		{
			return x*x+y*y+z*z-0.8*0.8;
		}		
	}
}