596 lines
19 KiB
C++
596 lines
19 KiB
C++
#include "PortModeSolver.h"
|
||
#include "../mesh/Mesh_Base.h"
|
||
#include "../material/Material_Base.h"
|
||
#include "../function/BF.h"
|
||
#include "../function/Gauss.h"
|
||
#include "../common/define.h"
|
||
#include "../nlohmann/json.hpp"
|
||
#include "../Eigen/Sparse"
|
||
#include "../Eigen/SparseLU"
|
||
#include "../Eigen/Eigenvalues"
|
||
|
||
#include <algorithm>
|
||
#include <array>
|
||
#include <cmath>
|
||
#include <complex>
|
||
#include <fstream>
|
||
#include <iostream>
|
||
#include <map>
|
||
#include <set>
|
||
#include <vector>
|
||
|
||
using namespace Eigen;
|
||
using json = nlohmann::json;
|
||
|
||
namespace {
|
||
|
||
struct PortBoundaryMesh
|
||
{
|
||
MatrixXd nodes; // nNodes x 3
|
||
MatrixXi newFaces; // nf x 3 (0-based local node)
|
||
MatrixXi facesConn; // nf x 2 (1-based tet, 1-based local face)
|
||
VectorXi volDomain; // nf (1-based volume domain for material)
|
||
MatrixXi edgeOfFace; // nf x 3 (0-based edge DOF)
|
||
MatrixXi edgeNodes; // nEdges x 2 (0-based local nodes)
|
||
int nNodes = 0;
|
||
int nEdges = 0;
|
||
int nFaces = 0;
|
||
};
|
||
|
||
bool domainInList(int domain, const VectorXi& list)
|
||
{
|
||
for (int i = 0; i < list.size(); ++i)
|
||
if (list(i) == domain)
|
||
return true;
|
||
return false;
|
||
}
|
||
|
||
bool extractPortMesh(Mesh_3D* mesh, const VectorXi& domains, PortBoundaryMesh& out)
|
||
{
|
||
std::vector<int> triIds;
|
||
const int nTri = mesh->GetNbrTri();
|
||
for (int t = 0; t < nTri; ++t)
|
||
{
|
||
if (domainInList(mesh->GetDomainOfTri(t), domains))
|
||
triIds.push_back(t);
|
||
}
|
||
if (triIds.empty())
|
||
{
|
||
std::cerr << "[PortModeSolver] no port triangles for domains" << std::endl;
|
||
return false;
|
||
}
|
||
|
||
const int nf = static_cast<int>(triIds.size());
|
||
MatrixXi oldFaces(nf, 3);
|
||
out.facesConn = MatrixXi::Zero(nf, 2);
|
||
out.volDomain = VectorXi::Zero(nf);
|
||
|
||
std::vector<int> globalNodes;
|
||
globalNodes.reserve(nf * 3);
|
||
for (int i = 0; i < nf; ++i)
|
||
{
|
||
const int tri = triIds[i];
|
||
for (int k = 0; k < 3; ++k)
|
||
{
|
||
oldFaces(i, k) = mesh->GetTri(tri, k);
|
||
globalNodes.push_back(oldFaces(i, k));
|
||
}
|
||
Vector3i conn;
|
||
mesh->GetCoonOfTri(tri, conn);
|
||
out.facesConn(i, 0) = conn(0) + 1;
|
||
out.facesConn(i, 1) = conn(1) + 1;
|
||
out.volDomain(i) = mesh->GetDomainOfTet(conn(0));
|
||
}
|
||
|
||
std::sort(globalNodes.begin(), globalNodes.end());
|
||
globalNodes.erase(std::unique(globalNodes.begin(), globalNodes.end()), globalNodes.end());
|
||
out.nNodes = static_cast<int>(globalNodes.size());
|
||
out.nodes = MatrixXd::Zero(out.nNodes, 3);
|
||
std::map<int, int> g2l;
|
||
for (int i = 0; i < out.nNodes; ++i)
|
||
{
|
||
g2l[globalNodes[i]] = i;
|
||
Vector3d v;
|
||
mesh->GetVertex(globalNodes[i], v);
|
||
out.nodes.row(i) = v.transpose();
|
||
}
|
||
|
||
out.newFaces = MatrixXi::Zero(nf, 3);
|
||
for (int i = 0; i < nf; ++i)
|
||
for (int k = 0; k < 3; ++k)
|
||
out.newFaces(i, k) = g2l[oldFaces(i, k)];
|
||
|
||
// Local edges like MATLAB: (n1,n2),(n1,n3),(n2,n3) without sorting pair; unique by lexicographic rows.
|
||
std::vector<std::array<int, 2>> edgeRows;
|
||
edgeRows.reserve(static_cast<size_t>(nf) * 3);
|
||
for (int i = 0; i < nf; ++i)
|
||
{
|
||
const int a = out.newFaces(i, 0);
|
||
const int b = out.newFaces(i, 1);
|
||
const int c = out.newFaces(i, 2);
|
||
edgeRows.push_back({ a, b });
|
||
edgeRows.push_back({ a, c });
|
||
edgeRows.push_back({ b, c });
|
||
}
|
||
std::vector<std::array<int, 2>> uniqueEdges = edgeRows;
|
||
std::sort(uniqueEdges.begin(), uniqueEdges.end());
|
||
uniqueEdges.erase(std::unique(uniqueEdges.begin(), uniqueEdges.end()), uniqueEdges.end());
|
||
out.nEdges = static_cast<int>(uniqueEdges.size());
|
||
std::map<std::array<int, 2>, int> e2i;
|
||
for (int i = 0; i < out.nEdges; ++i)
|
||
e2i[uniqueEdges[i]] = i;
|
||
|
||
out.edgeOfFace = MatrixXi::Zero(nf, 3);
|
||
out.nFaces = nf;
|
||
for (int i = 0; i < nf; ++i)
|
||
{
|
||
for (int k = 0; k < 3; ++k)
|
||
out.edgeOfFace(i, k) = e2i[edgeRows[static_cast<size_t>(i) * 3 + k]];
|
||
}
|
||
out.edgeNodes = MatrixXi::Zero(out.nEdges, 2);
|
||
for (int i = 0; i < out.nEdges; ++i)
|
||
{
|
||
out.edgeNodes(i, 0) = uniqueEdges[static_cast<size_t>(i)][0];
|
||
out.edgeNodes(i, 1) = uniqueEdges[static_cast<size_t>(i)][1];
|
||
}
|
||
return true;
|
||
}
|
||
|
||
void assemblePortEigen(
|
||
const PortBoundaryMesh& pm,
|
||
MaterialLib* matLib,
|
||
double lam0,
|
||
SparseMatrix<std::complex<double>, RowMajor>& Ared,
|
||
SparseMatrix<std::complex<double>, RowMajor>& Bred,
|
||
SparseMatrix<std::complex<double>, RowMajor>& P,
|
||
int& dofFull)
|
||
{
|
||
const double k0 = 2.0 * Pi / lam0;
|
||
const int nNodes = pm.nNodes;
|
||
const int nEdges = pm.nEdges;
|
||
dofFull = nNodes + nEdges;
|
||
|
||
Gauss gauss;
|
||
const int nGP = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 2);
|
||
std::vector<double> u(nGP), v(nGP), w(nGP), wght(nGP);
|
||
gauss.GetGaussPoints(TWODIM, TRIANGLE, u.data(), v.data(), w.data(), wght.data());
|
||
|
||
BF bfEt, bfCurlEt, bfEz, bfGradEz;
|
||
bfEt.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
|
||
bfCurlEt.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
|
||
bfEz.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
|
||
bfGradEz.GetNbrBF(TWODIM, TRIANGLE, BF_GRAD_LAGRANGE, BF_LINEFUNC);
|
||
|
||
using Trip = Triplet<std::complex<double>>;
|
||
std::vector<Trip> tripsA, tripsB;
|
||
tripsA.reserve(static_cast<size_t>(pm.nFaces) * 9);
|
||
tripsB.reserve(static_cast<size_t>(pm.nFaces) * 36);
|
||
|
||
std::vector<int> edgeCount(static_cast<size_t>(nEdges), 0);
|
||
for (int f = 0; f < pm.nFaces; ++f)
|
||
for (int k = 0; k < 3; ++k)
|
||
edgeCount[static_cast<size_t>(pm.edgeOfFace(f, k))]++;
|
||
|
||
for (int n = 0; n < pm.nFaces; ++n)
|
||
{
|
||
double x[3], y[3], l[3];
|
||
for (int i = 0; i < 3; ++i)
|
||
{
|
||
x[i] = pm.nodes(pm.newFaces(n, i), 0);
|
||
y[i] = pm.nodes(pm.newFaces(n, i), 1);
|
||
}
|
||
l[0] = std::sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
|
||
l[1] = std::sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
|
||
l[2] = std::sqrt((x[2] - x[1]) * (x[2] - x[1]) + (y[2] - y[1]) * (y[2] - y[1]));
|
||
|
||
Matrix3d Jac = Matrix3d::Zero();
|
||
Jac(0, 0) = x[1] - x[0]; Jac(0, 1) = y[1] - y[0];
|
||
Jac(1, 0) = x[2] - x[0]; Jac(1, 1) = y[2] - y[0];
|
||
Jac(2, 2) = 1.0;
|
||
const Matrix3d InvJac = Jac.inverse();
|
||
const double DetJac = std::abs(Jac.determinant());
|
||
const Matrix3d TJac = Jac.transpose() / Jac.determinant();
|
||
|
||
std::vector<Vector3d> Et(3 * nGP), curlEt(3 * nGP), Ez(3 * nGP), gradEz(3 * nGP);
|
||
auto at = [nGP](std::vector<Vector3d>& a, int j, int gp) -> Vector3d& {
|
||
return a[static_cast<size_t>(j * nGP + gp)];
|
||
};
|
||
for (int gp = 0; gp < nGP; ++gp)
|
||
{
|
||
for (int j = 0; j < 3; ++j)
|
||
{
|
||
Vector3d temp;
|
||
bfEt.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
|
||
at(Et, j, gp) = InvJac * temp * l[j];
|
||
bfCurlEt.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
|
||
at(curlEt, j, gp) = TJac * temp * l[j];
|
||
bfEz.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
|
||
at(Ez, j, gp) = temp;
|
||
bfGradEz.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
|
||
at(gradEz, j, gp) = InvJac * temp;
|
||
}
|
||
}
|
||
|
||
const int domain = pm.volDomain(n);
|
||
Matrix3cd epsr = matLib->GetEpsr(domain);
|
||
const Matrix3d sigma = matLib->GetSigma(domain);
|
||
epsr = epsr - sigma * std::complex<double>(0.0, 1.0 / k0 * 120.0 * Pi);
|
||
const Matrix3cd mur = matLib->GetMur(domain);
|
||
const std::complex<double> mu = mur(0, 0);
|
||
const std::complex<double> epsilon = epsr(0, 0);
|
||
|
||
Matrix3cd St = Matrix3cd::Zero(), Tte = Matrix3cd::Zero();
|
||
Matrix3cd Sz = Matrix3cd::Zero(), Tz = Matrix3cd::Zero();
|
||
Matrix3cd G = Matrix3cd::Zero(), Ttu = Matrix3cd::Zero();
|
||
for (int i = 0; i < 3; ++i)
|
||
{
|
||
for (int j = 0; j < 3; ++j)
|
||
{
|
||
for (int gp = 0; gp < nGP; ++gp)
|
||
{
|
||
const double wgt = wght[gp] * DetJac;
|
||
St(i, j) += wgt / mu * at(curlEt, i, gp).dot(at(curlEt, j, gp));
|
||
Tte(i, j) += wgt * k0 * k0 * epsilon * at(Et, i, gp).dot(at(Et, j, gp));
|
||
Sz(i, j) += wgt / mu * at(gradEz, i, gp).dot(at(gradEz, j, gp));
|
||
Tz(i, j) += wgt * k0 * k0 * epsilon * at(Ez, i, gp).dot(at(Ez, j, gp));
|
||
G(i, j) += wgt / mu * at(Et, i, gp).dot(at(gradEz, j, gp));
|
||
Ttu(i, j) += wgt / mu * at(Et, i, gp).dot(at(Et, j, gp));
|
||
}
|
||
}
|
||
}
|
||
const Matrix3cd Gt = G.transpose();
|
||
|
||
for (int i = 0; i < 3; ++i)
|
||
{
|
||
for (int j = 0; j < 3; ++j)
|
||
{
|
||
const int si = pm.newFaces(n, i);
|
||
const int sj = pm.newFaces(n, j);
|
||
const int vi = pm.edgeOfFace(n, i) + nNodes;
|
||
const int vj = pm.edgeOfFace(n, j) + nNodes;
|
||
tripsA.emplace_back(vi, vj, St(i, j) - Tte(i, j));
|
||
tripsB.emplace_back(vi, vj, Ttu(i, j));
|
||
tripsB.emplace_back(si, sj, Sz(i, j) - Tz(i, j));
|
||
tripsB.emplace_back(vi, sj, G(i, j));
|
||
tripsB.emplace_back(si, vj, Gt(i, j));
|
||
}
|
||
}
|
||
}
|
||
|
||
SparseMatrix<std::complex<double>, RowMajor> Afull(dofFull, dofFull), Bfull(dofFull, dofFull);
|
||
Afull.setFromTriplets(tripsA.begin(), tripsA.end());
|
||
Bfull.setFromTriplets(tripsB.begin(), tripsB.end());
|
||
|
||
// Auto PEC: boundary edges (count==1) and their nodes (MATLAB EigenModeMatrixAssembly)
|
||
std::set<int> pec;
|
||
for (int e = 0; e < nEdges; ++e)
|
||
{
|
||
if (edgeCount[static_cast<size_t>(e)] != 1)
|
||
continue;
|
||
pec.insert(pm.edgeNodes(e, 0));
|
||
pec.insert(pm.edgeNodes(e, 1));
|
||
pec.insert(e + nNodes);
|
||
}
|
||
|
||
std::vector<int> pecInd(pec.begin(), pec.end());
|
||
std::sort(pecInd.begin(), pecInd.end());
|
||
std::vector<char> isPec(static_cast<size_t>(dofFull), 0);
|
||
for (int idx : pecInd)
|
||
if (idx >= 0 && idx < dofFull)
|
||
isPec[static_cast<size_t>(idx)] = 1;
|
||
|
||
std::vector<int> freeOfFull;
|
||
freeOfFull.reserve(static_cast<size_t>(dofFull));
|
||
for (int i = 0; i < dofFull; ++i)
|
||
if (!isPec[static_cast<size_t>(i)])
|
||
freeOfFull.push_back(i);
|
||
const int nFree = static_cast<int>(freeOfFull.size());
|
||
|
||
std::vector<Trip> pTrips;
|
||
pTrips.reserve(static_cast<size_t>(nFree));
|
||
for (int j = 0; j < nFree; ++j)
|
||
pTrips.emplace_back(freeOfFull[j], j, 1.0);
|
||
P = SparseMatrix<std::complex<double>, RowMajor>(dofFull, nFree);
|
||
P.setFromTriplets(pTrips.begin(), pTrips.end());
|
||
|
||
Ared = SparseMatrix<std::complex<double>, RowMajor>(nFree, nFree);
|
||
Bred = SparseMatrix<std::complex<double>, RowMajor>(nFree, nFree);
|
||
Ared = SparseMatrix<std::complex<double>, RowMajor>(P.transpose() * Afull * P);
|
||
Bred = SparseMatrix<std::complex<double>, RowMajor>(P.transpose() * Bfull * P);
|
||
}
|
||
|
||
double computePowerCoef(
|
||
const PortBoundaryMesh& pm,
|
||
MaterialLib* matLib,
|
||
double lam0,
|
||
const std::complex<double>& gamma,
|
||
const VectorXcd& Ez,
|
||
const VectorXcd& Et)
|
||
{
|
||
const double k0 = 2.0 * Pi / lam0;
|
||
Gauss gauss;
|
||
const int nGP = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 2);
|
||
std::vector<double> u(nGP), v(nGP), w(nGP), wght(nGP);
|
||
gauss.GetGaussPoints(TWODIM, TRIANGLE, u.data(), v.data(), w.data(), wght.data());
|
||
|
||
BF bfEt, bfCurlEt, bfEz, bfCurlEz;
|
||
bfEt.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
|
||
bfCurlEt.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
|
||
bfEz.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
|
||
bfCurlEz.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_LAGRANGE, BF_LINEFUNC);
|
||
|
||
double power = 0.0;
|
||
for (int n = 0; n < pm.nFaces; ++n)
|
||
{
|
||
double x[3], y[3], l[3];
|
||
for (int i = 0; i < 3; ++i)
|
||
{
|
||
x[i] = pm.nodes(pm.newFaces(n, i), 0);
|
||
y[i] = pm.nodes(pm.newFaces(n, i), 1);
|
||
}
|
||
l[0] = std::sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
|
||
l[1] = std::sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
|
||
l[2] = std::sqrt((x[2] - x[1]) * (x[2] - x[1]) + (y[2] - y[1]) * (y[2] - y[1]));
|
||
|
||
Matrix3d Jac = Matrix3d::Zero();
|
||
Jac(0, 0) = x[1] - x[0]; Jac(0, 1) = y[1] - y[0];
|
||
Jac(1, 0) = x[2] - x[0]; Jac(1, 1) = y[2] - y[0];
|
||
Jac(2, 2) = 1.0;
|
||
const Matrix3d InvJac = Jac.inverse();
|
||
const double DetJac = std::abs(Jac.determinant());
|
||
const Matrix3d TJac = Jac.transpose() / Jac.determinant();
|
||
Matrix3d JacS = Matrix3d::Zero();
|
||
JacS(0, 0) = InvJac(1, 1); JacS(0, 1) = -InvJac(1, 0);
|
||
JacS(1, 0) = -InvJac(0, 1); JacS(1, 1) = InvJac(0, 0);
|
||
JacS(2, 2) = 1.0;
|
||
|
||
const int domain = pm.volDomain(n);
|
||
const double mu = std::real(matLib->GetMur(domain)(0, 0));
|
||
const std::complex<double> physicsIndex = std::complex<double>(0.0, 1.0) / (k0 * 120.0 * Pi) / mu;
|
||
|
||
for (int gp = 0; gp < nGP; ++gp)
|
||
{
|
||
Vector3cd EE = Vector3cd::Zero();
|
||
Vector3cd curlEE = Vector3cd::Zero();
|
||
for (int i = 0; i < 3; ++i)
|
||
{
|
||
Vector3d temp;
|
||
bfEt.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
|
||
const Vector3d et = InvJac * temp * l[i];
|
||
bfEz.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
|
||
const Vector3d ezVec(0, 0, temp(2));
|
||
bfCurlEt.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
|
||
const Vector3d curlEtVec = TJac * Vector3d(0, 0, temp(2)) * l[i];
|
||
bfCurlEz.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
|
||
const Vector3d curlEzVec = JacS * Vector3d(temp(0), temp(1), 0.0);
|
||
|
||
const int mapEt = pm.edgeOfFace(n, i);
|
||
const int mapEz = pm.newFaces(n, i);
|
||
EE += et.cast<std::complex<double>>() * Et(mapEt) + ezVec.cast<std::complex<double>>() * Ez(mapEz);
|
||
curlEE += curlEtVec.cast<std::complex<double>>() * Et(mapEt)
|
||
+ curlEzVec.cast<std::complex<double>>() * Ez(mapEz);
|
||
}
|
||
Vector3cd HH = curlEE;
|
||
HH(0) += EE(1) * gamma;
|
||
HH(1) -= EE(0) * gamma;
|
||
HH *= physicsIndex;
|
||
|
||
const double rEx = EE(0).real(), iEx = EE(0).imag();
|
||
const double rEy = EE(1).real(), iEy = EE(1).imag();
|
||
const double rHx = HH(0).real(), iHx = HH(0).imag();
|
||
const double rHy = HH(1).real(), iHy = HH(1).imag();
|
||
power += wght[gp] * 0.5 * DetJac * (-rEy * rHx - iEy * iHx + rEx * rHy + iEx * iHy);
|
||
}
|
||
}
|
||
if (!(power > 0.0))
|
||
{
|
||
std::cerr << "[PortModeSolver] power<=0 (" << power << "), powerCoef=1" << std::endl;
|
||
return 1.0;
|
||
}
|
||
return std::sqrt(1.0 / power);
|
||
}
|
||
|
||
bool solveOnePort(
|
||
Mesh_3D* mesh,
|
||
MaterialLib* matLib,
|
||
const VectorXi& domains,
|
||
double lambda0,
|
||
double targetNeff,
|
||
int modeNum,
|
||
int portType,
|
||
bool normalizePower,
|
||
PortNumericMode& mode)
|
||
{
|
||
PortBoundaryMesh pm;
|
||
if (!extractPortMesh(mesh, domains, pm))
|
||
return false;
|
||
|
||
SparseMatrix<std::complex<double>, RowMajor> A, B, P;
|
||
int dofFull = 0;
|
||
assemblePortEigen(pm, matLib, lambda0, A, B, P, dofFull);
|
||
|
||
const double k0 = 2.0 * Pi / lambda0;
|
||
// MATLAB: sigma = -(k0*targetNeff)^2, A x = λ B x, gamma=sqrt(λ), neff=(-i*gamma)/k0
|
||
const std::complex<double> sigma(-(k0 * targetNeff) * (k0 * targetNeff), 0.0);
|
||
|
||
const int nFree = static_cast<int>(A.rows());
|
||
if (nFree <= 0)
|
||
{
|
||
std::cerr << "[PortModeSolver] empty free DOF after PEC" << std::endl;
|
||
return false;
|
||
}
|
||
const MatrixXcd Ad = MatrixXcd(A);
|
||
const MatrixXcd Bd = MatrixXcd(B);
|
||
const MatrixXcd As = Ad - sigma * Bd;
|
||
Eigen::PartialPivLU<MatrixXcd> lu(As);
|
||
if (lu.determinant() == std::complex<double>(0, 0))
|
||
{
|
||
std::cerr << "[PortModeSolver] (A-σB) singular" << std::endl;
|
||
return false;
|
||
}
|
||
// M = (A-σB)^{-1} B ; eig(M)=ν=1/(λ-σ) ⇒ λ = σ + 1/ν
|
||
const MatrixXcd M = lu.solve(Bd);
|
||
Eigen::ComplexEigenSolver<MatrixXcd> ces(M);
|
||
if (ces.info() != Eigen::Success)
|
||
{
|
||
std::cerr << "[PortModeSolver] ComplexEigenSolver failed" << std::endl;
|
||
return false;
|
||
}
|
||
const VectorXcd nu = ces.eigenvalues();
|
||
const MatrixXcd V = ces.eigenvectors();
|
||
// Prefer largest |ν| (closest to σ), then pick among top few by Re(neff)~target
|
||
std::vector<int> order(static_cast<size_t>(nu.size()));
|
||
for (int i = 0; i < nu.size(); ++i)
|
||
order[static_cast<size_t>(i)] = i;
|
||
std::sort(order.begin(), order.end(), [&](int a, int b) {
|
||
return std::abs(nu(a)) > std::abs(nu(b));
|
||
});
|
||
int best = order[0];
|
||
double bestScore = 1e300;
|
||
const int nCand = std::min(modeNum + 8, static_cast<int>(order.size()));
|
||
for (int c = 0; c < nCand; ++c)
|
||
{
|
||
const int i = order[static_cast<size_t>(c)];
|
||
if (std::abs(nu(i)) < 1e-30)
|
||
continue;
|
||
const std::complex<double> lam = sigma + std::complex<double>(1.0, 0.0) / nu(i);
|
||
const std::complex<double> gam = std::sqrt(lam);
|
||
const std::complex<double> ne = -std::complex<double>(0.0, 1.0) * gam / k0;
|
||
const double score = std::abs(ne.real() - targetNeff) + 0.25 * std::abs(ne.imag());
|
||
if (score < bestScore)
|
||
{
|
||
bestScore = score;
|
||
best = i;
|
||
}
|
||
}
|
||
const std::complex<double> lamBest = sigma + std::complex<double>(1.0, 0.0) / nu(best);
|
||
std::complex<double> gamma = std::sqrt(lamBest);
|
||
// Branch: prefer Im(gamma)>0 for decaying envelope convention used with neff≈-i*gamma/k0
|
||
std::complex<double> neff = -std::complex<double>(0.0, 1.0) * gamma / k0;
|
||
if (neff.real() < 0.0)
|
||
{
|
||
gamma = -gamma;
|
||
neff = -std::complex<double>(0.0, 1.0) * gamma / k0;
|
||
}
|
||
|
||
VectorXcd xFree = V.col(best);
|
||
VectorXcd xFull = P * xFree;
|
||
|
||
std::cout << "[PortModeSolver] type=" << portType
|
||
<< " neff=" << neff << " gamma=" << gamma
|
||
<< " nFaces=" << pm.nFaces << " nNodes=" << pm.nNodes
|
||
<< " nEdges=" << pm.nEdges << " nFree=" << nFree << std::endl;
|
||
|
||
VectorXcd Ez = xFull.head(pm.nNodes);
|
||
VectorXcd Et = xFull.segment(pm.nNodes, pm.nEdges);
|
||
if (std::abs(gamma) > 0.0)
|
||
Et /= gamma;
|
||
|
||
// Phase gauge e^{iφ}=1: make the peak-|Et| DOF real and positive (COMSOL-like Mode phase 0).
|
||
{
|
||
Eigen::Index iPeak = 0;
|
||
Et.cwiseAbs().maxCoeff(&iPeak);
|
||
const std::complex<double> peak = Et(iPeak);
|
||
if (std::abs(peak) > 0.0)
|
||
{
|
||
const std::complex<double> fix = std::abs(peak) / peak;
|
||
Et *= fix;
|
||
Ez *= fix;
|
||
std::cout << "[PortModeSolver] phase gauge: arg(Et_peak)->0 (e^{iφ}=1)"
|
||
<< " iPeak=" << iPeak << std::endl;
|
||
}
|
||
}
|
||
|
||
double powerCoef = 1.0;
|
||
if (normalizePower && portType == 1)
|
||
powerCoef = computePowerCoef(pm, matLib, lambda0, gamma, Ez, Et);
|
||
|
||
mode.type = portType;
|
||
mode.gamma = gamma;
|
||
mode.powerCoef = powerCoef;
|
||
mode.normal = (portType == 1) ? Vector3d(0, 0, -1) : Vector3d(0, 0, 1);
|
||
mode.facesConn = pm.facesConn;
|
||
mode.portNewFaces = pm.newFaces;
|
||
mode.portEdgeOfFace = pm.edgeOfFace;
|
||
mode.portNodes = pm.nodes;
|
||
mode.Ez = Ez;
|
||
mode.Et = Et;
|
||
|
||
std::cout << "[PortModeSolver] powerCoef=" << powerCoef
|
||
<< " |Et|_max=" << Et.cwiseAbs().maxCoeff()
|
||
<< " |Ez|_max=" << Ez.cwiseAbs().maxCoeff() << std::endl;
|
||
return true;
|
||
}
|
||
|
||
} // namespace
|
||
|
||
bool ComputePortModesFromMesh(
|
||
Mesh_3D* mesh,
|
||
MaterialLib* matLib,
|
||
const VectorXi& domainsInc,
|
||
const VectorXi& domainsOut,
|
||
double lambda0,
|
||
double targetNeff,
|
||
int modeNum,
|
||
PortModeLibrary& outModes,
|
||
bool normalizeInputPower)
|
||
{
|
||
outModes.Clear();
|
||
PortNumericMode inMode, outMode;
|
||
if (!solveOnePort(mesh, matLib, domainsInc, lambda0, targetNeff, modeNum, 1, normalizeInputPower, inMode))
|
||
return false;
|
||
if (!solveOnePort(mesh, matLib, domainsOut, lambda0, targetNeff, modeNum, 2, false, outMode))
|
||
return false;
|
||
outModes.AddMode(inMode);
|
||
outModes.AddMode(outMode);
|
||
return true;
|
||
}
|
||
|
||
bool SavePortModesToFile(const PortModeLibrary& modes, const std::string& path)
|
||
{
|
||
if (modes.GetNbrModes() < 2)
|
||
return false;
|
||
auto pack = [](const PortNumericMode& m) {
|
||
json j;
|
||
j["type"] = m.type;
|
||
j["gamma"] = { m.gamma.real(), m.gamma.imag() };
|
||
j["powerCoef"] = m.powerCoef;
|
||
j["normal"] = { m.normal(0), m.normal(1), m.normal(2) };
|
||
j["facesConn"] = json::array();
|
||
for (int i = 0; i < m.facesConn.rows(); ++i)
|
||
j["facesConn"].push_back({ m.facesConn(i, 0), m.facesConn(i, 1) });
|
||
j["portNewFaces"] = json::array();
|
||
for (int i = 0; i < m.portNewFaces.rows(); ++i)
|
||
j["portNewFaces"].push_back({ m.portNewFaces(i, 0), m.portNewFaces(i, 1), m.portNewFaces(i, 2) });
|
||
j["portEdgeOfFace"] = json::array();
|
||
for (int i = 0; i < m.portEdgeOfFace.rows(); ++i)
|
||
j["portEdgeOfFace"].push_back({ m.portEdgeOfFace(i, 0), m.portEdgeOfFace(i, 1), m.portEdgeOfFace(i, 2) });
|
||
j["portNodes"] = json::array();
|
||
for (int i = 0; i < m.portNodes.rows(); ++i)
|
||
j["portNodes"].push_back({ m.portNodes(i, 0), m.portNodes(i, 1), m.portNodes(i, 2) });
|
||
j["Ez_re"] = json::array();
|
||
j["Ez_im"] = json::array();
|
||
for (int i = 0; i < m.Ez.size(); ++i)
|
||
{
|
||
j["Ez_re"].push_back(m.Ez(i).real());
|
||
j["Ez_im"].push_back(m.Ez(i).imag());
|
||
}
|
||
j["Et_re"] = json::array();
|
||
j["Et_im"] = json::array();
|
||
for (int i = 0; i < m.Et.size(); ++i)
|
||
{
|
||
j["Et_re"].push_back(m.Et(i).real());
|
||
j["Et_im"].push_back(m.Et(i).imag());
|
||
}
|
||
return j;
|
||
};
|
||
json js;
|
||
js["input"] = pack(modes.GetMode(0));
|
||
js["output"] = pack(modes.GetMode(1));
|
||
std::ofstream ofs(path);
|
||
if (!ofs.is_open())
|
||
return false;
|
||
ofs << js.dump();
|
||
return true;
|
||
}
|