using System;
using System.Collections.Generic;
namespace Science.Mathematics.PartialDifferentialEquation
{
///
/// LinearEllipticEquationTwoDimension
/// C_xxU_xx+C_xU_x+C_yyU_yy+C_yU_y+C_xyU_xy+C_propU=C_source
///
public class LinearEllipticEquation2D
{
public delegate double Coefficient(double x, double y);
public LinearEllipticEquation2D()
{
}
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 Coefficient cxx, cx, cyy, cy, cxy, 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 CoefficientOfDerivativeOfXandY
{
set{cxy=value;}
}
public Coefficient CoefficientOfProportional
{
set{cprop=value;}
}
public Coefficient CoefficientOfSourceTerm
{
set{csource=value;}
}
private List al = new List();
public void SetBoundaryValue(double xFrom, double yFrom,
double xTo, double yTo, double boundaryValue)
{
al.Add(xFrom);
al.Add(xTo);
al.Add(yFrom);
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];
double[,] f = new double[jx,jy];
double[,] g = 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]=cxx(x,y)+cx(x,y)*delta/2.0;
b[k,m]=cxx(x,y)-cx(x,y)*delta/2.0;
c[k,m]=cyy(x,y)+cy(x,y)*delta/2.0;
d[k,m]=cyy(x,y)-cy(x,y)*delta/2.0;
e[k,m]=-2.0*cxx(x,y)-2.0*cyy(x,y)+cprop(x,y)*delta*delta;
f[k,m]=csource(x,y)*delta*delta;
g[k,m]=cxy(x,y)/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+1])
&& Convert.ToDouble(al[5*kk+2]) < 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,f,g,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[,] f, double[,] g, double[,] u, double rjac)
{
int jmax = a.GetLength(0);
int ipass,j,jsw,l,lsw,n;
double anorm, omega=1.0,resid;
// double anormf=0.0;
// for (j=1;j