372 lines
10 KiB
C++
372 lines
10 KiB
C++
#include"Assemble_Base.h"
|
||
#include"../common/util.h"
|
||
#include<iostream>
|
||
|
||
void OpticsFEM_2D_EigenFreq::Assemble_PEC_PBC()
|
||
{
|
||
//find PEC index
|
||
Eigen::VectorXi PECDomain = Eigen::VectorXi::Zero(_mPhy->GetNbrPEC());
|
||
for (int i = 0; i < _mPhy->GetNbrPEC(); i++)
|
||
PECDomain(i) = _mPhy->GetPECDomain(i) + 1;
|
||
Eigen::VectorXi nodeIndexOfPEC, edgeIndexOfPEC;
|
||
_mMesh->GetIndexOfDomain(PECDomain, nodeIndexOfPEC, edgeIndexOfPEC);
|
||
|
||
//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);
|
||
//行 1:DOF
|
||
//列 1：(DOF-NbrPEC-NbrDst)
|
||
int NbrVertex = _mMesh->GetNbrVertex();
|
||
int NbrDstIndex = nodeIndexOfPEC.rows() + edgeIndexOfPEC.rows() + dstNodeIndex.rows() + dstEdgeIndex.rows();
|
||
//电场类别去除数目
|
||
int numEle = 0;
|
||
if (_mElectricType == 0)
|
||
{
|
||
numEle = _mMesh->GetNbrVertex();
|
||
}
|
||
else if (_mElectricType == 1)
|
||
{
|
||
numEle = _mMesh->GetNbrEdge();
|
||
}
|
||
NbrDstIndex = NbrDstIndex + numEle;
|
||
Eigen::VectorXi dstIndex = Eigen::VectorXi::Zero(NbrDstIndex);
|
||
Eigen::VectorXi srcIndex = Eigen::VectorXi::Zero(NbrDstIndex);
|
||
Eigen::VectorXcd IndexPhi = Eigen::VectorXcd::Zero(NbrDstIndex);
|
||
//电场类别去除索引添加
|
||
if (_mElectricType == 0)
|
||
{
|
||
for (int i = 0; i < numEle; i++)
|
||
{
|
||
dstIndex(i) = i;
|
||
srcIndex(i) = i;
|
||
IndexPhi(i) = 0.;
|
||
}
|
||
}
|
||
else if (_mElectricType == 1)
|
||
{
|
||
for (int i = 0; i < numEle; i++)
|
||
{
|
||
dstIndex(i) = i + _mMesh->GetNbrVertex();
|
||
srcIndex(i) = i + _mMesh->GetNbrVertex();
|
||
IndexPhi(i) = 0.;
|
||
}
|
||
}
|
||
for (int i = 0; i < nodeIndexOfPEC.rows(); i++)
|
||
{
|
||
int num = numEle;
|
||
dstIndex(i + num) = nodeIndexOfPEC(i);
|
||
srcIndex(i + num) = nodeIndexOfPEC(i);
|
||
IndexPhi(i + num) = 0;
|
||
}
|
||
for (int i = 0; i < edgeIndexOfPEC.rows(); i++)
|
||
{
|
||
int num = nodeIndexOfPEC.rows() + numEle;
|
||
dstIndex(i + num) = edgeIndexOfPEC(i) + NbrVertex;
|
||
srcIndex(i + num) = edgeIndexOfPEC(i) + NbrVertex;
|
||
IndexPhi(i + num) = 0.;
|
||
}
|
||
for (int i = 0; i < dstNodeIndex.rows(); i++)
|
||
{
|
||
int num = nodeIndexOfPEC.rows() + edgeIndexOfPEC.rows() + numEle;
|
||
dstIndex(i + num) = dstNodeIndex(i);
|
||
srcIndex(i + num) = srcNodeIndex(i);
|
||
IndexPhi(i + num) = nodePhi(i);
|
||
}
|
||
for (int i = 0; i < dstEdgeIndex.rows(); i++)
|
||
{
|
||
int num = nodeIndexOfPEC.rows() + edgeIndexOfPEC.rows() + dstNodeIndex.rows() + numEle;
|
||
dstIndex(i + num) = dstEdgeIndex(i) + NbrVertex;
|
||
srcIndex(i + num) = srcEdgeIndex(i) + NbrVertex;
|
||
IndexPhi(i + num) = edgePhi(i);
|
||
}
|
||
//排序
|
||
Eigen::VectorXi tempIndex = Eigen::VectorXi::Zero(NbrDstIndex);
|
||
for (int i = 0; i < NbrDstIndex; i++)
|
||
tempIndex(i) = i;
|
||
QuickSort(dstIndex, tempIndex, 0, NbrDstIndex - 1);
|
||
//去重
|
||
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;
|
||
|
||
//不同数目
|
||
int dof = _mMesh->GetNbrEdge() + _mMesh->GetNbrVertex();
|
||
int vdof = _mMesh->GetNbrEdge();
|
||
int sdof = _mMesh->GetNbrVertex();
|
||
|
||
//P矩阵
|
||
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()) //处理PEC或PBC
|
||
{
|
||
if (i == dstIndex(num))
|
||
{
|
||
temp_P.push_back(Eigen::Triplet<double>(dstIndex(num), srcIndex(num), IndexPhi(num).real()));
|
||
num++;
|
||
}
|
||
}
|
||
}
|
||
/*临时P矩阵排序*/
|
||
std::sort(temp_P.begin(), temp_P.end(),
|
||
[](const Eigen::Triplet<double>& a, const Eigen::Triplet<double>& b)
|
||
{
|
||
return a.col() < b.col();
|
||
});
|
||
///*排序后打印*/
|
||
//std::sort(temp_P.begin(), temp_P.end(),
|
||
// [](const Eigen::Triplet<double>& a, const Eigen::Triplet<double>& b) {
|
||
// return a.col() < b.col();
|
||
// });
|
||
// for (const auto& triplet : temp_P) {
|
||
// std::cout << "(" << triplet.row() << ", " << triplet.col() << "): " << triplet.value() << std::endl;
|
||
// }
|
||
|
||
/*删除dstIndex对应列*/
|
||
std::vector<Eigen::Triplet<double>> tripleP_real;
|
||
num = 0;// 用于跟踪新矩阵的列索引
|
||
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());
|
||
|
||
/*排序后打印*/
|
||
//std::sort(tripleP_real.begin(), tripleP_real.end(),
|
||
// [](const Eigen::Triplet<double>& a, const Eigen::Triplet<double>& b) {
|
||
// return a.col() < b.col();
|
||
// });
|
||
// for (const auto& triplet : tripleP_real) {
|
||
// std::cout << "(" << triplet.row() << ", " << triplet.col() << "): " << triplet.value() << std::endl;
|
||
// }
|
||
}
|
||
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()) //处理PEC或PBC
|
||
{
|
||
if (i == dstIndex(num))
|
||
{
|
||
temp_P_complex.push_back(Eigen::Triplet<complex<double>>(dstIndex(num), srcIndex(num), IndexPhi(num)));
|
||
num++;
|
||
}
|
||
}
|
||
}
|
||
/*临时P矩阵排序*/
|
||
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();
|
||
});
|
||
|
||
|
||
/*删除dstIndex对应列*/
|
||
std::vector<Eigen::Triplet<complex<double>>> tripleP_complex;
|
||
num = 0;// 用于跟踪新矩阵的列索引
|
||
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());;
|
||
///*排序后打印*/
|
||
//std::sort(tripleP_complex.begin(), tripleP_complex.end(),
|
||
//[](const Eigen::Triplet<complex<double>>& a, const Eigen::Triplet<complex<double>>& b)
|
||
// {
|
||
// return a.col() < b.col();
|
||
// });
|
||
//for (const auto& triplet : tripleP_complex) {
|
||
// std::cout << "(" << triplet.row() << ", " << triplet.col() << "): " << triplet.value() << std::endl;
|
||
//}
|
||
}
|
||
}
|
||
|
||
|
||
void OpticsFEM_3D_EigenFreq::Assemble_PEC_PBC()
|
||
{
|
||
//find PEC index
|
||
Eigen::VectorXi PECDomain = Eigen::VectorXi::Zero(_mPhy->GetNbrPEC());
|
||
for (int i = 0; i < _mPhy->GetNbrPEC(); i++)
|
||
PECDomain(i) = _mPhy->GetPECDomain(i);
|
||
Eigen::VectorXi edgeIndexOfPEC;
|
||
_mMesh->GetTriIndexOfDomain(PECDomain, edgeIndexOfPEC);
|
||
|
||
//find PBC index
|
||
Eigen::VectorXi srcEdgeIndex, dstEdgeIndex;
|
||
Eigen::VectorXcd 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, srcEdgeIndex, dstEdgeIndex, edgePhi);
|
||
|
||
//行 1:DOF
|
||
//列 1：(DOF-NbrPEC-NbrDst)
|
||
int NbrDstIndex = edgeIndexOfPEC.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 < edgeIndexOfPEC.rows(); i++)
|
||
{
|
||
dstIndex(i) = edgeIndexOfPEC(i);
|
||
srcIndex(i) = edgeIndexOfPEC(i);
|
||
IndexPhi(i) = 0;
|
||
}
|
||
for (int i = 0; i < dstEdgeIndex.rows(); i++)
|
||
{
|
||
int num = edgeIndexOfPEC.rows();
|
||
dstIndex(i + num) = dstEdgeIndex(i);
|
||
srcIndex(i + num) = srcEdgeIndex(i);
|
||
IndexPhi(i + num) = edgePhi(i);
|
||
}
|
||
//排序
|
||
Eigen::VectorXi tempIndex = Eigen::VectorXi::Zero(NbrDstIndex);
|
||
for (int i = 0; i < NbrDstIndex; i++)
|
||
tempIndex(i) = i;
|
||
QuickSort(dstIndex, tempIndex, 0, NbrDstIndex - 1);
|
||
Eigen::VectorXcd tempIndexPhi = Eigen::VectorXcd::Zero(NbrDstIndex);
|
||
for (int i = 0; i < NbrDstIndex; i++)
|
||
tempIndexPhi(i) = IndexPhi(tempIndex(i));
|
||
IndexPhi = tempIndexPhi;
|
||
Eigen::VectorXi tempSrcIndex = Eigen::VectorXi::Zero(NbrDstIndex);
|
||
for (int i = 0; i < NbrDstIndex; i++)
|
||
tempSrcIndex(i) = srcIndex(tempIndex(i));
|
||
srcIndex = tempSrcIndex;
|
||
//去重
|
||
Unique(dstIndex, tempIndex);
|
||
if (NbrDstIndex > tempIndex.rows())
|
||
{
|
||
for (int i = 0; i < tempIndex.rows(); i++)
|
||
IndexPhi(i) = IndexPhi(tempIndex(i));
|
||
IndexPhi = IndexPhi.head(tempIndex.rows());
|
||
|
||
for (int i = 0; i < tempIndex.rows(); i++)
|
||
srcIndex(i) = srcIndex(tempIndex(i));
|
||
srcIndex = srcIndex.head(tempIndex.rows());
|
||
}
|
||
|
||
//P矩阵
|
||
int num = 0;
|
||
if (_mIsReal)
|
||
{
|
||
std::vector<Eigen::Triplet<double>> tripleP_real;
|
||
for (int i = 0; i < _mDof; i++) // 1 4 5
|
||
{
|
||
if (i == dstIndex(num))
|
||
{
|
||
tripleP_real.push_back(Eigen::Triplet<double>(dstIndex(num), srcIndex(num) - num, IndexPhi(num).real()));
|
||
num++;
|
||
}
|
||
else
|
||
{
|
||
tripleP_real.push_back(Eigen::Triplet<double>(i, i - num, 1));
|
||
}
|
||
|
||
}
|
||
|
||
_mP_real.setFromTriplets(tripleP_real.begin(), tripleP_real.end());
|
||
}
|
||
else
|
||
{
|
||
std::vector<Eigen::Triplet<complex<double>>> tripleP_complex;
|
||
for (int i = 0; i < _mDof; i++) // 1 4 5
|
||
{
|
||
if (i == dstIndex(num))
|
||
{
|
||
tripleP_complex.push_back(Eigen::Triplet<complex<double>>(dstIndex(num), srcIndex(num) - num, IndexPhi(num)));
|
||
num++;
|
||
}
|
||
else
|
||
{
|
||
tripleP_complex.push_back(Eigen::Triplet<complex<double>>(i, i - num, 1));
|
||
}
|
||
|
||
}
|
||
|
||
_mP_complex.setFromTriplets(tripleP_complex.begin(), tripleP_complex.end());
|
||
}
|
||
|
||
}
|