using System; using System.Collections.Generic; namespace Science.Mathematics.PartialDifferentialEquation { /// /// LaplaceEquationThreeDimension /// u_xx+u_yy+u_zz=\rho /// public class PoissonEquation3D { public delegate double Source(double x, double y, double z); public PoissonEquation3D() { } 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 Source source; public Source SourceTerm { set{source=value;} } private List al = new List(); 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[,,] f = 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; f[k, m, n] = source(x, y, z) * delta * delta; 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(f,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[,,] f, 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 = u[j + 1, m, l] + u[j - 1, m, l] + u[j, m + 1, l] + u[j, m - 1, l] + u[j, m, l + 1] + u[j, m, l - 1] - 6.0 * u[j, m, l] - f[j,m,l]; 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"); } } } }