XIAN-FEM-2026June/opticsfem-master/kernel/Assemble_Scatter_2D_Boundary.cpp

#include"Assemble_Base.h"
#include"../function/BF.h"
#include"../function/Gauss.h"
#include"../common/define.h"
#include "../parser/mpParser.h"
#include"../Eigen/Sparse"
#include"../Eigen/Dense"

#include<cmath>

using namespace Eigen;

void OpticsFEM_2D_Scatter::Assemble_SBC()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();

	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(ONEDIM, LINE, BF_LINEFUNC * 2);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(ONEDIM, LINE, u, v, w, wght);

	//init of geo
	Vector3d* vertex = new Vector3d[5];
	Matrix3d Jac, InvJac;
	Vector3d normal;

	//init of basis function
	BF BF_Lagrange, BF_Nedelec;
	Vector3d* Et, ** Ez;
	BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	Et = new Vector3d[NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Ez[i] = new Vector3d[2];
	}

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	Vector3cd* Einc;
	complex<double>* EincPhi;
	Einc = new Vector3cd[NbrGuassPoints];
	EincPhi = new complex<double>[NbrGuassPoints];

	//overloop of SBC
	int NbrSBC = _mPhy->GetNbrSBC();
	for (int n = 0; n < NbrSBC; n++)
	{
		//SBC data
		int SBCDomain = _mPhy->GetSBCDomain(n) + 1;
		int SBCType = _mPhy->GetSBCType(n);
		string Eincx, Eincy, Eincz;
		Vector3d kdir;
		if (SBCType == 1)
		{
			_mPhy->GetSBDir(n, kdir);
			_mPhy->GetEinc(n, Eincx, Eincy, Eincz);
		}
		VectorXi edgesIndex;
		_mMesh->GetEdgesIndexOfDomain(SBCDomain, edgesIndex);

		//overloop of edgesOfSBC
		for (int m = 0; m < edgesIndex.rows(); m++)
		{
			//edge data
			Vector2i coonOfEdges;
			_mMesh->GetCoonOfEdges(edgesIndex(m), coonOfEdges);
			int numEle = coonOfEdges(0);
			int numEdge = coonOfEdges(1);

			//coordinate of vertex,integration coordinate, and mapping
			int MappingIndex[3];
			int index[3];
			for (int i = 0; i < 3; i++)
				_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
			if (numEdge == 0)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[1];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = 0.0;
				}

				index[0] = 1; index[1] = 2; index[2] = 1;
				MappingIndex[0] = _mMesh->GetTri(numEle, 0);
				MappingIndex[1] = _mMesh->GetTri(numEle, 1);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 0) + _mMesh->GetNbrVertex();
			}
			else if (numEdge == 1)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = 0.0;
					w[i] = (u[i] + 1.0) / 2.0;
				}

				index[0] = 1; index[1] = 3; index[2] = 2;
				MappingIndex[0] = _mMesh->GetTri(numEle, 0);
				MappingIndex[1] = _mMesh->GetTri(numEle, 2);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 1) + _mMesh->GetNbrVertex();
			}
			else if (numEdge == 2)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = (1.0 - u[i]) / 2.0;
				}

				index[0] = 2; index[1] = 3; index[2] = 3;
				MappingIndex[0] = _mMesh->GetTri(numEle, 1);
				MappingIndex[1] = _mMesh->GetTri(numEle, 2);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 2) + _mMesh->GetNbrVertex();
			}

			//Jac
			Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
			Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
			Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
			InvJac = Jac.inverse();

			//integration coe
			double phi1 = vertex[4](1) - vertex[3](1);
			double phi2 = vertex[3](0) - vertex[4](0);
			double integCoe = 0;
			if (vertex[3](0)!= vertex[4](0))
			{
				integCoe = phi1 / phi2;
				integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
			}
			else
			{
				integCoe = abs(phi1) / 2.0;
			}

			//normal
			_mMesh->GetNormOfEdges(_mMesh->GetDomainOfEdges(edgesIndex(m))-1, normal);

			//basis function
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				for (int j = 0; j < 2; j++)
				{
					BF_Lagrange.GetValueBF(index[j], v[i], w[i], 0.0, Ez[i][j]);
				}
				BF_Nedelec.GetValueBF(index[2], v[i], w[i], 0.0, Et[i]);
				Et[i] = InvJac * Et[i];
			}

			//material 
			double n = _mMatLib->GetN(_mMesh->GetDomainOfTri(numEle));
			complex<double> phySBC(0, k0 * n);

			//submatrix
			MatrixXcd Be;
			Be = MatrixXcd::Zero(2, 2);
			complex<double>Ae;
			for (int i = 0; i < 2; i++)
				for (int j = 0; j < 2; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Be(i, j) = Be(i, j) + wght[k] * integCoe * phySBC * (normal.cross(Ez[k][i])).dot(normal.cross(Ez[k][j]));
					}
			for (int k = 0; k < NbrGuassPoints; k++)
			{
				Ae = Ae + wght[k] * integCoe * phySBC * (normal.cross(Et[k])).dot(normal.cross(Et[k]));
			}

			//store in triple
			for (int i = 0; i < 2; i++)
				for (int j = 0; j < 2; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndex[i], MappingIndex[j], Be(i, j)));
				}
			_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndex[2], MappingIndex[2], Ae));

			//inc
			if (SBCType == 1)
			{
				//Einc
				for (int i = 0; i < NbrGuassPoints; i++)
				{
					double tempX = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
					double tempY = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);
					xx = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
					yy = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);

					p.SetExpr(_T(Eincx));
					Einc[i](0).real(p.Eval().GetFloat());
					Einc[i](0).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eincy));
					Einc[i](1).real(p.Eval().GetFloat());
					Einc[i](1).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eincz));
					Einc[i](2).real(p.Eval().GetFloat());
					Einc[i](2).imag(p.Eval().GetImag());

					complex<double> phi(0, -k0 * n * (kdir(0) * tempX + kdir(1) * tempY));
					EincPhi[i] = exp(
						phi);
				}

				//submatrix
				complex<double> be[3];
				for (int k = 0; k < NbrGuassPoints; k++)
				{
					be[0] = be[0] + wght[k] * integCoe * phySBC * /* EincPhi[k]* */ Ez[k][0].dot(normal.cross(Einc[k].cross(normal - kdir)));
					be[1] = be[1] + wght[k] * integCoe * phySBC * /* EincPhi[k]* */ Ez[k][1].dot(normal.cross(Einc[k].cross(normal - kdir)));
					be[2] = be[2] + wght[k] * integCoe * phySBC * /* EincPhi[k]* */Et[k].dot(normal.cross(Einc[k].cross(normal - kdir)));
				}

				//store in vector
				for (int i = 0; i < 3; i++)
					_mB_complex(MappingIndex[i]) = _mB_complex(MappingIndex[i]) + be[i];

			}
		}
	}

	delete[] u, v, w, wght;

	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] Ez[i];
	}
	delete[] Ez, Et, Einc, EincPhi;
}

void OpticsFEM_2D_Scatter::Assemble_BELE()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();

	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 2);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(TWODIM, TRIANGLE, u, v, w, wght);

	//init of geo
	Vector3d* vertex = new Vector3d[5];
	Matrix3d Jac, InvJac, TJac, JacS;

	//init of basis function
	BF BF_Lagrange, BF_Curl_Lagrange, BF_Nedelec, BF_Curl_Nedelec;
	int sdof, vdof;
	Vector3d** Et, ** curlEt, ** Ez, ** curlEz;
	/*Vector3d** bEt, ** bcurlEt, ** bEz, ** bcurlEz;*/
	sdof = BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	BF_Curl_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_LAGRANGE, BF_LINEFUNC);
	vdof = BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	BF_Curl_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
	Et = new Vector3d * [NbrGuassPoints];
	curlEt = new Vector3d * [NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	curlEz = new Vector3d * [NbrGuassPoints];

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Et[i] = new Vector3d[vdof];
		curlEt[i] = new Vector3d[vdof];
		Ez[i] = new Vector3d[sdof];
		curlEz[i] = new Vector3d[sdof];
	}

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	Vector3cd* bEt, * bEz, * bcurlEt, * bcurlEz;
	bEt = new Vector3cd[NbrGuassPoints];
	bEz = new Vector3cd[NbrGuassPoints];
	bcurlEt = new Vector3cd[NbrGuassPoints];
	bcurlEz = new Vector3cd[NbrGuassPoints];

	string Ebx, Eby, Ebz;
	_mPhy->GetEb(Ebx, Eby, Ebz);

	string curlEbx, curlEby, curlEbz;
	_mPhy->GetCurlEb(curlEbx, curlEby, curlEbz);

	//loop over tri
	int NbrTri = _mMesh->GetNbrTri();
	for (int n = 0; n < NbrTri; n++)
	{
		int domain = _mMesh->GetDomainOfTri(n);
		Eigen::VectorXi BELE;
		_mPhy->GetBELE(BELE);
		if (std::find(BELE.begin(), BELE.end(), domain) == BELE.end())
		{
			continue;
		}

		//coordinate of vertex
		for (int i = 0; i < 3; i++)
			_mMesh->GetVertex(_mMesh->GetTri(n, i), vertex[i]);

		//Jac
		Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
		Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
		Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
		InvJac = Jac.inverse();
		JacS(0, 0) = InvJac(1, 1);  JacS(0, 1) = -InvJac(1, 0); JacS(0, 2) = 0.;
		JacS(1, 0) = -InvJac(0, 1); JacS(1, 1) = InvJac(0, 0);  JacS(1, 2) = 0.;
		JacS(2, 0) = 0.;			JacS(2, 1) = 0.;		    JacS(2, 2) = 1.;
		double DetJac = fabs(Jac.determinant());
		TJac = Jac.transpose() / Jac.determinant();

		//basis function
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			for (int j = 0; j < sdof; j++)
			{
				BF_Lagrange.GetValueBF(j + 1, u[i], v[i], w[i], Ez[i][j]);
				BF_Curl_Lagrange.GetValueBF(j + 1, u[i], v[i], w[i], curlEz[i][j]);
				curlEz[i][j] = JacS * curlEz[i][j];
			}
			for (int j = 0; j < vdof; j++)
			{
				BF_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], Et[i][j]);
				Et[i][j] = InvJac * Et[i][j];
				BF_Curl_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], curlEt[i][j]);
				curlEt[i][j] = TJac * curlEt[i][j];
			}
		}

		//background field
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			xx = vertex[0][0] + Jac(0, 0) * u[i] + Jac(1, 0) * v[i];
			yy = vertex[0][1] + Jac(0, 0) * u[i] + Jac(1, 1) * v[i];

			p.SetExpr(_T(Ebx));
			bEt[i](0).real(p.Eval().GetFloat());
			bEt[i](0).imag(p.Eval().GetImag());
			p.SetExpr(_T(Eby));
			bEt[i](1).real(p.Eval().GetFloat());
			bEt[i](1).imag(p.Eval().GetImag());
			p.SetExpr(_T(Ebz));
			bEz[i](2).real(p.Eval().GetFloat());
			bEz[i](2).imag(p.Eval().GetImag());

			p.SetExpr(_T(curlEbx));
			bcurlEz[i](0).real(p.Eval().GetFloat());
			bcurlEz[i](0).imag(p.Eval().GetImag());
			p.SetExpr(_T(curlEby));
			bcurlEz[i](1).real(p.Eval().GetFloat());
			bcurlEz[i](1).imag(p.Eval().GetImag());
			p.SetExpr(_T(curlEbz));
			bcurlEt[i](2).real(p.Eval().GetFloat());
			bcurlEt[i](2).imag(p.Eval().GetImag());

		}

		//material
		Matrix3cd epsr = _mMatLib->GetEpsr(domain);
		Matrix3cd Mur = _mMatLib->GetMur(domain);
		Matrix3cd invMur;
		invMur = Mur.inverse();

		//submatrix
		MatrixXcd St, Sz, Tt, Tz;
		St = MatrixXcd::Zero(vdof, vdof);
		Sz = MatrixXcd::Zero(sdof, sdof);
		Tt = MatrixXcd::Zero(vdof, vdof);
		Tz = MatrixXcd::Zero(sdof, sdof);

		for (int i = 0; i < sdof; i++)
			for (int k = 0; k < NbrGuassPoints; k++)
			{
				Sz(i, 0) = Sz(i, 0) + wght[k] * DetJac * curlEz[k][i].dot(invMur * bcurlEz[k]);
				Tz(i, 0) = Tz(i, 0) + wght[k] * DetJac * k0 * k0 * Ez[k][i].dot(epsr * bEz[k]);
			}
		for (int i = 0; i < vdof; i++)
			for (int k = 0; k < NbrGuassPoints; k++)
			{
				St(i, 0) = St(i, 0) + wght[k] * DetJac * curlEt[k][i].dot(invMur * bcurlEt[k]);
				Tt(i, 0) = Tt(i, 0) + wght[k] * DetJac * k0 * k0 * Et[k][i].dot(epsr * bEt[k]);
			}

		//Mapping
		VectorXi MappingIndexS = VectorXi::Zero(sdof);
		VectorXi MappingIndexV = VectorXi::Zero(vdof);

		for (int i = 0; i < 3; i++)
		{
			MappingIndexS(i) = _mMesh->GetTri(n, i);
			MappingIndexV(i) = _mMesh->GetEdgeOfTri(n, i) + _mMesh->GetNbrVertex();
			_mB_complex(MappingIndexS[i]) = _mB_complex(MappingIndexS[i]) + Tz(i, 0) - Sz(i, 0);
			_mB_complex(MappingIndexV[i]) = _mB_complex(MappingIndexV[i]) + Tt(i, 0) - St(i, 0);
		}
	}

	//Line
	//init of Gauss point
	NbrGuassPoints = gauss.GetNbrGaussPoints(ONEDIM, LINE, BF_LINEFUNC * 4);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(ONEDIM, LINE, u, v, w, wght);

	Eigen::VectorXi BELE;
	_mPhy->GetBELE(BELE);
	int NbrBELE = BELE.rows();
	//find PML index
	int NbrPML = _mPhy->GetNbrPML();
	Eigen::VectorXi PMLEdge = Eigen::VectorXi::Zero(NbrPML);
	for (int i = 0; i < NbrPML; i++)
		PMLEdge(i) = _mPhy->GetPMLEdge(i) + 1;
	Eigen::VectorXi edgeIndexOfPML, edgeNumOfPML;//edges<65><73><EFBFBD><EFBFBD> <20><>ӦELE<4C><45><EFBFBD><EFBFBD>
	_mMesh->GetIndexOfDomain2(PMLEdge, edgeIndexOfPML, edgeNumOfPML);

	//overloop of edgesOfBELE
	for (int m = 0; m < edgeNumOfPML.rows(); m++)
	{
		//edge data
		Vector2i coonOfEdges;
		_mMesh->GetCoonOfEdges(edgeIndexOfPML(m), coonOfEdges);
		int numEle = coonOfEdges(0);
		int numEdge = coonOfEdges(1);

		//coordinate of vertex,integration coordinate, and mapping
		int index[3];
		for (int i = 0; i < 3; i++)
			_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
		if (numEdge == 0)
		{
			vertex[3] = vertex[0];
			vertex[4] = vertex[1];

			for (int i = 0; i < NbrGuassPoints; i++)
			{
				v[i] = (u[i] + 1.0) / 2.0;
				w[i] = 0.0;
			}
		}
		else if (numEdge == 1)
		{
			vertex[3] = vertex[0];
			vertex[4] = vertex[2];

			for (int i = 0; i < NbrGuassPoints; i++)
			{
				v[i] = 0.0;
				w[i] = (u[i] + 1.0) / 2.0;
			}
		}
		else if (numEdge == 2)
		{
			vertex[3] = vertex[1];
			vertex[4] = vertex[2];

			for (int i = 0; i < NbrGuassPoints; i++)
			{
				v[i] = (u[i] + 1.0) / 2.0;
				w[i] = (1.0 - u[i]) / 2.0;
			}
		}
		//Jac
		Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
		Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
		Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
		InvJac = Jac.inverse();
		//integration coe
		double phi1 = vertex[4](1) - vertex[3](1);
		double phi2 = vertex[3](0) - vertex[4](0);
		double integCoe = 0;
		if (vertex[3](0) != vertex[4](0))
		{
			integCoe = phi1 / phi2;
			integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
		}
		else
		{
			integCoe = abs(phi1) / 2.0;
		}

		//normal
		Vector3d normal;
		int edge = _mMesh->GetDomainOfEdges(edgeNumOfPML(m));
		_mMesh->GetNormOfEdges(_mMesh->GetDomainOfEdges(edgeNumOfPML(m)) - 1, normal);


		if (normal.isZero())
		{
			normal[0] = (vertex[3](0) + vertex[4](0)) / 2 - normal[0];
			normal[1] = (vertex[3](1) + vertex[4](1)) / 2 - normal[1];
			normal[2] = 0;
			normal.normalize();
		}

		//basis function
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			for (int j = 0; j < sdof; j++)
			{
				BF_Lagrange.GetValueBF(j + 1, v[i], w[i], 0.0, Ez[i][j]);

			}
			for (int j = 0; j < vdof; j++)
			{
				BF_Nedelec.GetValueBF(j + 1, v[i], w[i], 0.0, Et[i][j]);
				Et[i][j] = InvJac * Et[i][j];
			}
		}
		//back ground field
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			xx = vertex[0][0] + Jac(0, 0) * v[i] + Jac(1, 0) * w[i];
			yy = vertex[0][1] + Jac(0, 0) * v[i] + Jac(1, 1) * w[i];

			p.SetExpr(_T(curlEbx));
			bcurlEz[i](0).real(p.Eval().GetFloat());
			bcurlEz[i](0).imag(p.Eval().GetImag());
			p.SetExpr(_T(curlEby));
			bcurlEz[i](1).real(p.Eval().GetFloat());
			bcurlEz[i](1).imag(p.Eval().GetImag());
			p.SetExpr(_T(curlEbz));
			bcurlEt[i](2).real(p.Eval().GetFloat());
			bcurlEt[i](2).imag(p.Eval().GetImag());
		}
		//submatrix
		MatrixXcd bz, bt;
		bz = MatrixXcd::Zero(3, 1);
		bt = MatrixXcd::Zero(3, 1);
		for (int i = 0; i < 3; i++)
		{
			for (int j = 0; j < NbrGuassPoints; j++)
			{
				bz(i) = bz(i) - wght[j] * integCoe * Ez[j][i].dot(normal.cross(bcurlEz[j]));
				bt(i) = bt(i) - wght[j] * integCoe * Et[j][i].dot(normal.cross(bcurlEt[j]));
			}

		}
		//Mapping
		VectorXi MappingIndexS = VectorXi::Zero(sdof);
		VectorXi MappingIndexV = VectorXi::Zero(vdof);
		for (int i = 0; i < 3; i++)
		{
			MappingIndexS(i) = _mMesh->GetTri(numEle, i);
			MappingIndexV(i) = _mMesh->GetEdgeOfTri(numEle, i) + _mMesh->GetNbrVertex();
			_mB_complex(MappingIndexS[i]) = _mB_complex(MappingIndexS[i]) + bz(i);
			_mB_complex(MappingIndexV[i]) = _mB_complex(MappingIndexV[i]) + bt(i);
		}
	}

	delete[] u, v, w, wght;
	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] Et[i], curlEt[i], Ez[i], curlEz[i];
	}
	delete[] Et, curlEt, Ez, curlEz, bEt, bcurlEt, bEz, bcurlEz;
}

void OpticsFEM_2D_Scatter::Assemble_MAG()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();

	//Line
	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(ONEDIM, LINE, BF_LINEFUNC * 4);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(ONEDIM, LINE, u, v, w, wght);

	//init of basis function
	Vector3d* vertex = new Vector3d[5];
	Matrix3d Jac, InvJac, TJac, JacS;

	//init of basis function
	BF BF_Lagrange, BF_Curl_Lagrange, BF_Nedelec, BF_Curl_Nedelec;
	Vector3d** Et, ** Ez;
	int sdof, vdof;
	Et = new Vector3d * [NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	sdof = BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	vdof = BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Et[i] = new Vector3d[vdof];
		Ez[i] = new Vector3d[sdof];
	}

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	Vector3cd * bcurlEt, * bcurlEz;
	bcurlEt = new Vector3cd[NbrGuassPoints];
	bcurlEz = new Vector3cd[NbrGuassPoints];
	
	//find MAG index
	int NbrMAG = _mPhy->GetNbrMAG();
	string H0x, H0y, H0z;
	Eigen::VectorXi MAGDomain = Eigen::VectorXi::Zero(1);
	for (int n = 0; n < NbrMAG; n++)
	{
		_mPhy->GetH0(H0x, H0y, H0z, n);
		MAGDomain(0) = _mPhy->GetMAGDomain(n) + 1;
		Eigen::VectorXi edgeIndexOfMAG, edgeNumOfMAG;//edges<65><73><EFBFBD><EFBFBD> <20><>ӦMAG<41><47><EFBFBD><EFBFBD>
		_mMesh->GetIndexOfDomain2(MAGDomain, edgeIndexOfMAG, edgeNumOfMAG);

		//overloop of edgesOfMAG
		for (int m = 0; m < edgeNumOfMAG.rows(); m++)
		{
			//edge data
			Vector2i coonOfEdges;
			_mMesh->GetCoonOfEdges(edgeIndexOfMAG(m), coonOfEdges);
			int numEle = coonOfEdges(0);
			int numEdge = coonOfEdges(1);

			//coordinate of vertex,integration coordinate, and mapping
			int index[3];
			for (int i = 0; i < 3; i++)
				_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
			if (numEdge == 0)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[1];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = 0.0;
				}
			}
			else if (numEdge == 1)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = 0.0;
					w[i] = (u[i] + 1.0) / 2.0;
				}
			}
			else if (numEdge == 2)
			{
				vertex[3] = vertex[1];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = (1.0 - u[i]) / 2.0;
				}
			}
			//Jac
			Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
			Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
			Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
			InvJac = Jac.inverse();
			//integration coe
			double phi1 = vertex[4](1) - vertex[3](1);
			double phi2 = vertex[3](0) - vertex[4](0);
			double integCoe = 0;
			if (vertex[3](0) != vertex[4](0))
			{
				integCoe = phi1 / phi2;
				integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
			}
			else
			{
				integCoe = abs(phi1) / 2.0;
			}

			//normal
			Vector3d normal;
			int edge = _mMesh->GetDomainOfEdges(edgeNumOfMAG(m));
			_mMesh->GetNormOfEdges(_mMesh->GetDomainOfEdges(edgeNumOfMAG(m)) - 1, normal);

			if (normal.isZero())
			{
				normal[0] = (vertex[3](0) + vertex[4](0)) / 2 - normal[0];
				normal[1] = (vertex[3](1) + vertex[4](1)) / 2 - normal[1];
				normal[2] = 0;
				normal.normalize();
			}

			//basis function
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				for (int j = 0; j < sdof; j++)
				{
					BF_Lagrange.GetValueBF(j + 1, v[i], w[i], 0.0, Ez[i][j]);

				}
				for (int j = 0; j < vdof; j++)
				{
					BF_Nedelec.GetValueBF(j + 1, v[i], w[i], 0.0, Et[i][j]);
					Et[i][j] = InvJac * Et[i][j];
				}
			}
			//back ground field
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				xx = vertex[0][0] + Jac(0, 0) * v[i] + Jac(1, 0) * w[i];
				yy = vertex[0][1] + Jac(0, 1) * v[i] + Jac(1, 1) * w[i];

				p.SetExpr(_T(H0x));
				bcurlEz[i](0).real(p.Eval().GetFloat());
				bcurlEz[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(H0y));
				bcurlEz[i](1).real(p.Eval().GetFloat());
				bcurlEz[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(H0z));
				bcurlEt[i](2).real(p.Eval().GetFloat());
				bcurlEt[i](2).imag(p.Eval().GetImag());
			}
			//submatrix
			MatrixXcd bz, bt;
			bz = MatrixXcd::Zero(3, 1);
			bt = MatrixXcd::Zero(3, 1);
			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < NbrGuassPoints; j++)
				{
					bz(i) = bz(i) - complex<double>(0.0, 1.0) * k0 * wght[j] * integCoe * Ez[j][i].dot(normal.cross(bcurlEz[j]));
					bt(i) = bt(i) - complex<double>(0.0, 1.0) * k0 * wght[j] * integCoe * Et[j][i].dot(normal.cross(bcurlEt[j]));
				}

			}
			//Mapping
			VectorXi MappingIndexS = VectorXi::Zero(sdof);
			VectorXi MappingIndexV = VectorXi::Zero(vdof);
			for (int i = 0; i < 3; i++)
			{
				MappingIndexS(i) = _mMesh->GetTri(numEle, i);
				MappingIndexV(i) = _mMesh->GetEdgeOfTri(numEle, i) + _mMesh->GetNbrVertex();
				_mB_complex(MappingIndexS[i]) = _mB_complex(MappingIndexS[i]) + bz(i);
				_mB_complex(MappingIndexV[i]) = _mB_complex(MappingIndexV[i]) + bt(i);
			}
		}
	}

	delete[] u, v, w, wght;
	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] Et[i], Ez[i];
	}
	delete[] Et, Ez, bcurlEt, bcurlEz;
}

void OpticsFEM_2D_Scatter::Assemble_SCD()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();

	//Line
	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(ONEDIM, LINE, BF_LINEFUNC * 4);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(ONEDIM, LINE, u, v, w, wght);

	//init of basis function
	Vector3d* vertex = new Vector3d[5];
	Matrix3d Jac, InvJac, TJac, JacS;

	//init of basis function
	BF BF_Lagrange, BF_Curl_Lagrange, BF_Nedelec, BF_Curl_Nedelec;
	Vector3d** Et, ** Ez;
	int sdof, vdof;
	Et = new Vector3d * [NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	sdof = BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	vdof = BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Et[i] = new Vector3d[vdof];
		Ez[i] = new Vector3d[sdof];
	}

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	Vector3cd* bEt, * bEz;
	bEt = new Vector3cd[NbrGuassPoints];
	bEz = new Vector3cd[NbrGuassPoints];

	//find SCD index
	int NbrSCD = _mPhy->GetNbrSCD();
	string J0x, J0y, J0z;
	Eigen::VectorXi SCDDomain = Eigen::VectorXi::Zero(1);
	for (int n = 0; n < NbrSCD; n++)
	{
		_mPhy->GetJ0(J0x, J0y, J0z, n);
		SCDDomain(0) = _mPhy->GetSCDDomain(n) + 1;
		Eigen::VectorXi edgeIndexOfSCD, edgeNumOfSCD;//edges<65><73><EFBFBD><EFBFBD> <20><>ӦSCD<43><44><EFBFBD><EFBFBD>
		_mMesh->GetIndexOfDomain2(SCDDomain, edgeIndexOfSCD, edgeNumOfSCD);
		//overloop of edgesOfSCD
		for (int m = 0; m < edgeNumOfSCD.rows(); m++)
		{
			//edge data
			Vector2i coonOfEdges;
			_mMesh->GetCoonOfEdges(edgeIndexOfSCD(m), coonOfEdges);
			int numEle = coonOfEdges(0);
			int numEdge = coonOfEdges(1);
			int index[3];
			for (int i = 0; i < 3; i++)
				_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
			if (numEdge == 0)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[1];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = 0.0;
				}
			}
			else if (numEdge == 1)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = 0.0;
					w[i] = (u[i] + 1.0) / 2.0;
				}
			}
			else if (numEdge == 2)
			{
				vertex[3] = vertex[1];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = (1.0 - u[i]) / 2.0;
				}
			}
			//Jac
			Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
			Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
			Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
			InvJac = Jac.inverse();
			//integration coe
			double phi1 = vertex[4](1) - vertex[3](1);
			double phi2 = vertex[3](0) - vertex[4](0);
			double integCoe = 0;
			if (vertex[3](0) != vertex[4](0))
			{
				integCoe = phi1 / phi2;
				integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
			}
			else
			{
				integCoe = abs(phi1) / 2.0;
			}

			//basis function
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				for (int j = 0; j < sdof; j++)
				{
					BF_Lagrange.GetValueBF(j + 1, v[i], w[i], 0.0, Ez[i][j]);

				}
				for (int j = 0; j < vdof; j++)
				{
					BF_Nedelec.GetValueBF(j + 1, v[i], w[i], 0.0, Et[i][j]);
					Et[i][j] = InvJac * Et[i][j];
				}
			}
			//back ground field
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				xx = vertex[0][0] + Jac(0, 0) * v[i] + Jac(1, 0) * w[i];
				yy = vertex[0][1] + Jac(0, 1) * v[i] + Jac(1, 1) * w[i];

				p.SetExpr(_T(J0x));
				bEt[i](0).real(p.Eval().GetFloat());
				bEt[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(J0y));
				bEt[i](1).real(p.Eval().GetFloat());
				bEt[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(J0z));
				bEz[i](2).real(p.Eval().GetFloat());
				bEz[i](2).imag(p.Eval().GetImag());
			}
			//submatrix
			MatrixXcd bz, bt;
			bz = MatrixXcd::Zero(3, 1);
			bt = MatrixXcd::Zero(3, 1);
			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; j < NbrGuassPoints; j++)
				{
					bz(i) = bz(i) - complex<double>(0.0, 1.0) * sqrt(mu0 / epsilon0) * k0 * wght[j] * integCoe * Ez[j][i].dot(bEz[j]);
					bt(i) = bt(i) - complex<double>(0.0, 1.0) * sqrt(mu0 / epsilon0) * k0 * wght[j] * integCoe * Et[j][i].dot(bEt[j]);
				}

			}
			//Mapping
			VectorXi MappingIndexS = VectorXi::Zero(sdof);
			VectorXi MappingIndexV = VectorXi::Zero(vdof);
			for (int i = 0; i < 3; i++)
			{
				MappingIndexS(i) = _mMesh->GetTri(numEle, i);
				MappingIndexV(i) = _mMesh->GetEdgeOfTri(numEle, i) + _mMesh->GetNbrVertex();
				_mB_complex(MappingIndexS[i]) = _mB_complex(MappingIndexS[i]) + bz(i);
				_mB_complex(MappingIndexV[i]) = _mB_complex(MappingIndexV[i]) + bt(i);
			}
		}
	}
	delete[] u, v, w, wght;
	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] Et[i], Ez[i];
	}
	delete[] Et, Ez, bEt, bEz;
}

void OpticsFEM_2D_Scatter::Assemble_MPD()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();
	Eigen::MatrixXd  m;
	Eigen::VectorXi MPD;
	_mPhy->GetMPD(MPD);
	_mPhy->GetMPDData(m);
	
	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 4);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(TWODIM, TRIANGLE, u, v, w, wght);

	//init of geo
	Vector3d* vertex = new Vector3d[3];
	Matrix3d Jac, InvJac, TJac, JacS;

	double xx, yy;
	Vector3cd* bEt, * bEz;
	bEt = new Vector3cd[NbrGuassPoints];
	bEz = new Vector3cd[NbrGuassPoints];

	
	//init of basis function
	BF BF_Lagrange, BF_Curl_Lagrange, BF_Nedelec, BF_Curl_Nedelec;
	int sdof, vdof;
	Vector3d** curlEt, ** curlEz;
	sdof = BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	BF_Curl_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_LAGRANGE, BF_LINEFUNC);
	vdof = BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	BF_Curl_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
	curlEt = new Vector3d * [NbrGuassPoints];
	curlEz = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		curlEt[i] = new Vector3d[vdof];
		curlEz[i] = new Vector3d[sdof];
	}
	//loop over tri
	int NbrTri = _mMesh->GetNbrTri();
	int NbrMPD = MPD.size();
	int num = 0;
	for (int i = 0; i < NbrMPD; i++)
	{
		for (int n = 0; n < NbrTri; n++)
		{
			bool  found = false;
			for (int j = 0; j < 3; j++)
			{

				if (MPD[i] == _mMesh->GetTri(n, j))
				{
					found = true;
					break;
				}
			}

			if (found) {
				//coordinate of vertex
				num = num + 1;
				for (int j = 0; j < 3; j++)
					_mMesh->GetVertex(_mMesh->GetTri(n, j), vertex[j]);

				//Jac
				Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
				Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
				Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
				InvJac = Jac.inverse();
				JacS(0, 0) = InvJac(1, 1);  JacS(0, 1) = -InvJac(1, 0); JacS(0, 2) = 0.;
				JacS(1, 0) = -InvJac(0, 1); JacS(1, 1) = InvJac(0, 0);  JacS(1, 2) = 0.;
				JacS(2, 0) = 0.;			JacS(2, 1) = 0.;		    JacS(2, 2) = 1.;
				double DetJac = fabs(Jac.determinant());
				TJac = Jac.transpose() / Jac.determinant();

				//basis function
				for (int j = 0; j < NbrGuassPoints; j++)
				{
					for (int k = 0; k < sdof; k++)
					{
						BF_Curl_Lagrange.GetValueBF(k + 1, u[j], v[j], w[j], curlEz[j][k]);
						curlEz[j][k] = JacS * curlEz[j][k];
					}
					for (int k = 0; k < vdof; k++)
					{
						BF_Curl_Nedelec.GetValueBF(k + 1, u[j], v[j], w[j], curlEt[j][k]);
						curlEt[j][k] = TJac * curlEt[j][k];
					}
				}

				//submatrix
				Eigen::Vector3d m_current = m.row(i);
				MatrixXcd Tt, Tz;
				Tt = MatrixXcd::Zero(vdof, 1);
				Tz = MatrixXcd::Zero(sdof, 1);
				for (int j = 0; j < sdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Tz(j) = Tz(j) + complex<double>(0.0, -1.0) * sqrt(mu0 / epsilon0) * k0 * wght[k] * DetJac * curlEz[k][j].dot(m_current);
					}
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Tt(j) = Tt(j) + complex<double>(0.0, -1.0) * sqrt(mu0 / epsilon0) * k0 * wght[k] * DetJac * curlEt[k][j].dot(m_current);
					}

				//Mapping
				VectorXi MappingIndexS = VectorXi::Zero(sdof);
				VectorXi MappingIndexV = VectorXi::Zero(vdof);
				for (int j = 0; j < 3; j++)
				{
					MappingIndexS(j) = _mMesh->GetTri(n, j);
					MappingIndexV(j) = _mMesh->GetEdgeOfTri(n, j) + _mMesh->GetNbrVertex();
					_mB_complex(MappingIndexS[j]) = _mB_complex(MappingIndexS[j]) + Tz(j);
					_mB_complex(MappingIndexV[j]) = _mB_complex(MappingIndexV[j]) + Tt(j);
				}
			}
		}		
	}
	delete[] u, v, w, wght;
	delete[] vertex;		
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] curlEt[i], curlEz[i];
	}
	delete[] bEt, bEz;
	
}

void OpticsFEM_2D_Scatter::Assemble_EPD()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();
	Eigen::MatrixXd  p;
	Eigen::VectorXi EPD;
	_mPhy->GetEPD(EPD);
	_mPhy->GetEPDData(p);

	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	double* u, * v, * w, * wght;
	NbrGuassPoints = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 4);
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(TWODIM, TRIANGLE, u, v, w, wght);

	//init of geo
	Vector3d* vertex = new Vector3d[3];
	Matrix3d Jac, InvJac, TJac, JacS;

	double xx, yy;
	Vector3cd* bEt, * bEz;
	bEt = new Vector3cd[NbrGuassPoints];
	bEz = new Vector3cd[NbrGuassPoints];


	//init of basis function
	BF BF_Lagrange, BF_Curl_Lagrange, BF_Nedelec, BF_Curl_Nedelec;
	int sdof, vdof;
	Vector3d** Et, ** Ez;
	sdof = BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	BF_Curl_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_LAGRANGE, BF_LINEFUNC);
	vdof = BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	BF_Curl_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
	Et = new Vector3d * [NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Et[i] = new Vector3d[vdof];
		Ez[i] = new Vector3d[sdof];
	}
	//loop over tri
	int NbrTri = _mMesh->GetNbrTri();
	int NbrEPD = EPD.size();
	int num = 0;
	for (int i = 0; i < NbrEPD; i++)
	{
		for (int n = 0; n < NbrTri; n++)
		{
			bool  found = false;
			for (int j = 0; j < 3; j++)
			{

				if (EPD[i] == _mMesh->GetTri(n, j))
				{
					found = true;
					break;
				}
			}

			if (found) {
				//coordinate of vertex
				num = num + 1;
				for (int j = 0; j < 3; j++)
					_mMesh->GetVertex(_mMesh->GetTri(n, j), vertex[j]);

				//Jac
				Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
				Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
				Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
				InvJac = Jac.inverse();
				JacS(0, 0) = InvJac(1, 1);  JacS(0, 1) = -InvJac(1, 0); JacS(0, 2) = 0.;
				JacS(1, 0) = -InvJac(0, 1); JacS(1, 1) = InvJac(0, 0);  JacS(1, 2) = 0.;
				JacS(2, 0) = 0.;			JacS(2, 1) = 0.;		    JacS(2, 2) = 1.;
				double DetJac = fabs(Jac.determinant());
				TJac = Jac.transpose() / Jac.determinant();

				//basis function
				for (int j = 0; j < NbrGuassPoints; j++)
				{
					for (int k = 0; k < sdof; k++)
					{
						BF_Lagrange.GetValueBF(k + 1, u[j], v[j], w[j], Ez[j][k]);
					}
					for (int k = 0; k < vdof; k++)
					{
						BF_Nedelec.GetValueBF(k + 1, u[j], v[j], w[j], Et[j][k]);
						Et[j][k] = InvJac * Et[j][k];
					}
				}

				//submatrix
				Eigen::Vector3d p_current = p.row(i);
				MatrixXcd Tt, Tz;
				Tt = MatrixXcd::Zero(vdof, 1);
				Tz = MatrixXcd::Zero(sdof, 1);
				for (int j = 0; j < sdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Tz(j) = Tz(j) + complex<double>(0.0, -1.0) * sqrt(mu0/epsilon0) * k0 * wght[k] * DetJac * Ez[k][j].dot(p_current);
					}
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Tt(j) = Tt(j) + complex<double>(0.0, -1.0) * sqrt(mu0 / epsilon0) * k0 * wght[k] * DetJac * Et[k][j].dot(p_current);
					}
				//Mapping
				VectorXi MappingIndexS = VectorXi::Zero(sdof);
				VectorXi MappingIndexV = VectorXi::Zero(vdof);
				for (int j = 0; j < 3; j++)
				{
					MappingIndexS(j) = _mMesh->GetTri(n, j);
					MappingIndexV(j) = _mMesh->GetEdgeOfTri(n, j) + _mMesh->GetNbrVertex();
					_mB_complex(MappingIndexS[j]) = _mB_complex(MappingIndexS[j]) + Tz(j);
					_mB_complex(MappingIndexV[j]) = _mB_complex(MappingIndexV[j]) + Tt(j);
				}
			}
		}
	}
	delete[] u, v, w, wght;
	delete[] vertex;
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] Et[i], Ez[i];
	}
	delete[] bEt, bEz;

}

void OpticsFEM_2D_Scatter::Assemble_PEC_ELE()
{
	//find PEC index
	int NbrPEC = _mPhy->GetNbrPEC();
	Eigen::VectorXi PECDomain = Eigen::VectorXi::Zero(NbrPEC);
	for (int i = 0; i < NbrPEC; i++)
		PECDomain(i) = _mPhy->GetPECDomain(i) + 1;
	Eigen::VectorXi nodeIndexOfPEC, edgeIndexOfPEC;
	_mMesh->GetIndexOfDomain(PECDomain, nodeIndexOfPEC, edgeIndexOfPEC);

	//find ELE index
	int NbrELE = _mPhy->GetNbrElE();
	Eigen::VectorXi ELEDomain = Eigen::VectorXi::Zero(NbrELE);
	for (int i = 0; i < NbrELE; i++)
		ELEDomain(i) = _mPhy->GetELEDomain(i) + 1;
	Eigen::VectorXi edgeIndexOfELE, edgeNumOfELE;//edges<65><73><EFBFBD><EFBFBD> <20><>ӦELE<4C><45><EFBFBD><EFBFBD>
	_mMesh->GetIndexOfDomain2(ELEDomain, edgeIndexOfELE, edgeNumOfELE);

	
	//init BCIndex BCValue
	int NbrBC = edgeIndexOfPEC.rows() + nodeIndexOfPEC.rows() + edgeNumOfELE.rows() * 3; //<2F><><EFBFBD><EFBFBD>+һ<><D2BB><EFBFBD><EFBFBD>
	Eigen::VectorXi BCIndex = Eigen::VectorXi::Zero(NbrBC);// 1:3 ELE  4:6 PEC
	Eigen::VectorXcd BCValue = Eigen::VectorXcd::Zero(NbrBC);

	//init of geo
	Vector3d* vertex = new Vector3d[5];

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy, zz;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	p.DefineVar(_T("z"), mup::Variable(&zz));
	Vector3cd E0;

	//overloop edge of ELE
	int bcOffset = 0;
	int NbrVertex = _mMesh->GetNbrVertex();
	for (int m = 0; m < NbrELE; m++)
	{
		Eigen::VectorXi ELEDomain = Eigen::VectorXi::Zero(NbrELE);
		ELEDomain(0) = _mPhy->GetELEDomain(m) + 1;
		Eigen::VectorXi edgeIndexOfELE, currentEdgeNum;//edges<65><73><EFBFBD><EFBFBD> <20><>ӦELE<4C><45><EFBFBD><EFBFBD>
		_mMesh->GetIndexOfDomain2(ELEDomain, edgeIndexOfELE, currentEdgeNum);
		string E0x, E0y, E0z;
		_mPhy->GetE0(E0x, E0y, E0z,m);
		for (int n = 0; n <currentEdgeNum.rows(); n++)
		{
			//edge data
			Eigen::Vector2i coonOfEdges;
			_mMesh->GetCoonOfEdges(edgeIndexOfELE(n), coonOfEdges);
			int numEle = coonOfEdges(0);
			int numEdge = coonOfEdges(1);

			for (int i = 0; i < 3; i++)
				_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);

			Eigen::Vector3d ee = Eigen::Vector3d::Zero();
			
			if (numEdge == 0)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[1];
				ee(0) = vertex[1](0) - vertex[0](0);
				ee(1) = vertex[1](1) - vertex[0](1);

				BCIndex(n + bcOffset)= _mMesh->GetTri(numEle, 0);
				BCIndex(n + bcOffset + edgeNumOfELE.rows()) = _mMesh->GetTri(numEle, 1);
				BCIndex(n + bcOffset + edgeNumOfELE.rows() * 2) = _mMesh->GetEdgeOfTri(numEle, 0) + NbrVertex;
			}
			else if (numEdge == 1)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];
				ee(0) = vertex[2](0) - vertex[0](0);
				ee(1) = vertex[2](1) - vertex[0](1);

				BCIndex(n + bcOffset) = _mMesh->GetTri(numEle, 0);
				BCIndex(n + bcOffset + edgeNumOfELE.rows()) = _mMesh->GetTri(numEle, 2);
				BCIndex(n + bcOffset + edgeNumOfELE.rows() * 2) = _mMesh->GetEdgeOfTri(numEle, 1) + NbrVertex;
			}
			else if (numEdge == 2)
			{
				vertex[3] = vertex[1];
				vertex[4] = vertex[2];
				ee(0) = vertex[2](0) - vertex[1](0);
				ee(1) = vertex[2](1) - vertex[1](1);

				BCIndex(n + bcOffset) = _mMesh->GetTri(numEle, 1);
				BCIndex(n + bcOffset + edgeNumOfELE.rows()) = _mMesh->GetTri(numEle, 2);
				BCIndex(n + bcOffset + edgeNumOfELE.rows() * 2) = _mMesh->GetEdgeOfTri(numEle, 2) + NbrVertex;
			}
			xx = (vertex[3](0) + vertex[4](0)) / 2;
			yy = (vertex[3](1) + vertex[4](1)) / 2;
			if (_mPhy->GetNbrBELE())
			{
				string Ebx, Eby, Ebz;
				string ExFunc, EyFunc, EzFunc;
				_mPhy->GetEb(Ebx, Eby, Ebz);
				ExFunc = E0x + "-" + Ebx;
				EyFunc = E0y + "-" + Eby;
				EzFunc = E0z + "-" + Ebz;
				p.SetExpr(_T(ExFunc));
				E0(0).real(p.Eval().GetFloat());
				E0(0).imag(p.Eval().GetImag());
				p.SetExpr(_T(EyFunc));
				E0(1).real(p.Eval().GetFloat());
				E0(1).imag(p.Eval().GetImag());

				xx = vertex[3](0); yy = vertex[4](0);
				p.SetExpr(_T(EzFunc));
				E0(2).real(p.Eval().GetFloat());
				E0(2).imag(p.Eval().GetImag());
				BCValue(n + bcOffset) = E0(2);

				xx = vertex[4](0); yy = vertex[5](0);
				p.SetExpr(_T(EzFunc));
				E0(2).real(p.Eval().GetFloat());
				E0(2).imag(p.Eval().GetImag());
				BCValue(n + bcOffset + edgeNumOfELE.rows()) = E0(2);
			}
			else
			{
				p.SetExpr(_T(E0x));
				E0(0).real(p.Eval().GetFloat());
				E0(0).imag(p.Eval().GetImag());
				p.SetExpr(_T(E0y));
				E0(1).real(p.Eval().GetFloat());
				E0(1).imag(p.Eval().GetImag());

				xx = vertex[3](0); yy = vertex[4](0);
				p.SetExpr(_T(E0z));
				E0(2).real(p.Eval().GetFloat());
				E0(2).imag(p.Eval().GetImag());
				BCValue(n + bcOffset) = E0(2);

				xx = vertex[4](0); yy = vertex[5](0);
				p.SetExpr(_T(E0z));
				E0(2).real(p.Eval().GetFloat());
				E0(2).imag(p.Eval().GetImag());
				BCValue(n + bcOffset + edgeNumOfELE.rows()) = E0(2);
			}

			BCValue(n + bcOffset + edgeNumOfELE.rows() * 2) = E0.dot(ee);
		}
		bcOffset += currentEdgeNum.rows();
	}
	delete[] vertex;

	
	//store PEC data
	for (int i = 0; i < nodeIndexOfPEC.rows(); i++)
	{
		BCIndex(i + bcOffset * 3) = nodeIndexOfPEC(i);
		BCValue(i + bcOffset * 3) = 0.;
	}
	for (int i = 0; i < edgeIndexOfPEC.rows(); i++)
	{
		BCIndex(i + bcOffset * 3 + nodeIndexOfPEC.rows()) = edgeIndexOfPEC(i) + NbrVertex;
		BCValue(i + bcOffset * 3 + nodeIndexOfPEC.rows()) = 0.;
	}

	//<2F><><EFBFBD><EFBFBD>
	Eigen::VectorXi tempIndex = Eigen::VectorXi::Zero(NbrBC);
	for (int i = 0; i < NbrBC; i++)
		tempIndex(i) = i;
	QuickSort(BCIndex, tempIndex, 0, NbrBC - 1);
	
	//ȥ<><C8A5>
	Unique(BCIndex, tempIndex);
	Eigen::VectorXcd tempValue = Eigen::VectorXcd::Zero(BCIndex.rows());
	for (int i = 0; i < BCIndex.rows(); i++)
		tempValue(i) = BCValue(tempIndex(i));
	BCValue = tempValue;

	
	//<2F><><EFBFBD><EFBFBD>
	if (_mIsReal)
	{
		Eigen::VectorXd tempValue;
		for (int i = 0; i < BCValue.rows(); i++)
		{
			tempValue = Eigen::VectorXd::Zero(_mA_real.rows());
			tempValue(BCIndex(i)) = 1.;
			_mB_real = _mB_real - _mA_real * tempValue * BCValue(i).real();
		}
		
		Eigen::SparseMatrix<double, Eigen::RowMajor> P(_mA_real.rows(), _mA_real.rows());
		std::vector<Eigen::Triplet<double>> tempTriple;
		int num = 0;
		for (int i = 0; i < _mA_real.rows(); i++)
		{
			if (num < BCIndex.rows())
			{
				if (i == BCIndex(num))
					num++;
				else
					tempTriple.push_back(Eigen::Triplet<double>(i, i, 1.));
			}
			else
				tempTriple.push_back(Eigen::Triplet<double>(i, i, 1.));
		}
		P.setFromTriplets(tempTriple.begin(), tempTriple.end());
		_mA_real = P * _mA_real * P;
		
		P.setZero();
		tempTriple.clear();
		
		for (int i = 0; i < BCIndex.rows(); i++)
		{
			tempTriple.push_back(Eigen::Triplet<double>(BCIndex(i), BCIndex(i), 1.));
			_mB_real(BCIndex(i)) = BCValue(i).real();
		}
		P.setFromTriplets(tempTriple.begin(), tempTriple.end());
		_mA_real = _mA_real + P;
		
	}
	else
	{
		Eigen::VectorXcd tempValue;
		for (int i = 0; i < BCValue.rows(); i++)
		{
			tempValue = Eigen::VectorXcd::Zero(_mA_complex.rows());
			tempValue(BCIndex(i)) = 1;
			_mB_complex = _mB_complex - _mA_complex * tempValue * BCValue(i);
		}

		Eigen::SparseMatrix<complex<double>, Eigen::RowMajor> P(_mA_complex.rows(), _mA_complex.rows());
		std::vector<Eigen::Triplet<complex<double>>> tempTriple;
		int num = 0;
		for (int i = 0; i < _mA_complex.rows(); i++)
		{
			if (num < BCIndex.rows())
			{
				if (i == BCIndex(num))
					num++;
				else
					tempTriple.push_back(Eigen::Triplet<complex<double>>(i, i, 1.));
			}
			else
				tempTriple.push_back(Eigen::Triplet<complex<double>>(i, i, 1.));
		}
		P.setFromTriplets(tempTriple.begin(), tempTriple.end());
		_mA_complex = P * _mA_complex * P;

		P.setZero();
		tempTriple.clear();
		for (int i = 0; i < BCIndex.rows(); i++)
		{
			tempTriple.push_back(Eigen::Triplet<complex<double>>(BCIndex(i), BCIndex(i), 1.));
			_mB_complex(BCIndex(i)) = BCValue(i);
		}
		P.setFromTriplets(tempTriple.begin(), tempTriple.end());
		_mA_complex = _mA_complex + P;
	}
	
}

void OpticsFEM_2D_Scatter::Assemble_PBC()
{
	//find PBC index
	Eigen::VectorXi srcNodeIndex, dstNodeIndex, srcEdgeIndex, dstEdgeIndex;
	Eigen::VectorXcd nodePhi, edgePhi, PBCData;
	Eigen::MatrixXi PBC;
	PBC = Eigen::MatrixXi::Zero(_mPhy->GetNbrPBC(), 2);
	PBCData = Eigen::VectorXcd::Zero(_mPhy->GetNbrPBC());
	for (int i = 0; i < _mPhy->GetNbrPBC(); i++)
	{
		PBC(i, 0) = _mPhy->GetSrcDomain(i);
		PBC(i, 1) = _mPhy->GetDstDomain(i);
		PBCData(i) = _mPhy->GetPBCPhi(i);
	}
	_mMesh->GetIndexOfPBC(PBC, PBCData, srcNodeIndex, dstNodeIndex, nodePhi, srcEdgeIndex, dstEdgeIndex, edgePhi);

	//<2F><> 1:DOF
	//<2F><> 1<><31>(DOF-NbrPEC-NbrDst)
	int NbrVertex = _mMesh->GetNbrVertex();
	int NbrDstIndex = dstNodeIndex.rows() + dstEdgeIndex.rows();

	Eigen::VectorXi dstIndex = Eigen::VectorXi::Zero(NbrDstIndex);
	Eigen::VectorXi srcIndex = Eigen::VectorXi::Zero(NbrDstIndex);
	Eigen::VectorXcd IndexPhi = Eigen::VectorXcd::Zero(NbrDstIndex);

	for (int i = 0; i < dstNodeIndex.rows(); i++)
	{
		dstIndex(i) = dstNodeIndex(i);
		srcIndex(i) = srcNodeIndex(i);
		IndexPhi(i) = nodePhi(i);
	}
	for (int i = 0; i < dstEdgeIndex.rows(); i++)
	{
		int num = dstNodeIndex.rows();
		dstIndex(i + num) = dstEdgeIndex(i) + NbrVertex;
		srcIndex(i + num) = srcEdgeIndex(i) + NbrVertex;
		IndexPhi(i + num) = edgePhi(i);
	}
	//<2F><><EFBFBD><EFBFBD>
	Eigen::VectorXi tempIndex = Eigen::VectorXi::Zero(NbrDstIndex);
	for (int i = 0; i < NbrDstIndex; i++)
		tempIndex(i) = i;
	QuickSort(dstIndex, tempIndex, 0, NbrDstIndex - 1);

	//ȥ<><C8A5>
	Unique(dstIndex, tempIndex);
	Eigen::VectorXcd tempIndexPhi = Eigen::VectorXcd::Zero(dstIndex.rows());
	Eigen::VectorXi tempSrcIndex = Eigen::VectorXi::Zero(dstIndex.rows());
	for (int i = 0; i < dstIndex.rows(); i++)
	{
		tempIndexPhi(i) = IndexPhi(tempIndex(i));
		tempSrcIndex(i) = srcIndex(tempIndex(i));
	}
	IndexPhi = tempIndexPhi;
	srcIndex = tempSrcIndex;

	//<2F><>ͬ<EFBFBD><CDAC>Ŀ
	int dof = _mMesh->GetNbrEdge() + _mMesh->GetNbrVertex();
	int vdof = _mMesh->GetNbrEdge();
	int sdof = _mMesh->GetNbrVertex();

	//P<><50><EFBFBD><EFBFBD>
	int num = 0;
	if (_mIsReal)
	{
		std::vector<Eigen::Triplet<double>> temp_P;
		for (int i = 0; i < dof; i++) // 
		{
			temp_P.push_back(Eigen::Triplet<double>(i, i, 1));
			if (num < dstIndex.size()) //<2F><><EFBFBD><EFBFBD>PBC
			{
				if (i == dstIndex(num))
				{
					temp_P.push_back(Eigen::Triplet<double>(dstIndex(num), srcIndex(num), IndexPhi(num).real()));
					num++;
				}
			}
		}
		/*<2A><>ʱP<CAB1><50><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>*/
		std::sort(temp_P.begin(), temp_P.end(),
			[](const Eigen::Triplet<double>& a, const Eigen::Triplet<double>& b)
			{
				return a.col() < b.col();
			});


		/*ɾ<><C9BE>dstIndex<65><78>Ӧ<EFBFBD><D3A6>*/
		std::vector<Eigen::Triplet<double>> tripleP_real;
		num = 0;// <20><><EFBFBD>ڸ<EFBFBD><DAB8><EFBFBD><EFBFBD>¾<EFBFBD><C2BE><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
		int	DeleteIndex = 0;
		for (int i = 0; i < dof; i++)
		{
			DeleteIndex = 0;
			for (int j = 0; j < dstIndex.rows(); j++)
			{
				if (i == dstIndex[j])
				{
					DeleteIndex = 1;
				}
			}
			if (DeleteIndex != 1)
			{
				for (const auto& triplet : temp_P)
				{
					if (triplet.col() == i)
					{
						tripleP_real.push_back(Eigen::Triplet<double>(triplet.row(), triplet.col() - num, triplet.value()));
					}
				}
			}
			for (int j = 0; j < dstIndex.rows(); j++)
			{
				if (i == dstIndex[j])
				{
					num = num + 1;
				}
			}

		}
		_mP_real = Eigen::SparseMatrix<double, Eigen::RowMajor>(dof, dof - dstIndex.rows());
		_mP_real.setFromTriplets(tripleP_real.begin(), tripleP_real.end());

	}
	else
	{
		std::vector<Eigen::Triplet<complex<double>>> temp_P_complex;
		for (int i = 0; i < dof; i++) // 
		{
			temp_P_complex.push_back(Eigen::Triplet<complex<double>>(i, i, 1));
			if (num < dstIndex.size()) //<2F><><EFBFBD><EFBFBD>PEC<45><43>PBC
			{
				if (i == dstIndex(num))
				{
					temp_P_complex.push_back(Eigen::Triplet<complex<double>>(dstIndex(num), srcIndex(num), IndexPhi(num)));
					num++;
				}
			}
		}
		/*<2A><>ʱP<CAB1><50><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>*/
		std::sort(temp_P_complex.begin(), temp_P_complex.end(),
			[](const Eigen::Triplet<complex<double>>& a, const Eigen::Triplet<complex<double>>& b)
			{
				return a.col() < b.col();
			});


		/*ɾ<><C9BE>dstIndex<65><78>Ӧ<EFBFBD><D3A6>*/
		std::vector<Eigen::Triplet<complex<double>>> tripleP_complex;
		num = 0;// <20><><EFBFBD>ڸ<EFBFBD><DAB8><EFBFBD><EFBFBD>¾<EFBFBD><C2BE><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
		int	DeleteIndex = 0;
		for (int i = 0; i < dof; i++)
		{
			DeleteIndex = 0;
			for (int j = 0; j < dstIndex.rows(); j++)
			{
				if (i == dstIndex[j])
				{
					DeleteIndex = 1;
				}
			}
			if (DeleteIndex != 1)
			{
				for (const auto& triplet : temp_P_complex)
				{
					if (triplet.col() == i)
					{
						tripleP_complex.push_back(Eigen::Triplet<complex<double>>(triplet.row(), triplet.col() - num, triplet.value()));
					}
				}
			}
			for (int j = 0; j < dstIndex.rows(); j++)
			{
				if (i == dstIndex[j])
				{
					num = num + 1;
				}
			}

		}
		_mP_complex = Eigen::SparseMatrix<complex<double>, Eigen::RowMajor>(dof, dof - dstIndex.rows());
		_mP_complex.setFromTriplets(tripleP_complex.begin(), tripleP_complex.end());;
	}
}

void OpticsFEM_2D_Scatter::Assemble_Port()
{
	//physic
	double k0 = 2 * Pi / _mSolver->GetLda0();

	//init of Gauss point
	Gauss gauss;
	int NbrGuassPoints;
	NbrGuassPoints = gauss.GetNbrGaussPoints(ONEDIM, LINE, BF_LINEFUNC * 4);
	double* u, * v, * w, * wght;
	u = new double[NbrGuassPoints];
	v = new double[NbrGuassPoints];
	w = new double[NbrGuassPoints];
	wght = new double[NbrGuassPoints];
	gauss.GetGaussPoints(ONEDIM, LINE, u, v, w, wght);
	//init of geo
	Vector3d* vertex = new Vector3d[5];
	Matrix3d Jac, InvJac;
	Vector3d normal;

	//init of basis function
	BF BF_Lagrange, BF_Nedelec;
	Vector3d* Et, ** Ez;
	BF_Lagrange.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	BF_Nedelec.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	Et = new Vector3d[NbrGuassPoints];
	Ez = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		Ez[i] = new Vector3d[2];
	}

	//inite parser
	mup::ParserX p(mup::pckALL_COMPLEX);
	p.EnableAutoCreateVar(true);
	mup::Value xx, yy;
	p.DefineVar(_T("x"), mup::Variable(&xx));
	p.DefineVar(_T("y"), mup::Variable(&yy));
	Vector3cd* N0, * N2, * curlN0, * curlN1, * curlN2;
	N0 = new Vector3cd[NbrGuassPoints];
	N2 = new Vector3cd[NbrGuassPoints];
	curlN0 = new Vector3cd[NbrGuassPoints];
	curlN1 = new Vector3cd[NbrGuassPoints];
	curlN2 = new Vector3cd[NbrGuassPoints];

	int NbrPortinc = _mPhy->GetNbrPortinc();
	for (int n = 0; n < NbrPortinc; n++)
	{
		int PortincDomain = _mPhy->GetPortincDomain(n) + 1;
		Eigen::VectorXi edgesIndex;
		_mMesh->GetEdgesIndexOfDomain(PortincDomain, edgesIndex);
		string Eincx, Eincy, Eincz, Eoutx, Eouty, Eoutz, EincCurlx, EincCurly, EincCurlz, EoutCurlx, EoutCurly, EoutCurlz;
		_mPhy->GetPortEinc(Eincx, Eincy, Eincz);
		_mPhy->GetPortEout(Eoutx, Eouty, Eoutz);
		_mPhy->GetPortEincCurl(EincCurlx, EincCurly, EincCurlz);
		_mPhy->GetPortEoutCurl(EoutCurlx, EoutCurly, EoutCurlz);
		//overloop of edgesOfPortinc
		for (int m = 0; m < edgesIndex.rows(); m++)
		{
			// edge data
			Vector2i coonOfEdges;
			_mMesh->GetCoonOfEdges(edgesIndex(m), coonOfEdges);
			int numEle = coonOfEdges(0);
			int numEdge = coonOfEdges(1);

			//coordinate of vertex,integration coordinate, and mapping
			int MappingIndex[3];
			int index[3];
			for (int i = 0; i < 3; i++)
				_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
			if (numEdge == 0)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[1];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = 0.0;
				}

				index[0] = 1; index[1] = 2; index[2] = 1;
				MappingIndex[0] = _mMesh->GetTri(numEle, 0);
				MappingIndex[1] = _mMesh->GetTri(numEle, 1);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 0) + _mMesh->GetNbrVertex();
			}
			else if (numEdge == 1)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = 0.0;
					w[i] = (u[i] + 1.0) / 2.0;
				}

				index[0] = 1; index[1] = 3; index[2] = 2;
				MappingIndex[0] = _mMesh->GetTri(numEle, 0);
				MappingIndex[1] = _mMesh->GetTri(numEle, 2);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 1) + _mMesh->GetNbrVertex();
			}
			else if (numEdge == 2)
			{
				vertex[3] = vertex[0];
				vertex[4] = vertex[2];

				for (int i = 0; i < NbrGuassPoints; i++)
				{
					v[i] = (u[i] + 1.0) / 2.0;
					w[i] = (1.0 - u[i]) / 2.0;
				}

				index[0] = 2; index[1] = 3; index[2] = 3;
				MappingIndex[0] = _mMesh->GetTri(numEle, 1);
				MappingIndex[1] = _mMesh->GetTri(numEle, 2);
				MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 2) + _mMesh->GetNbrVertex();
			}

			//Jac
			Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
			Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
			Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
			InvJac = Jac.inverse();

			//integration coe
			double phi1 = vertex[4](1) - vertex[3](1);
			double phi2 = vertex[3](0) - vertex[4](0);
			double integCoe = 0;
			if (vertex[3](0) != vertex[4](0))
			{
				integCoe = phi1 / phi2;
				integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
			}
			else
			{
				integCoe = abs(phi1) / 2.0;
			}

			//normal
			_mMesh->GetNormOfEdges(_mMesh->GetDomainOfEdges(edgesIndex(m)) - 1, normal);

			//basis function
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				for (int j = 0; j < 2; j++)
				{
					BF_Lagrange.GetValueBF(index[j], v[i], w[i], 0.0, Ez[i][j]);
				}
				BF_Nedelec.GetValueBF(index[2], v[i], w[i], 0.0, Et[i]);
				Et[i] = InvJac * Et[i];
			}
			for (int i = 0; i < NbrGuassPoints; i++)
			{
				double tempX = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
				double tempY = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);
				xx = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
				yy = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);
				//<2F><><EFBFBD>䳡N0
				p.SetExpr(_T(Eincx));
				N0[i](0).real(p.Eval().GetFloat());
				N0[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(Eincy));
				N0[i](1).real(p.Eval().GetFloat());
				N0[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(Eincz));
				N0[i](2).real(p.Eval().GetFloat());
				N0[i](2).imag(p.Eval().GetImag());

				//<2F><><EFBFBD>䳡N2
				p.SetExpr(_T(Eoutx));
				N2[i](0).real(p.Eval().GetFloat());
				N2[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(Eincy));
				N0[i](1).real(p.Eval().GetFloat());
				N0[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(Eincz));
				N0[i](2).real(p.Eval().GetFloat());
				N0[i](2).imag(p.Eval().GetImag());

				//<2F><><EFBFBD>䳡N0<4E><30><EFBFBD><EFBFBD>
				p.SetExpr(_T(EincCurlx));
				curlN0[i](0).real(p.Eval().GetFloat());
				curlN0[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(EincCurly));
				curlN0[i](1).real(p.Eval().GetFloat());
				curlN0[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(EincCurlz));
				curlN0[i](2).real(p.Eval().GetFloat());
				curlN0[i](2).imag(p.Eval().GetImag());

				//<2F><><EFBFBD>䳡N1<4E><31><EFBFBD><EFBFBD>
				curlN1[i](0) = std::conj(curlN0[i](0));
				curlN1[i](1) = std::conj(curlN0[i](1));
				curlN1[i](2) = std::conj(curlN0[i](2));
				//<2F><><EFBFBD>䳡N2
				p.SetExpr(_T(EoutCurlx));
				curlN2[i](0).real(p.Eval().GetFloat());
				curlN2[i](0).imag(p.Eval().GetImag());
				p.SetExpr(_T(EoutCurly));
				curlN2[i](1).real(p.Eval().GetFloat());
				curlN2[i](1).imag(p.Eval().GetImag());
				p.SetExpr(_T(EoutCurlz));
				curlN2[i](2).real(p.Eval().GetFloat());
				curlN2[i](2).imag(p.Eval().GetImag());
			}

			//material 
			int domain = _mMesh->GetDomainOfTri(n);
			Matrix3cd epsr = _mMatLib->GetEpsr(domain);
			Matrix3cd refra = epsr.cwiseSqrt();

			//submatrix
			MatrixXcd Se1, Ge1, be1;
			Se1 = MatrixXcd::Zero(2, 1);
			Ge1 = MatrixXcd::Zero(2, 1);
			be1 = MatrixXcd::Zero(2, 1);
			complex<double>Se2, Ge2, be2;
			for (int i = 0; i < 2; i++)
			{
				for (int k = 0; k < NbrGuassPoints; k++)
				{
					be1(i, 0) = be1(i, 0) + wght[k] * integCoe * Ez[k][i].dot(curlN0[k].cross(refra * normal));
					Se1(i, 0) = Se1(i, 0) + wght[k] * integCoe * Ez[k][i].dot(curlN1[k].cross(refra * normal));
					Ge1(i, 0) = Ge1(i, 0) + wght[k] * integCoe * N0[k].dot(normal.cross(Ez[k][i].cross(normal)));
				}
			}
			for (int k = 0; k < NbrGuassPoints; k++)
			{
				be2 = be2 + wght[k] * integCoe * Et[k].dot(curlN0[k].cross(refra * normal));
				Se2 = Se2 + wght[k] * integCoe * Et[k].dot(curlN1[k].cross(refra * normal));
				Ge2 = Ge2 + wght[k] * integCoe * N0[k].dot(normal.cross(Et[k].cross(normal)));
				P1 = P1 + wght[k] * integCoe * N0[k].dot(normal.cross(N0[k].cross(normal)));
				b1 = b1 + wght[k] * integCoe * N0[k].dot(normal.cross(N0[k].cross(normal)));
			}
			//<2F><>װ<EFBFBD><D7B0><EFBFBD><EFBFBD>
			for (int i = 0; i < 2; i++)
			{
				_mPort_B_complex(MappingIndex[i]) = _mPort_B_complex(MappingIndex[i]) - be1(i, 0);
				_mS_complex(MappingIndex[i]) = _mS_complex(MappingIndex[i]) + Se1(i, 0);
				_mG1_complex(MappingIndex[i]) = _mG1_complex(MappingIndex[i]) + Ge1(i, 0);
			}
			_mPort_B_complex(MappingIndex[2]) = _mPort_B_complex(MappingIndex[2]) - be2;
			_mS_complex(MappingIndex[2]) = _mS_complex(MappingIndex[2]) + Se2;
			_mG1_complex(MappingIndex[2]) = _mG1_complex(MappingIndex[2]) + Ge2;


			//<2F><><EFBFBD><EFBFBD><EFBFBD>˿<EFBFBD>
			int PortoutDomain = _mPhy->GetPortoutDomain(n) + 1;
			Eigen::VectorXi edgesIndex;
			_mMesh->GetEdgesIndexOfDomain(PortoutDomain, edgesIndex);
			for (int m = 0; m < edgesIndex.rows(); m++)
			{
				// edge data
				Vector2i coonOfEdges;
				_mMesh->GetCoonOfEdges(edgesIndex(m), coonOfEdges);
				int numEle = coonOfEdges(0);
				int numEdge = coonOfEdges(1);

				//coordinate of vertex,integration coordinate, and mapping
				int MappingIndex[3];
				int index[3];
				for (int i = 0; i < 3; i++)
					_mMesh->GetVertex(_mMesh->GetTri(numEle, i), vertex[i]);
				if (numEdge == 0)
				{
					vertex[3] = vertex[0];
					vertex[4] = vertex[1];

					for (int i = 0; i < NbrGuassPoints; i++)
					{
						v[i] = (u[i] + 1.0) / 2.0;
						w[i] = 0.0;
					}

					index[0] = 1; index[1] = 2; index[2] = 1;
					MappingIndex[0] = _mMesh->GetTri(numEle, 0);
					MappingIndex[1] = _mMesh->GetTri(numEle, 1);
					MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 0) + _mMesh->GetNbrVertex();
				}
				else if (numEdge == 1)
				{
					vertex[3] = vertex[0];
					vertex[4] = vertex[2];

					for (int i = 0; i < NbrGuassPoints; i++)
					{
						v[i] = 0.0;
						w[i] = (u[i] + 1.0) / 2.0;
					}

					index[0] = 1; index[1] = 3; index[2] = 2;
					MappingIndex[0] = _mMesh->GetTri(numEle, 0);
					MappingIndex[1] = _mMesh->GetTri(numEle, 2);
					MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 1) + _mMesh->GetNbrVertex();
				}
				else if (numEdge == 2)
				{
					vertex[3] = vertex[0];
					vertex[4] = vertex[2];

					for (int i = 0; i < NbrGuassPoints; i++)
					{
						v[i] = (u[i] + 1.0) / 2.0;
						w[i] = (1.0 - u[i]) / 2.0;
					}

					index[0] = 2; index[1] = 3; index[2] = 3;
					MappingIndex[0] = _mMesh->GetTri(numEle, 1);
					MappingIndex[1] = _mMesh->GetTri(numEle, 2);
					MappingIndex[2] = _mMesh->GetEdgeOfTri(numEle, 2) + _mMesh->GetNbrVertex();
				}

				//Jac
				Jac(0, 0) = vertex[1](0) - vertex[0](0); Jac(0, 1) = vertex[1](1) - vertex[0](1); Jac(0, 2) = 0.;
				Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = 0.;
				Jac(2, 0) = 0.;			Jac(2, 1) = 0.;			Jac(2, 2) = 1.;
				InvJac = Jac.inverse();

				//integration coe
				double phi1 = vertex[4](1) - vertex[3](1);
				double phi2 = vertex[3](0) - vertex[4](0);
				double integCoe = 0;
				if (vertex[3](0) != vertex[4](0))
				{
					integCoe = phi1 / phi2;
					integCoe = fabs(sqrtf(1 + integCoe * integCoe) * phi2 / 2.0);
				}
				else
				{
					integCoe = abs(phi1) / 2.0;
				}

				//normal
				_mMesh->GetNormOfEdges(_mMesh->GetDomainOfEdges(edgesIndex(m)) - 1, normal);

				//basis function
				for (int i = 0; i < NbrGuassPoints; i++)
				{
					for (int j = 0; j < 2; j++)
					{
						BF_Lagrange.GetValueBF(index[j], v[i], w[i], 0.0, Ez[i][j]);
					}
					BF_Nedelec.GetValueBF(index[2], v[i], w[i], 0.0, Et[i]);
					Et[i] = InvJac * Et[i];
				}
				for (int i = 0; i < NbrGuassPoints; i++)
				{
					double tempX = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
					double tempY = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);
					xx = Jac(0, 0) * v[i] + Jac(1, 0) * w[i] + vertex[0](0);
					yy = Jac(0, 1) * v[i] + Jac(1, 1) * w[i] + vertex[0](1);
					//<2F><><EFBFBD>䳡N0
					p.SetExpr(_T(Eincx));
					N0[i](0).real(p.Eval().GetFloat());
					N0[i](0).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eincy));
					N0[i](1).real(p.Eval().GetFloat());
					N0[i](1).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eincz));
					N0[i](2).real(p.Eval().GetFloat());
					N0[i](2).imag(p.Eval().GetImag());

					//<2F><><EFBFBD>䳡N2
					p.SetExpr(_T(Eoutx));
					N2[i](0).real(p.Eval().GetFloat());
					N2[i](0).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eouty));
					N2[i](1).real(p.Eval().GetFloat());
					N2[i](1).imag(p.Eval().GetImag());
					p.SetExpr(_T(Eoutz));
					N2[i](2).real(p.Eval().GetFloat());
					N2[i](2).imag(p.Eval().GetImag());

					//<2F><><EFBFBD>䳡N0<4E><30><EFBFBD><EFBFBD>
					p.SetExpr(_T(EincCurlx));
					curlN0[i](0).real(p.Eval().GetFloat());
					curlN0[i](0).imag(p.Eval().GetImag());
					p.SetExpr(_T(EincCurly));
					curlN0[i](1).real(p.Eval().GetFloat());
					curlN0[i](1).imag(p.Eval().GetImag());
					p.SetExpr(_T(EincCurlz));
					curlN0[i](2).real(p.Eval().GetFloat());
					curlN0[i](2).imag(p.Eval().GetImag());

					//<2F><><EFBFBD>䳡N1<4E><31><EFBFBD><EFBFBD>
					curlN1[i](0) = std::conj(curlN0[i](0));
					curlN1[i](1) = std::conj(curlN0[i](1));
					curlN1[i](2) = std::conj(curlN0[i](2));
					//<2F><><EFBFBD>䳡N2
					p.SetExpr(_T(EoutCurlx));
					curlN2[i](0).real(p.Eval().GetFloat());
					curlN2[i](0).imag(p.Eval().GetImag());
					p.SetExpr(_T(EoutCurly));
					curlN2[i](1).real(p.Eval().GetFloat());
					curlN2[i](1).imag(p.Eval().GetImag());
					p.SetExpr(_T(EoutCurlz));
					curlN2[i](2).real(p.Eval().GetFloat());
					curlN2[i](2).imag(p.Eval().GetImag());
				}

				//material 
				int domain = _mMesh->GetDomainOfTri(n);
				Matrix3cd epsr = _mMatLib->GetEpsr(domain);
				Matrix3cd refra = epsr.cwiseSqrt();

				//submatrix
				MatrixXcd Te1, Ge1;
				Te1 = MatrixXcd::Zero(2, 1);
				Ge1 = MatrixXcd::Zero(2, 1);
				complex<double>P2;
				for (int i = 0; i < 2; i++)
				{
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Te1(i, 0) = Te1(i, 0) + wght[k] * integCoe * Ez[k][i].dot(curlN2[k].cross(refra * normal));
						Ge1(i, 0) = Ge1(i, 0) + wght[k] * integCoe * N2[k].dot(normal.cross(Ez[k][i].cross(normal)));
						P2 = P2 + wght[k] * integCoe * N2[k].dot(normal.cross(N2[k].cross(normal)));
					}
				}
				complex<double>Te2, Ge2;
				for (int k = 0; k < NbrGuassPoints; k++)
				{
					Te2 = Te2 + wght[k] * integCoe * Et[k].dot(curlN2[k].cross(refra * normal));
					Ge2 = Ge2 + wght[k] * integCoe * N2[k].dot(normal.cross(Et[k].cross(normal)));
				}

				//<2F><>װ<EFBFBD><D7B0><EFBFBD><EFBFBD>
				for (int i = 0; i < 2; i++)
				{
					_mT_complex(MappingIndex[i]) = _mT_complex(MappingIndex[i]) + Te1(i, 0);
					_mG2_complex(MappingIndex[i]) = _mG2_complex(MappingIndex[i]) + Ge1(i, 0);
				}
				_mT_complex(MappingIndex[2]) = _mT_complex(MappingIndex[2]) + Te2;
				_mG2_complex(MappingIndex[2]) = _mG2_complex(MappingIndex[2]) + Ge2;
			}
		}
	}

	_mS_complex = _mP_complex.adjoint() * _mS_complex;
	_mT_complex = _mP_complex.adjoint() * _mT_complex;
	_mG1_complex = _mP_complex.adjoint() * _mG1_complex;
	_mG2_complex = _mP_complex.adjoint() * _mG2_complex;

	for (const auto& triplet : _mTripleA_complex)
	{
		_mTriplePort_A_complex.push_back(Triplet<complex<double>>(triplet.row(), triplet.col(), triplet.value()));
		/*std::cout << "col: " << triplet.col()  << "row: " << triplet.row() << ", Value: " << triplet.value() << std::endl;*/
	}
	for (int k = 0; k < _mA_complex.outerSize(); k++)
	{
		for (Eigen::SparseMatrix<complex<double>, RowMajor>::InnerIterator it(_mA_complex, k); it; ++it)
		{
			_mTriplePort_A_complex.push_back(Triplet<complex<double>>(it.row(), it.col(), it.value()));
		}
	}
	for (int i = 0; i < _mS_complex.size(); i++) {
		_mTriplePort_A_complex.push_back(Triplet<complex<double>>(i, _mDof, _mS_complex[i]));
	}
	for (int i = 0; i < _mT_complex.size(); i++) {
		_mTriplePort_A_complex.push_back(Triplet<complex<double>>(i, _mDof + 1, _mT_complex[i]));
	}
	for (int i = 0; i < _mG1_complex.size(); i++) {
		_mTriplePort_A_complex.push_back(Triplet<complex<double>>(_mDof, i, _mG1_complex[i]));
	}
	for (int i = 0; i < _mG2_complex.size(); i++) {
		_mTriplePort_A_complex.push_back(Triplet<complex<double>>(_mDof + 1, i, _mG2_complex[i]));
	}
	_mTriplePort_A_complex.push_back(Triplet<complex<double>>(_mDof, _mDof, -P1));
	_mTriplePort_A_complex.push_back(Triplet<complex<double>>(_mDof + 1, _mDof + 1, -P2));


	_mPort_B_complex(_mDof) = b1;
	for (int k = 0; k < _mP_complex.outerSize(); k++)
	{
		for (Eigen::SparseMatrix<complex<double>, RowMajor>::InnerIterator it(_mP_complex, k); it; ++it)
		{
			_mPort_P_complex.coeffRef(it.row(), it.col()) = it.value();
		}
	}
}