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

#include <vector>
#include <algorithm>
#include <string>

using namespace std;
using namespace Eigen;


void OpticsFEM_2D_Scatter::Assemble_WaveEquation()
{
	//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[3];
	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;
	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];
	}

	//loop over tri
	complex<double> iUnit{ 0,1 };
	int NbrTri = _mMesh->GetNbrTri();
	for (int n = 0; n < NbrTri; n++)
	{
		//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];
			}
		}

		//material
		int domain = _mMesh->GetDomainOfTri(n);
		Matrix3cd epsr = _mMatLib->GetEpsr(domain);
		Matrix3d sigma = _mMatLib->GetSigma(domain);
		epsr = epsr - sigma * complex<double>(0.0, 1.0 / k0 * 120.0 * Pi);
		Matrix3cd Mur = _mMatLib->GetMur(domain);
		Matrix3cd chihe = _mMatLib->GetChihe(domain);
		Matrix3cd chieh = _mMatLib->GetChieh(domain);
		//PML
		Matrix3cd invMur;
		Eigen::VectorXd PMLData;
		int PMLType;
		if (_mPhy->GetPML(domain, PMLType, PMLData))
		{
			int R0 = 10;
			double averX = (vertex[0](0) + vertex[1](0) + vertex[2](0)) / 3.0;
			double averY = (vertex[0](1) + vertex[1](1) + vertex[2](1)) / 3.0;
			Matrix3cd Lambda = Matrix3cd::Zero();
			if (PMLType == 0)
			{
				complex<double>  sx{ 1,-fabs((averX - PMLData(0)) / PMLData(1)) * R0 / PMLData(1) / k0 };
				complex<double>  sy{ 1,-fabs((averY - PMLData(2)) / PMLData(3)) * R0 / PMLData(3) / k0 };
				Lambda(0, 0) = sy / sx; Lambda(1, 1) = sx / sy; Lambda(2, 2) = sx * sy;
			}
			else if (PMLType == 1)
			{
				double rho = sqrt(averX * averX + averY * averY);
				double sigma = pow((rho - PMLData(1)) / PMLData(3), 2) * R0 / PMLData(3) / k0;
				complex<double> s1{ 1, -PMLData(3) / 2 / rho * sigma };
				complex<double> s2{ 1, -sigma };
				complex<double> aa = s1 / s2; complex<double> bb = s2 / s1; complex<double> cc = s1 * s2;
				Lambda(0, 0) = (aa * averX * averX + bb * averY * averY) / rho / rho;
				Lambda(0, 1) = (aa - bb) * averX * averY / rho / rho;
				Lambda(1, 0) = (aa - bb) * averX * averY / rho / rho;
				Lambda(1, 1) = (bb * averX * averX + aa * averY * averY) / rho / rho;
				Lambda(2, 2) = cc;
			}
			epsr = epsr * Lambda;
			invMur = (Mur * Lambda).inverse();
		}
		else
		{
			invMur = Mur.inverse();
		}

		//mapping 
		VectorXi MappingIndexS = VectorXi::Zero(sdof);
		VectorXi MappingIndexV = VectorXi::Zero(vdof);
		for (int i = 0; i < sdof; i++)
			MappingIndexS(i) = _mMesh->GetTri(n, i);
		for (int i = 0; i < vdof; i++)
			MappingIndexV(i) = _mMesh->GetEdgeOfTri(n, i) + _mMesh->GetNbrVertex();

		//submatrix
		MatrixXcd St, Sz, Tt, Tz, Yt, Yz, Ft, Fz, Gt, Gz;
		St = MatrixXcd::Zero(vdof, vdof);
		Sz = MatrixXcd::Zero(sdof, sdof);
		Tt = MatrixXcd::Zero(vdof, vdof);
		Tz = MatrixXcd::Zero(sdof, sdof);
		Yt = MatrixXcd::Zero(vdof, vdof);
		Yz = MatrixXcd::Zero(sdof, sdof);
		Ft = MatrixXcd::Zero(sdof, vdof);
		Fz = MatrixXcd::Zero(vdof, sdof);
		Gt = MatrixXcd::Zero(sdof, vdof);
		Gz = MatrixXcd::Zero(vdof, sdof);
		
		if (_mPhy->GetBeamState() == 1)
		{
			double kx, ky;
			Vector3cd curlEEz, curlWWz, curlkir;
			_mPhy->GetBeamDir(kx, ky);
			curlkir = Vector3cd::Zero();
			curlkir[0].imag(ky * k0);
			curlkir[1].imag(-kx * k0);

			//????????? ??????Ez????
			for (int i = 0; i < sdof; i++)
				for (int j = 0; j < sdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						curlEEz = curlEz[k][j] + curlkir * Ez[k][j][2];
						curlWWz = curlEz[k][i] - curlkir * Ez[k][i][2];
						Sz(i, j) = Sz(i, j) + wght[k] * DetJac * curlWWz.dot(invMur * curlEEz.conjugate());
						Tz(i, j) = Tz(i, j) + wght[k] * DetJac * k0 * k0 * Ez[k][i].dot(epsr * Ez[k][j]);
						Yz(i, j) = 0.;
					}
			for (int i = 0; i < vdof; i++)
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						St(i, j) = 0.;
						Tt(i, j) = 0.;
						Yt(i, j) = 0.;
					}
			for (int i = 0; i < sdof; i++)
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Ft(i, j) = 0.;
						Fz(j, i) = 0.;
						Gt(i, j) = 0.;
						Gz(j, i) = 0.;
					}
		}
		else
		{
			for (int i = 0; i < sdof; i++)
				for (int j = 0; j < sdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Sz(i, j) = Sz(i, j) + wght[k] * DetJac * curlEz[k][i].dot(invMur * curlEz[k][j]);
						Tz(i, j) = Tz(i, j) + wght[k] * DetJac * k0 * k0 * Ez[k][i].dot(epsr * Ez[k][j]);
						Yz(i, j) = Yz(i, j) - wght[k] * DetJac * (iUnit * k0 * chieh * Ez[k][i]).dot(iUnit * k0 * invMur * chihe * Ez[k][j]);
					}
			for (int i = 0; i < vdof; i++)
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						St(i, j) = St(i, j) + wght[k] * DetJac * curlEt[k][i].dot(invMur * curlEt[k][j]);
						Tt(i, j) = Tt(i, j) + wght[k] * DetJac * k0 * k0 * Et[k][i].dot(epsr * Et[k][j]);
						Yt(i, j) = Yt(i, j) - wght[k] * DetJac * (iUnit * k0 * chieh * Et[k][i]).dot(iUnit * k0 * invMur * chihe * Et[k][j]);
					}
			for (int i = 0; i < sdof; i++)
				for (int j = 0; j < vdof; j++)
					for (int k = 0; k < NbrGuassPoints; k++)
					{
						Ft(i, j) = Ft(i, j) - wght[k] * DetJac * (iUnit * k0 * chieh * Ez[k][i]).dot(invMur * curlEt[k][j]);
						Fz(j, i) = Fz(j, i) - wght[k] * DetJac * (iUnit * k0 * chieh * Et[k][j]).dot(invMur * curlEz[k][i]);
						Gt(i, j) = Gt(i, j) + wght[k] * DetJac * curlEz[k][i].dot(iUnit * k0 * invMur * chihe * Et[k][j]);
						Gz(j, i) = Gz(j, i) + wght[k] * DetJac * curlEt[k][j].dot(iUnit * k0 * invMur * chihe * Ez[k][i]);
					}
		}
		


		//store in triplet
		if (_mIsReal)
		{
			for (int i = 0; i < sdof; i++)
			{
				for (int j = 0; j < sdof; j++)
				{
					_mTripleA_real.push_back(Triplet<double>(MappingIndexS(i), MappingIndexS(j), Sz(i, j).real() - Tz(i, j).real()));
				}
			}
			for (int i = 0; i < vdof; i++)
			{
				for (int j = 0; j < vdof; j++)
				{
					_mTripleA_real.push_back(Triplet<double>(MappingIndexV(i), MappingIndexV(j), St(i, j).real() - Tt(i, j).real()));
				}
			}

		}
		else
		{
			for (int i = 0; i < sdof; i++)
			{
				for (int j = 0; j < sdof; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndexS(i), MappingIndexS(j), Sz(i, j) + Yz(i, j) - Tz(i, j)));
				}
			}
			for (int i = 0; i < vdof; i++)
			{
				for (int j = 0; j < vdof; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndexV(i), MappingIndexV(j), St(i, j) + Yt(i, j) - Tt(i, j)));
				}
			}
			for (int i = 0; i < sdof; i++)
				for (int j = 0; j < vdof; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndexS(i), MappingIndexV(j), Ft(i, j) + Gt(i, j)));
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndexV(j), MappingIndexS(i), Fz(j, i) + Gz(j, 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;
}

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

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

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

	//init of basis function
	BF BF_Nedelec, BF_Curl_Nedelec;
	int dof;
	Vector3d** E, ** curlE;
	dof = BF_Nedelec.GetNbrBF(THREEDIM, TETRAHEDRON, BF_NEDELEC, BF_LINEFUNC);
	BF_Curl_Nedelec.GetNbrBF(THREEDIM, TETRAHEDRON, BF_CURL_NEDELEC, BF_LINEFUNC);
	E = new Vector3d * [NbrGuassPoints];
	curlE = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		E[i] = new Vector3d[dof];
		curlE[i] = new Vector3d[dof];
	}

	//loop over tri
	int NbrTet = _mMesh->GetNbrTet();
	for (int n = 0; n < NbrTet; n++)
	{
		//coordinate of vertex
		for (int i = 0; i < 4; i++)
			_mMesh->GetVertex(_mMesh->GetTet(n, i), vertex[i]);

		//Jac - 3D reference tetrahedron mapping (vertex 4 = index 3)
		Jac(0, 0) = vertex[0](0) - vertex[3](0); Jac(0, 1) = vertex[0](1) - vertex[3](1); Jac(0, 2) = vertex[0](2) - vertex[3](2);
		Jac(1, 0) = vertex[1](0) - vertex[3](0); Jac(1, 1) = vertex[1](1) - vertex[3](1); Jac(1, 2) = vertex[1](2) - vertex[3](2);
		Jac(2, 0) = vertex[2](0) - vertex[3](0); Jac(2, 1) = vertex[2](1) - vertex[3](1); Jac(2, 2) = vertex[2](2) - vertex[3](2);
		double DetJac = fabs(Jac.determinant());
		InvJac = Jac.inverse();
		TJac = Jac.transpose() / Jac.determinant();

		//basis function
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			for (int j = 0; j < dof; j++)
			{
				BF_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], E[i][j]);
				E[i][j] = InvJac * E[i][j];
				BF_Curl_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], curlE[i][j]);
				curlE[i][j] = TJac * curlE[i][j];
			}
		}

		//material
		int domain = _mMesh->GetDomainOfTet(n);
		Matrix3cd epsr = _mMatLib->GetEpsr(domain);
		Matrix3d sigma = _mMatLib->GetSigma(domain);
		epsr = epsr - sigma * complex<double>(0.0, 1.0 / k0 * 120.0 * Pi);
		Matrix3cd Mur = _mMatLib->GetMur(domain);
		//PML
		Matrix3cd invMur;
		Eigen::VectorXd PMLData;
		int PMLType;
		if (_mPhy->GetPML(domain, PMLType, PMLData))
		{
			int R0 = 10;
			double averX = (vertex[0](0) + vertex[1](0) + vertex[2](0)) / 3.0;
			double averY = (vertex[0](1) + vertex[1](1) + vertex[2](1)) / 3.0;
			double averZ = (vertex[0](2) + vertex[1](2) + vertex[2](2)) / 3.0;
			complex<double>  sx{ 1,-fabs((averX - PMLData(0)) * PMLData(1)) * R0 * PMLData(1) / k0 };
			complex<double>  sy{ 1,-fabs((averY - PMLData(2)) * PMLData(3)) * R0 * PMLData(3) / k0 };
			complex<double>  sz{ 1,-fabs((averZ - PMLData(4)) * PMLData(5)) * R0 * PMLData(5) / k0 };
			Matrix3cd Lambda = Matrix3cd::Zero();
			Lambda(0, 0) = sy * sz / sx; Lambda(1, 1) = sx * sz / sy; Lambda(2, 2) = sx * sy / sz;
			epsr = epsr * Lambda;
			invMur = (Mur * Lambda).inverse();
		}
		else
		{
			invMur = Mur.inverse();
		}

		//mapping 
		VectorXi MappingIndex = VectorXi::Zero(dof);
		for (int i = 0; i < dof; i++)
			MappingIndex(i) = _mMesh->GetEdgeOfTet(n, i);

		//submatrix
		MatrixXcd Se, Te;
		Se = MatrixXcd::Zero(dof, dof);
		Te = MatrixXcd::Zero(dof, dof);
		for (int i = 0; i < dof; i++)
			for (int j = 0; j < dof; j++)
				for (int k = 0; k < NbrGuassPoints; k++)
				{
					Se(i, j) = Se(i, j) + wght[k] * DetJac * curlE[k][i].dot(invMur * curlE[k][j]);
					Te(i, j) = Te(i, j) + wght[k] * DetJac * k0 * k0 * E[k][i].dot(epsr * E[k][j]);
				}


		//store in triplet
		if (_mIsReal)
		{
			for (int i = 0; i < dof; i++)
			{
				for (int j = 0; j < dof; j++)
				{
					_mTripleA_real.push_back(Triplet<double>(MappingIndex(i), MappingIndex(j), Se(i, j).real() - Te(i, j).real()));
				}
			}

		}
		else
		{
			for (int i = 0; i < dof; i++)
			{
				for (int j = 0; j < dof; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndex(i), MappingIndex(j), Se(i, j) - Te(i, j)));
				}
			}

		}
	}

	delete[] u, v, w, wght;

	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] E[i], curlE[i];
	}
	delete[] E, curlE;
}

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

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

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

	//init of basis function
	BF BF_Nedelec, BF_Curl_Nedelec;
	int dof;
	Vector3d** E, ** curlE;
	dof = BF_Nedelec.GetNbrBF(THREEDIM, PRISM, BF_NEDELEC, BF_LINEFUNC);
	BF_Curl_Nedelec.GetNbrBF(THREEDIM, PRISM, BF_CURL_NEDELEC, BF_LINEFUNC);
	E = new Vector3d * [NbrGuassPoints];
	curlE = new Vector3d * [NbrGuassPoints];
	for (int i = 0; i < NbrGuassPoints; i++)
	{
		E[i] = new Vector3d[dof];
		curlE[i] = new Vector3d[dof];
	}

	//loop over tri
	int NbrPrism = _mMesh->GetNbrPrism();
	for (int n = 0; n < NbrPrism; n++)
	{
		//coordinate of vertex
		for (int i = 0; i < 6; i++)
			_mMesh->GetVertex(_mMesh->GetPrism(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) = vertex[1](2)-vertex[0](2);
		Jac(1, 0) = vertex[2](0) - vertex[0](0); Jac(1, 1) = vertex[2](1) - vertex[0](1); Jac(1, 2) = vertex[2](2) - vertex[0](2);
		Jac(2, 0) = (vertex[3](0) - vertex[0](0)) / 2.; Jac(2, 1) = (vertex[3](1) - vertex[0](1)) / 2.; Jac(2, 2) = (vertex[3](2) - vertex[0](2)) / 2.;
		double d[3] = { (vertex[3](0) + vertex[0](0)) / 2.,(vertex[3](1) + vertex[0](1)) / 2.,(vertex[3](2) + vertex[0](2)) / 2. };
		double DetJac = fabs(Jac.determinant());
		InvJac = Jac.inverse();
		TJac = Jac.transpose() / Jac.determinant();

		//basis function
		for (int i = 0; i < NbrGuassPoints; i++)
		{
			for (int j = 0; j < dof; j++)
			{
				BF_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], E[i][j]);
				E[i][j] = InvJac * E[i][j];
				BF_Curl_Nedelec.GetValueBF(j + 1, u[i], v[i], w[i], curlE[i][j]);
				curlE[i][j] = TJac * curlE[i][j];
			}
		}

		//material
		int domain = _mMesh->GetDomainOfPrism(n);
		Matrix3cd epsr = _mMatLib->GetEpsr(domain);
		Matrix3d sigma = _mMatLib->GetSigma(domain);
		epsr = epsr - sigma * complex<double>(0.0, 1.0 / k0 * 120.0 * Pi);
		Matrix3cd Mur = _mMatLib->GetMur(domain);
		//PML
		Matrix3cd invMur;
		Eigen::VectorXd PMLData;
		int PMLType;
		if (_mPhy->GetPML(domain,PMLType,PMLData))
		{
			int R0 = 10;
			double averX = (vertex[0](0) + vertex[1](0) + vertex[2](0)) / 3.0;
			double averY = (vertex[0](1) + vertex[1](1) + vertex[2](1)) / 3.0;
			double averZ = (vertex[0](2) + vertex[1](2) + vertex[2](2)) / 3.0;
			complex<double>  sx{ 1,-fabs((averX - PMLData(0)) * PMLData(1)) * R0 * PMLData(1) / k0 };
			complex<double>  sy{ 1,-fabs((averY - PMLData(2)) * PMLData(3)) * R0 * PMLData(3) / k0 };
			complex<double>  sz{ 1,-fabs((averZ - PMLData(4)) * PMLData(5)) * R0 * PMLData(5) / k0 };
			Matrix3cd Lambda = Matrix3cd::Zero();
			Lambda(0, 0) = sy * sz / sx; Lambda(1, 1) = sx * sz / sy; Lambda(2, 2) = sx * sy / sz;
			epsr = epsr * Lambda;
			invMur = (Mur * Lambda).inverse();
		}
		else
		{
			invMur = Mur.inverse();
		}

		//mapping 
		VectorXi MappingIndex = VectorXi::Zero(dof);
		for (int i = 0; i < dof; i++)
			MappingIndex(i) = _mMesh->GetEdgeOfPrism(n, i);

		//submatrix
		MatrixXcd Se, Te;
		Se = MatrixXcd::Zero(dof, dof);
		Te = MatrixXcd::Zero(dof, dof);
		for (int i = 0; i < dof; i++)
			for (int j = 0; j < dof; j++)
				for (int k = 0; k < NbrGuassPoints; k++)
				{
					Se(i, j) = Se(i, j) + wght[k] * DetJac * curlE[k][i].dot(invMur * curlE[k][j]);
					Te(i, j) = Te(i, j) + wght[k] * DetJac * k0 * k0 * E[k][i].dot(epsr * E[k][j]);
				}


		//store in triplet
		if (_mIsReal)
		{
			for (int i = 0; i < dof; i++)
			{
				for (int j = 0; j < dof; j++)
				{
					_mTripleA_real.push_back(Triplet<double>(MappingIndex(i), MappingIndex(j), Se(i, j).real() - Te(i, j).real()));
				}
			}

		}
		else
		{
			for (int i = 0; i < dof; i++)
			{
				for (int j = 0; j < dof; j++)
				{
					_mTripleA_complex.push_back(Triplet<complex<double>>(MappingIndex(i), MappingIndex(j), Se(i, j) - Te(i, j)));
				}
			}

		}
	}

	delete[] u, v, w, wght;

	delete[] vertex;

	for (int i = 0; i < NbrGuassPoints; i++)
	{
		delete[] E[i], curlE[i];
	}
	delete[] E, curlE;
}