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");
}
}
}
}