yuhui
/
XIAN-FEM-2026June
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
XIAN-FEM-2026June/opticsfem-master/kernel/Assemble_Scatter_Equation.cpp

573 lines
18 KiB
C++
Raw Permalink Blame History

#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);
//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ez<45><7A><EFBFBD><EFBFBD>
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
Jac(0, 0) = vertex[0](0) - vertex[3](0); Jac(0, 1) = vertex[0](1) - vertex[3](1); Jac(0, 2) = 0.;
Jac(1, 0) = vertex[1](0) - vertex[3](0); Jac(1, 1) = vertex[1](1) - vertex[3](1); Jac(1, 2) = 0.;
Jac(2, 0) = vertex[2](0) - vertex[3](0); Jac(2, 1) = vertex[2](1) - vertex[3](1); Jac(2, 2) = 1.;
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->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);
//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->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);
//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;
}
Reference in New Issue View Git Blame Copy Permalink