using System;
using System.Collections.Generic;
namespace Science.Mathematics.PartialDifferentialEquation
{
///
/// LaplaceEquationTwoDimension
/// u_xx+u_yy=0
///
public class LaplaceEquation2D
{
public LaplaceEquation2D()
{
}
private double delta, rj;
public double MeshSize
{
get{return delta;}
set{delta=value;}
}
public double JacobiRadius
{
get{return rj;}
set{rj=value;}
}
private Domain2D dtd;
public Domain2D Domain
{
set{dtd=value;}
}
private List al = new List();
public void SetBoundaryValue(double xFrom, double yFrom,
double xTo, double yTo, double boundaryValue)
{
al.Add(xFrom);
al.Add(yFrom);
al.Add(xTo);
al.Add(yTo);
al.Add(boundaryValue);
}
private bool[,] isItInDomain;
private double[,] bv;
public void Solve()
{
double x,y;
int nbv = al.Count/5;
double xf = dtd.LowerBoundOfX;
double xt = dtd.UpperBoundOfX;
double yf = dtd.LowerBoundOfY;
double yt = dtd.UpperBoundOfY;
int jx = (int) ((xt - xf) / delta) + 1;
int jy = (int) ((yt - yf) / delta) + 1;
double[,] a = new double[jx,jy];
double[,] b = new double[jx,jy];
double[,] c = new double[jx,jy];
double[,] d = new double[jx,jy];
double[,] e = new double[jx,jy];
isItInDomain = new bool[jx,jy];
bv = new double[jx,jy];
solution = new double[jx,jy];
for(int k = 0; k < jx; k++)
{
x = xf + k*delta;
for(int m = 0; m < jy; m++)
{
y = yf + m*delta;
a[k,m]=1.0;
b[k,m]=1.0;
c[k,m]=1.0;
d[k,m]=1.0;
e[k,m]=-4.0;
isItInDomain[k,m]=dtd.Check(x,y);
if(!isItInDomain[k,m])
{
for(int kk = 0; kk < nbv; kk++)
{
if(Convert.ToDouble(al[5*kk]) < x && x < Convert.ToDouble(al[5*kk+2])
&& Convert.ToDouble(al[5*kk+1]) < y && y < Convert.ToDouble(al[5*kk+3]))
bv[k,m] = Convert.ToDouble(al[5*kk+4]);
}
}
}
}
for(int k = 0; k < jx; k++)
{
for(int m = 0; m < jy; m++)
{
solution[k,m]=bv[k,m];
}
}
sor(a,b,c,d,e,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[,] a, double[,] b, double[,] c, double[,] d, double[,] e,
double[,] u, double rjac)
{
int jmax = a.GetLength(0);
int ipass,j,jsw,l,lsw,n;
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