XIAN-FEM-2026June/3D opticsfem-master/mesh/Mesh_Find.cpp

#include"Mesh_Base.h"
#include<algorithm>
#include<complex>
#include<iostream>

//get edges with edge's flag of domain
void Mesh_2D::GetIndexOfDomain(Eigen::VectorXi domain, Eigen::VectorXi& nodeIndex, Eigen::VectorXi& edgeIndex)
{
	int NbrEdges = 0;

	int NbrDomain = domain.rows();
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j) == domain(i))
				NbrEdges++;
		}
	}

	//store edgeIndex and nodeIndex
	Eigen::Vector2i coonOfEdges;
	int* tempEdgeIndex, * tempNodeIndex;
	tempEdgeIndex = new int[NbrEdges];
	tempNodeIndex = new int[NbrEdges * 2];
	int numNode = 0;
	int numEdge = 0;
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j) == domain(i))
			{
				this->GetCoonOfEdges(j, coonOfEdges); //coonOfEdges 0:tri 1:numOfEdge
				tempEdgeIndex[numEdge] = GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1));
				numEdge++;
				tempNodeIndex[numNode] = GetEdge(GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1)), 0);
				numNode++;
				tempNodeIndex[numNode] = GetEdge(GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1)), 1);
				numNode++;
			}
		}
	}



	//sort and unique nodeIndex
	std::sort(tempNodeIndex, tempNodeIndex + NbrEdges * 2);
	int NbrNodes = (std::unique(tempNodeIndex, tempNodeIndex + NbrEdges * 2) - tempNodeIndex);
	nodeIndex = Eigen::VectorXi::Zero(NbrNodes);
	for (int i = 0; i < NbrNodes; i++)
		nodeIndex(i) = tempNodeIndex[i];

	//sort edgeIndex
	edgeIndex = Eigen::VectorXi::Zero(NbrEdges);
	std::sort(tempEdgeIndex, tempEdgeIndex + NbrEdges);
	for (int i = 0; i < NbrEdges; i++)
		edgeIndex(i) = tempEdgeIndex[i];

	delete[] tempEdgeIndex, tempNodeIndex;


}

void Mesh_2D::GetIndexOfDomain2(Eigen::VectorXi domain, Eigen::VectorXi& edgeIndex, Eigen::VectorXi& edgeNum)
{
	int NbrEdges = 0;

	int NbrDomain = domain.rows();
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j) == domain(i))
				NbrEdges++;
		}
	}

	//store edgeIndex
	edgeIndex = Eigen::VectorXi::Zero(NbrEdges);
	edgeNum = Eigen::VectorXi::Zero(NbrEdges);
	int numEdge = 0;
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j) == domain(i))
			{
				edgeIndex(numEdge) = j;
				edgeNum(numEdge) = i;
				numEdge++;
			}
		}
	}

}

void Mesh_2D::GetIndexOfPBC(Eigen::MatrixXi PBC, Eigen::VectorXcd PBCData,
	Eigen::VectorXi& srcNodeIndex, Eigen::VectorXi& dstNodeIndex, Eigen::VectorXcd& nodePhi,
	Eigen::VectorXi& srcEdgeIndex, Eigen::VectorXi& dstEdgeIndex, Eigen::VectorXcd& edgePhi)
{
	int NbrEdges = 0;

	int NbrDomain = PBC.rows();
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j)-1 == PBC(i, 0))
				NbrEdges++;
		}
	}
	//store srcEdgeIndex and srcNodeIndex
	Eigen::Vector2i coonOfEdges;
	int* tempSrcEdgeIndex, * tempSrcNodeIndex;
	tempSrcEdgeIndex = new int[NbrEdges];
	tempSrcNodeIndex = new int[NbrEdges * 2];
	std::complex<double>* tempEdgePhi, * tempNodePhi;
	tempEdgePhi = new std::complex<double>[NbrEdges];
	tempNodePhi = new std::complex<double>[NbrEdges * 2];
	int numNode = 0;
	int numEdge = 0;
	int Tage_numEdge = 0;
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrEdges; j++)
		{
			if (GetDomainOfEdges(j) - 1 == PBC(i, 0))
			{
				tempSrcEdgeIndex[numEdge] = j;
				tempEdgePhi[numEdge] = PBCData(i);
				numEdge++;
				this->GetCoonOfEdges(j, coonOfEdges);
				Tage_numEdge = GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1));
				tempSrcNodeIndex[numNode] = GetEdge(Tage_numEdge, 0);
				tempNodePhi[numNode] = PBCData(i);
				numNode++;
				tempSrcNodeIndex[numNode] = GetEdge(Tage_numEdge, 1);
				tempNodePhi[numNode] = PBCData(i);
				numNode++;
			}
		}
	}
	//store dstEdgeIndex and dstNodeIndex
	int* tempDstEdgeIndex, * tempDstNodeIndex;
	tempDstEdgeIndex = new int[NbrEdges];
	tempDstNodeIndex = new int[NbrEdges * 2];
	for (int i = 0; i < NbrEdges; i++)
	{
		for (int j = 0; j < _mCopyOfEdges.rows(); j++)
		{
			if (tempSrcEdgeIndex[i] == _mCopyOfEdges(j, 0))
			{
				tempDstEdgeIndex[i] = _mCopyOfEdges(j, 1);
				break;
			}

		}
	}
	for (int i = 0; i < NbrEdges; i++)
	{
		this->GetCoonOfEdges(tempDstEdgeIndex[i], coonOfEdges);
		Tage_numEdge = GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1));
		tempDstNodeIndex[i * 2] = GetEdge(Tage_numEdge, 0);
		tempDstNodeIndex[i * 2 + 1] = GetEdge(Tage_numEdge, 1);
	}
	//<2F><><EFBFBD>ݷ<EFBFBD><DDB7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> DstNodeIndex <20><> edgePhi
	//δ<><CEB4><EFBFBD><EFBFBD>rotational bloch bounadry
	for (int i = 0; i < NbrEdges; i++)
	{
		for (int j = 0; j < _mCopyOfVertex.rows(); j++)
		{
			if (tempSrcNodeIndex[2 * i] == _mCopyOfVertex(j, 0) && tempDstNodeIndex[2 * i + 1] == _mCopyOfVertex(j, 1))
			{
				int temp = tempDstNodeIndex[i * 2];
				tempDstNodeIndex[i * 2] = tempDstNodeIndex[i * 2 + 1];
				tempDstNodeIndex[i * 2 + 1] = temp;
                         tempEdgePhi[i] = - tempEdgePhi[i];
			}
		}
	}
	srcEdgeIndex = Eigen::VectorXi::Zero(NbrEdges);
	dstEdgeIndex = Eigen::VectorXi::Zero(NbrEdges);
	edgePhi = Eigen::VectorXcd::Zero(NbrEdges);
	for (int i = 0; i < NbrEdges; i++)
	{
		this->GetCoonOfEdges(tempSrcEdgeIndex[i], coonOfEdges);
		Tage_numEdge = GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1));
		srcEdgeIndex(i) = Tage_numEdge;
		this->GetCoonOfEdges(tempDstEdgeIndex[i], coonOfEdges);
		Tage_numEdge = GetEdgeOfTri(coonOfEdges(0), coonOfEdges(1));
		dstEdgeIndex(i) = Tage_numEdge;
		edgePhi(i) = tempEdgePhi[i];
	}
	//<2F><><EFBFBD><EFBFBD><EFBFBD>ص<EFBFBD>
	Eigen::VectorXi overlapIndex1, overlapIndex2;
	overlapIndex1 = Eigen::VectorXi::Zero(NbrEdges * 2);
	overlapIndex2 = Eigen::VectorXi::Zero(NbrEdges * 2);
	int Nbr_overlap;
	Nbr_overlap = 0;
	for (int i = 0; i < NbrEdges * 2; i++)
	{
		for (int j = i + 1; j < NbrEdges * 2; j++)
		{

			if (tempSrcNodeIndex[i] == tempSrcNodeIndex[j])
			{
				if (tempDstNodeIndex[i] == tempDstNodeIndex[j])
				{
					overlapIndex1[Nbr_overlap] = i;
					overlapIndex2[Nbr_overlap] = j;
					Nbr_overlap = Nbr_overlap + 1;
					break;
				}
			}

		}
	}
	overlapIndex2.conservativeResize(Nbr_overlap);
	int tempNum = 0;
	srcNodeIndex = Eigen::VectorXi::Zero(NbrEdges * 2 - Nbr_overlap);
	dstNodeIndex = Eigen::VectorXi::Zero(NbrEdges * 2 - Nbr_overlap);
	nodePhi = Eigen::VectorXcd::Zero(NbrEdges * 2 - Nbr_overlap);
	for (int i = 0; i < NbrEdges * 2; i++)
	{
		if ((overlapIndex2.array() == i).any())
		{
			tempNum++;
			continue;
		}
		srcNodeIndex(i - tempNum) = tempSrcNodeIndex[i];
		dstNodeIndex(i - tempNum) = tempDstNodeIndex[i];
		nodePhi(i - tempNum) = tempNodePhi[i];
	}
	srcNodeIndex.conservativeResize(NbrEdges * 2 - Nbr_overlap);
	dstNodeIndex.conservativeResize(NbrEdges * 2 - Nbr_overlap);
	int blochIndex, blochtemp = -1;
	for (int i = 0; i < NbrEdges * 2 - Nbr_overlap; i++)
	{
		for (int j = i + 1; j < NbrEdges * 2 - Nbr_overlap; j++)
		{
			if (srcNodeIndex(i) == srcNodeIndex(j))
			{
				blochtemp = srcNodeIndex[i];
			}
		}
	}
	for (int i = 0; i < NbrEdges * 2 - Nbr_overlap; i++)
	{
		for (int j = i + 1; j < NbrEdges * 2 - Nbr_overlap; j++)
		{
			if (dstNodeIndex(i) == dstNodeIndex(j))
			{
				srcNodeIndex(i) = blochtemp;
				nodePhi(i) = nodePhi(i) * nodePhi(j);
				if (j != srcNodeIndex.size() - 1)
				{
					srcNodeIndex.segment(j, NbrEdges * 2 - Nbr_overlap - j - 1) = srcNodeIndex.segment(j + 1, NbrEdges * 2 - Nbr_overlap - j - 1);
					dstNodeIndex.segment(j, NbrEdges * 2 - Nbr_overlap - j - 1) = dstNodeIndex.segment(j + 1, NbrEdges * 2 - Nbr_overlap - j - 1);
				}
			}
		}
	}
	if (blochtemp != -1)
	{
		srcNodeIndex.conservativeResize(NbrEdges * 2 - Nbr_overlap - 1);
		dstNodeIndex.conservativeResize(NbrEdges * 2 - Nbr_overlap - 1);
	}
	delete[] tempSrcEdgeIndex, tempSrcNodeIndex, tempEdgePhi, tempNodePhi, tempDstEdgeIndex, tempDstNodeIndex;
}

void Mesh_2D::GetEdgesIndexOfDomain(int domain, Eigen::VectorXi& edgesIndex)
{
	int NbrEdgesIndex = 0;
	for (int j = 0; j < _mNbrEdges; j++)
	{
		if (GetDomainOfEdges(j) == domain)
			NbrEdgesIndex++;
	}

	edgesIndex = Eigen::VectorXi::Zero(NbrEdgesIndex);
	int NumEdgesIndex = 0;
	for (int j = 0; j < _mNbrEdges; j++)
	{
		if (GetDomainOfEdges(j) == domain)
		{
			edgesIndex(NumEdgesIndex) = j;
			NumEdgesIndex++;
		}
	}
}

void Mesh_3D::GetTriIndexOfDomain(Eigen::VectorXi domain, Eigen::VectorXi& edgeIndex)
{
	int NbrTri = 0;

	int NbrDomain = domain.rows();
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrTri; j++)
		{
			if (GetDomainOfTri(j) == domain(i))
				NbrTri++;
		}
	}

	//store edgeIndex
	Eigen::Vector3i coonOfTri;
	int* tempEdgeIndex;
	tempEdgeIndex = new int[NbrTri * 3];
	int numTri = 0;
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrTri; j++)
		{
			if (GetDomainOfTri(j) == domain(i))
			{
				Eigen::Vector3i connOfTri;
				this->GetCoonOfTri(j, connOfTri);
				if (connOfTri(0) == 0)
				{
					if (connOfTri(2) == 0)
					{
						tempEdgeIndex[numTri * 3] = this->GetEdgeOfTet(connOfTri(1), 0);
						tempEdgeIndex[numTri * 3 + 1] = this->GetEdgeOfTet(connOfTri(1), 1);
						tempEdgeIndex[numTri * 3 + 2] = this->GetEdgeOfTet(connOfTri(1), 2);
					}
					else if (connOfTri(2) == 1)
					{
						tempEdgeIndex[numTri * 3] = this->GetEdgeOfTet(connOfTri(1), 0);
						tempEdgeIndex[numTri * 3 + 1] = this->GetEdgeOfTet(connOfTri(1), 2);
						tempEdgeIndex[numTri * 3 + 2] = this->GetEdgeOfTet(connOfTri(1), 4);
					}
					else if (connOfTri(2) == 2)
					{
						tempEdgeIndex[numTri * 3] = this->GetEdgeOfTet(connOfTri(1), 1);
						tempEdgeIndex[numTri * 3 + 1] = this->GetEdgeOfTet(connOfTri(1), 2);
						tempEdgeIndex[numTri * 3 + 2] = this->GetEdgeOfTet(connOfTri(1), 5);
					}
					else if (connOfTri(2) == 3)
					{
						tempEdgeIndex[numTri * 3] = this->GetEdgeOfTet(connOfTri(1), 3);
						tempEdgeIndex[numTri * 3 + 1] = this->GetEdgeOfTet(connOfTri(1), 4);
						tempEdgeIndex[numTri * 3 + 2] = this->GetEdgeOfTet(connOfTri(1), 5);
					}
				}
				else
				{

				}
			}
		}
	}

	//sort and unique nodeIndex
	std::sort(tempEdgeIndex, tempEdgeIndex + NbrTri * 3);
	int NbrEdge = (std::unique(tempEdgeIndex, tempEdgeIndex + NbrTri * 3) - tempEdgeIndex);
	edgeIndex = Eigen::VectorXi::Zero(NbrEdge);
	for (int i = 0; i < NbrEdge; i++)
		edgeIndex(i) = tempEdgeIndex[i];

	delete[] tempEdgeIndex;
}

void Mesh_3D::GetTriIndexOfDomain2(Eigen::VectorXi domain, Eigen::VectorXi indexNum,
	Eigen::VectorXi& triIndex, Eigen::VectorXi& triNum)
{
	int NbrTri = 0;

	int NbrDomain = domain.rows();
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < _mNbrTri; j++)
		{
			if (GetDomainOfTri(j) == domain(i))
				NbrTri++;
		}
	}

	//store triIndex
	triIndex = Eigen::VectorXi::Zero(NbrTri);
	triNum = Eigen::VectorXi::Zero(NbrTri);
	int numTri = 0;
	for (int i = 0; i < NbrDomain; i++)
	{
		for (int j = 0; j < NbrTri; j++)
		{
			if (GetDomainOfTri(j) == domain(i))
			{
				triIndex(numTri) = j;
				triNum(numTri) = i;
				numTri++;
			}
		}
	}
}

void Mesh_3D::GetIndexOfPBC(Eigen::MatrixXi PBC, Eigen::VectorXcd PBCData,
	Eigen::VectorXi& srcEdgeIndex, Eigen::VectorXi& dstEdgeIndex, Eigen::VectorXcd& edgePhi)
{
	//δд
}