Major Change: Moved CartOperator into Operator
No more seperate Operator for cartesian and future cylindrical operator.
This commit is contained in:
Thorsten Liebig
parent
161bdf9c28
commit
5208d3f9a2
@ -1,529 +0,0 @@
|
|||||||
#include "cartoperator.h"
|
|
||||||
#include "tools/array_ops.h"
|
|
||||||
|
|
||||||
CartOperator::CartOperator() : Operator()
|
|
||||||
{
|
|
||||||
Init();
|
|
||||||
}
|
|
||||||
|
|
||||||
CartOperator::~CartOperator()
|
|
||||||
{
|
|
||||||
Reset();
|
|
||||||
}
|
|
||||||
|
|
||||||
void CartOperator::Init()
|
|
||||||
{
|
|
||||||
MainOp=NULL;
|
|
||||||
DualOp=NULL;
|
|
||||||
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
EC_C[n]=NULL;
|
|
||||||
EC_G[n]=NULL;
|
|
||||||
EC_L[n]=NULL;
|
|
||||||
EC_R[n]=NULL;
|
|
||||||
}
|
|
||||||
Operator::Init();
|
|
||||||
}
|
|
||||||
|
|
||||||
void CartOperator::Reset()
|
|
||||||
{
|
|
||||||
delete MainOp;
|
|
||||||
delete DualOp;
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
delete[] EC_C[n];
|
|
||||||
delete[] EC_G[n];
|
|
||||||
delete[] EC_L[n];
|
|
||||||
delete[] EC_R[n];
|
|
||||||
}
|
|
||||||
Operator::Reset();
|
|
||||||
Init();
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
void CartOperator::SetGeometryCSX(ContinuousStructure* geo)
|
|
||||||
{
|
|
||||||
if (geo==NULL) return;
|
|
||||||
|
|
||||||
Reset();
|
|
||||||
CSX = geo;
|
|
||||||
|
|
||||||
CSRectGrid* grid=CSX->GetGrid();
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
|
|
||||||
if (numLines[n]<3) {cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; Reset(); return;}
|
|
||||||
}
|
|
||||||
MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
|
|
||||||
MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
|
|
||||||
if (grid->GetDeltaUnit()<=0) {cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl; Reset(); return;}
|
|
||||||
else gridDelta=grid->GetDeltaUnit();
|
|
||||||
MainOp->SetGridDelta(1);
|
|
||||||
MainOp->AddCellAdrOp();
|
|
||||||
}
|
|
||||||
|
|
||||||
void CartOperator::InitOperator()
|
|
||||||
{
|
|
||||||
Delete_N_3DArray(vv,numLines);
|
|
||||||
Delete_N_3DArray(vi,numLines);
|
|
||||||
Delete_N_3DArray(iv,numLines);
|
|
||||||
Delete_N_3DArray(ii,numLines);
|
|
||||||
vv = Create_N_3DArray(numLines);
|
|
||||||
vi = Create_N_3DArray(numLines);
|
|
||||||
iv = Create_N_3DArray(numLines);
|
|
||||||
ii = Create_N_3DArray(numLines);
|
|
||||||
}
|
|
||||||
|
|
||||||
int CartOperator::CalcECOperator()
|
|
||||||
{
|
|
||||||
if (Calc_EC()==0)
|
|
||||||
return -1;
|
|
||||||
|
|
||||||
CalcTimestep();
|
|
||||||
|
|
||||||
InitOperator();
|
|
||||||
|
|
||||||
unsigned int i=0;
|
|
||||||
unsigned int pos[3];
|
|
||||||
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
|
||||||
{
|
|
||||||
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
|
||||||
{
|
|
||||||
i = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
|
|
||||||
vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
|
|
||||||
|
|
||||||
ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
|
|
||||||
iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
//cleanup
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
delete[] EC_C[n];EC_C[n]=NULL;
|
|
||||||
delete[] EC_G[n];EC_G[n]=NULL;
|
|
||||||
delete[] EC_L[n];EC_L[n]=NULL;
|
|
||||||
delete[] EC_R[n];EC_R[n]=NULL;
|
|
||||||
}
|
|
||||||
|
|
||||||
//Always apply PEC to all boundary's
|
|
||||||
bool PEC[6]={1,1,1,1,1,1};
|
|
||||||
ApplyElectricBC(PEC);
|
|
||||||
|
|
||||||
if (CalcEFieldExcitation()==false) return -1;
|
|
||||||
CalcPEC();
|
|
||||||
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
void CartOperator::ApplyElectricBC(bool* dirs)
|
|
||||||
{
|
|
||||||
if (dirs==NULL) return;
|
|
||||||
unsigned int pos[3];
|
|
||||||
unsigned int ipos;
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
|
|
||||||
{
|
|
||||||
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
|
|
||||||
{
|
|
||||||
pos[n]=0;
|
|
||||||
vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
pos[n]=numLines[n]-1;
|
|
||||||
vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void CartOperator::ApplyMagneticBC(bool* dirs)
|
|
||||||
{
|
|
||||||
if (dirs==NULL) return;
|
|
||||||
unsigned int pos[3];
|
|
||||||
unsigned int ipos;
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
|
|
||||||
{
|
|
||||||
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
|
|
||||||
{
|
|
||||||
pos[n]=0;
|
|
||||||
ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
|
||||||
|
|
||||||
pos[n]=numLines[n]-2;
|
|
||||||
ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
|
|
||||||
{
|
|
||||||
double coord[3];
|
|
||||||
double shiftCoord[3];
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
coord[0] = discLines[0][pos[0]];
|
|
||||||
coord[1] = discLines[1][pos[1]];
|
|
||||||
coord[2] = discLines[2][pos[2]];
|
|
||||||
unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
double delta=MainOp->GetIndexDelta(n,pos[n]);
|
|
||||||
double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
|
|
||||||
double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
|
|
||||||
double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
|
|
||||||
double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
|
|
||||||
double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
|
|
||||||
|
|
||||||
//******************************* epsilon,kappa averaging *****************************//
|
|
||||||
//shift up-right
|
|
||||||
shiftCoord[n] = coord[n]+delta*0.5;
|
|
||||||
shiftCoord[nP] = coord[nP]+deltaP*0.25;
|
|
||||||
shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
|
|
||||||
CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[0] = 1*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] = 0;
|
|
||||||
}
|
|
||||||
//shift up-left
|
|
||||||
shiftCoord[n] = coord[n]+delta*0.5;
|
|
||||||
shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
|
|
||||||
shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
|
|
||||||
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[0] += 1*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
//shift down-right
|
|
||||||
shiftCoord[n] = coord[n]+delta*0.5;
|
|
||||||
shiftCoord[nP] = coord[nP]+deltaP*0.25;
|
|
||||||
shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
|
|
||||||
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[0] += 1*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
//shift down-left
|
|
||||||
shiftCoord[n] = coord[n]+delta*0.5;
|
|
||||||
shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
|
|
||||||
shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
|
|
||||||
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[0] += 1*fabs(deltaP*deltaPP);
|
|
||||||
inEC[1] += 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
inEC[0]*=gridDelta/fabs(delta)/4*__EPS0__;
|
|
||||||
inEC[1]*=gridDelta/fabs(delta)/4;
|
|
||||||
|
|
||||||
//******************************* mu,sigma averaging *****************************//
|
|
||||||
//shift down
|
|
||||||
shiftCoord[n] = coord[n]-delta_M*0.25;
|
|
||||||
shiftCoord[nP] = coord[nP]+deltaP*0.5;
|
|
||||||
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
|
|
||||||
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[2] = fabs(delta_M) / mat->GetMueWeighted(n,shiftCoord);
|
|
||||||
if (mat->GetSigma(n))
|
|
||||||
inEC[3] = fabs(delta_M) / mat->GetSigmaWeighted(n,shiftCoord);
|
|
||||||
else
|
|
||||||
inEC[3] = 0;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[2] = fabs(delta_M);
|
|
||||||
inEC[3] = 0;
|
|
||||||
}
|
|
||||||
//shift up
|
|
||||||
shiftCoord[n] = coord[n]+delta*0.25;
|
|
||||||
shiftCoord[nP] = coord[nP]+deltaP*0.5;
|
|
||||||
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
|
|
||||||
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropMaterial* mat = prop->ToMaterial();
|
|
||||||
inEC[2] += mat->GetMue(n)*fabs(delta);
|
|
||||||
if (mat->GetSigmaWeighted(n,shiftCoord))
|
|
||||||
inEC[3] += fabs(delta)/mat->GetSigmaWeighted(n,shiftCoord);
|
|
||||||
else
|
|
||||||
inEC[3] = 0;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
inEC[2] += 1*fabs(delta);
|
|
||||||
inEC[3] = 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2 * __MUE0__ / inEC[2];
|
|
||||||
if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2 / inEC[3];
|
|
||||||
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool CartOperator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
|
|
||||||
{
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
|
|
||||||
unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
double delta=MainOp->GetIndexDelta(n,pos[n]);
|
|
||||||
double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
|
|
||||||
double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
|
|
||||||
|
|
||||||
double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
|
|
||||||
double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
|
|
||||||
double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
|
|
||||||
|
|
||||||
this->Calc_ECPos(n,pos,inMat);
|
|
||||||
|
|
||||||
inMat[0] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
|
|
||||||
inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
|
|
||||||
|
|
||||||
inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
|
|
||||||
inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
|
|
||||||
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
bool CartOperator::Calc_EC()
|
|
||||||
{
|
|
||||||
if (CSX==NULL) {cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false;}
|
|
||||||
|
|
||||||
unsigned int ipos;
|
|
||||||
unsigned int pos[3];
|
|
||||||
double inEC[4];
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
//init x-cell-array
|
|
||||||
delete[] EC_C[n];
|
|
||||||
delete[] EC_G[n];
|
|
||||||
delete[] EC_L[n];
|
|
||||||
delete[] EC_R[n];
|
|
||||||
EC_C[n] = new double[MainOp->GetSize()];
|
|
||||||
EC_G[n] = new double[MainOp->GetSize()];
|
|
||||||
EC_L[n] = new double[MainOp->GetSize()];
|
|
||||||
EC_R[n] = new double[MainOp->GetSize()];
|
|
||||||
for (unsigned int i=0;i<MainOp->GetSize();i++) //init all
|
|
||||||
{
|
|
||||||
EC_C[n][i]=0;
|
|
||||||
EC_G[n][i]=0;
|
|
||||||
EC_L[n][i]=0;
|
|
||||||
EC_R[n][i]=0;
|
|
||||||
}
|
|
||||||
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
|
||||||
{
|
|
||||||
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
|
||||||
{
|
|
||||||
Calc_ECPos(n,pos,inEC);
|
|
||||||
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
EC_C[n][ipos]=inEC[0];
|
|
||||||
EC_G[n][ipos]=inEC[1];
|
|
||||||
EC_L[n][ipos]=inEC[2];
|
|
||||||
EC_R[n][ipos]=inEC[3];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
double CartOperator::CalcTimestep()
|
|
||||||
{
|
|
||||||
dT=1e200;
|
|
||||||
double newT;
|
|
||||||
unsigned int pos[3];
|
|
||||||
unsigned int ipos;
|
|
||||||
unsigned int ipos_PM;
|
|
||||||
unsigned int ipos_PPM;
|
|
||||||
MainOp->SetReflection2Cell();
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
|
|
||||||
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
|
||||||
{
|
|
||||||
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
|
||||||
{
|
|
||||||
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
ipos_PM = MainOp->Shift(nP,-1);
|
|
||||||
MainOp->ResetShift();
|
|
||||||
ipos_PPM= MainOp->Shift(nPP,-1);
|
|
||||||
MainOp->ResetShift();
|
|
||||||
newT = 2/sqrt( ( 4/EC_L[nP][ipos] + 4/EC_L[nP][ipos_PPM] + 4/EC_L[nPP][ipos] + 4/EC_L[nPP][ipos_PM]) / EC_C[n][ipos] );
|
|
||||||
if (newT<dT) dT=newT;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// cerr << "Timestep: " << dT << endl;
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool CartOperator::CalcEFieldExcitation()
|
|
||||||
{
|
|
||||||
if (dT==0) return false;
|
|
||||||
vector<unsigned int> vIndex[3];
|
|
||||||
vector<FDTD_FLOAT> vExcit[3];
|
|
||||||
vector<unsigned int> vDelay;
|
|
||||||
unsigned int ipos;
|
|
||||||
unsigned int pos[3];
|
|
||||||
double coord[3];
|
|
||||||
|
|
||||||
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
|
||||||
{
|
|
||||||
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
|
||||||
{
|
|
||||||
coord[0] = discLines[0][pos[0]];
|
|
||||||
coord[1] = discLines[1][pos[1]];
|
|
||||||
coord[2] = discLines[2][pos[2]];
|
|
||||||
CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::ELECTRODE);
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
CSPropElectrode* elec = prop->ToElectrode();
|
|
||||||
if ((elec->GetExcitType()==0) || (elec->GetExcitType()==1)) //soft or hard E-Field excite!
|
|
||||||
{
|
|
||||||
vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
vIndex[n].push_back(pos[n]);
|
|
||||||
double delta=MainOp->GetIndexDelta(n,pos[n])*gridDelta;
|
|
||||||
if (elec->GetActiveDir(n))
|
|
||||||
vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta);
|
|
||||||
else
|
|
||||||
vExcit[n].push_back(0);
|
|
||||||
if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite
|
|
||||||
{
|
|
||||||
vv[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
vi[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
vv[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
vi[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
E_Ex_Count = vIndex[0].size();
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
delete[] E_Ex_index[n];
|
|
||||||
E_Ex_index[n] = new unsigned int[E_Ex_Count];
|
|
||||||
for (unsigned int i=0;i<E_Ex_Count;++i)
|
|
||||||
E_Ex_index[n][i]=vIndex[n].at(i);
|
|
||||||
}
|
|
||||||
delete[] E_Ex_delay;
|
|
||||||
E_Ex_delay = new unsigned int[E_Ex_Count];
|
|
||||||
for (unsigned int i=0;i<E_Ex_Count;++i)
|
|
||||||
E_Ex_delay[i]=vDelay.at(i);
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
delete[] E_Ex_amp[n];
|
|
||||||
E_Ex_amp[n] = new FDTD_FLOAT[E_Ex_Count];
|
|
||||||
for (unsigned int i=0;i<E_Ex_Count;++i)
|
|
||||||
E_Ex_amp[n][i]=vExcit[n].at(i);
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool CartOperator::CalcPEC()
|
|
||||||
{
|
|
||||||
unsigned int pos[3];
|
|
||||||
double coord[3];
|
|
||||||
double delta;
|
|
||||||
|
|
||||||
for (int n=0;n<3;++n)
|
|
||||||
{
|
|
||||||
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
|
||||||
{
|
|
||||||
coord[2] = discLines[2][pos[2]];
|
|
||||||
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
|
||||||
{
|
|
||||||
coord[1] = discLines[1][pos[1]];
|
|
||||||
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
|
||||||
{
|
|
||||||
coord[0] = discLines[0][pos[0]];
|
|
||||||
MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
||||||
delta=MainOp->GetIndexDelta(n,pos[n]);
|
|
||||||
coord[n]= discLines[n][pos[n]] + delta*0.5;
|
|
||||||
CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
|
|
||||||
if (prop)
|
|
||||||
{
|
|
||||||
if (prop->GetType()==CSProperties::METAL) //set to PEC
|
|
||||||
{
|
|
||||||
vv[n][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
vi[n][pos[0]][pos[1]][pos[2]] = 0;
|
|
||||||
// cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; " << pos[1] << " ; " << pos[2] << endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@ -1,42 +0,0 @@
|
|||||||
#ifndef CARTOPERATOR_H
|
|
||||||
#define CARTOPERATOR_H
|
|
||||||
|
|
||||||
#include "operator.h"
|
|
||||||
|
|
||||||
class CartOperator : public Operator
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
CartOperator();
|
|
||||||
virtual ~CartOperator();
|
|
||||||
|
|
||||||
virtual void SetGeometryCSX(ContinuousStructure* geo);
|
|
||||||
|
|
||||||
virtual int CalcECOperator();
|
|
||||||
|
|
||||||
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
|
|
||||||
virtual void ApplyMagneticBC(bool* dirs);
|
|
||||||
|
|
||||||
virtual void Reset();
|
|
||||||
|
|
||||||
protected:
|
|
||||||
virtual void Init();
|
|
||||||
virtual void InitOperator();
|
|
||||||
|
|
||||||
AdrOp* MainOp;
|
|
||||||
AdrOp* DualOp;
|
|
||||||
|
|
||||||
virtual bool CalcEFieldExcitation();
|
|
||||||
virtual bool CalcPEC();
|
|
||||||
virtual double CalcTimestep();
|
|
||||||
|
|
||||||
//EC elements, internal only!
|
|
||||||
bool Calc_EC();
|
|
||||||
bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
|
|
||||||
bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
|
|
||||||
double* EC_C[3];
|
|
||||||
double* EC_G[3];
|
|
||||||
double* EC_L[3];
|
|
||||||
double* EC_R[3];
|
|
||||||
};
|
|
||||||
|
|
||||||
#endif // CARTOPERATOR_H
|
|
||||||
@ -28,6 +28,17 @@ void Operator::Init()
|
|||||||
E_Ex_amp[n]=NULL;
|
E_Ex_amp[n]=NULL;
|
||||||
E_Ex_index[n]=NULL;
|
E_Ex_index[n]=NULL;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
MainOp=NULL;
|
||||||
|
DualOp=NULL;
|
||||||
|
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
EC_C[n]=NULL;
|
||||||
|
EC_G[n]=NULL;
|
||||||
|
EC_L[n]=NULL;
|
||||||
|
EC_R[n]=NULL;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void Operator::Reset()
|
void Operator::Reset()
|
||||||
@ -44,6 +55,16 @@ void Operator::Reset()
|
|||||||
delete[] E_Ex_amp[n];
|
delete[] E_Ex_amp[n];
|
||||||
delete[] E_Ex_index[n];
|
delete[] E_Ex_index[n];
|
||||||
}
|
}
|
||||||
|
delete MainOp;
|
||||||
|
delete DualOp;
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
delete[] EC_C[n];
|
||||||
|
delete[] EC_G[n];
|
||||||
|
delete[] EC_L[n];
|
||||||
|
delete[] EC_R[n];
|
||||||
|
}
|
||||||
|
|
||||||
Operator::Init();
|
Operator::Init();
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -84,14 +105,6 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
|
|||||||
return ok;
|
return ok;
|
||||||
}
|
}
|
||||||
|
|
||||||
void Operator::SetGeometryCSX(ContinuousStructure* geo)
|
|
||||||
{
|
|
||||||
if (geo==NULL) return;
|
|
||||||
|
|
||||||
Reset();
|
|
||||||
CSX = geo;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Operator::GetNumberCells()
|
double Operator::GetNumberCells()
|
||||||
{
|
{
|
||||||
if (numLines)
|
if (numLines)
|
||||||
@ -155,3 +168,490 @@ void Operator::DumpOperator2File(string filename)
|
|||||||
|
|
||||||
file.close();
|
file.close();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void Operator::SetGeometryCSX(ContinuousStructure* geo)
|
||||||
|
{
|
||||||
|
if (geo==NULL) return;
|
||||||
|
|
||||||
|
Reset();
|
||||||
|
CSX = geo;
|
||||||
|
|
||||||
|
CSRectGrid* grid=CSX->GetGrid();
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
|
||||||
|
if (numLines[n]<3) {cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; Reset(); return;}
|
||||||
|
}
|
||||||
|
MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
|
||||||
|
MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
|
||||||
|
if (grid->GetDeltaUnit()<=0) {cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl; Reset(); return;}
|
||||||
|
else gridDelta=grid->GetDeltaUnit();
|
||||||
|
MainOp->SetGridDelta(1);
|
||||||
|
MainOp->AddCellAdrOp();
|
||||||
|
}
|
||||||
|
|
||||||
|
void Operator::InitOperator()
|
||||||
|
{
|
||||||
|
Delete_N_3DArray(vv,numLines);
|
||||||
|
Delete_N_3DArray(vi,numLines);
|
||||||
|
Delete_N_3DArray(iv,numLines);
|
||||||
|
Delete_N_3DArray(ii,numLines);
|
||||||
|
vv = Create_N_3DArray(numLines);
|
||||||
|
vi = Create_N_3DArray(numLines);
|
||||||
|
iv = Create_N_3DArray(numLines);
|
||||||
|
ii = Create_N_3DArray(numLines);
|
||||||
|
}
|
||||||
|
|
||||||
|
int Operator::CalcECOperator()
|
||||||
|
{
|
||||||
|
if (Calc_EC()==0)
|
||||||
|
return -1;
|
||||||
|
|
||||||
|
CalcTimestep();
|
||||||
|
|
||||||
|
InitOperator();
|
||||||
|
|
||||||
|
unsigned int i=0;
|
||||||
|
unsigned int pos[3];
|
||||||
|
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
||||||
|
{
|
||||||
|
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
||||||
|
{
|
||||||
|
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
||||||
|
{
|
||||||
|
i = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
|
||||||
|
vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
|
||||||
|
|
||||||
|
ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
|
||||||
|
iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
//cleanup
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
delete[] EC_C[n];EC_C[n]=NULL;
|
||||||
|
delete[] EC_G[n];EC_G[n]=NULL;
|
||||||
|
delete[] EC_L[n];EC_L[n]=NULL;
|
||||||
|
delete[] EC_R[n];EC_R[n]=NULL;
|
||||||
|
}
|
||||||
|
|
||||||
|
//Always apply PEC to all boundary's
|
||||||
|
bool PEC[6]={1,1,1,1,1,1};
|
||||||
|
ApplyElectricBC(PEC);
|
||||||
|
|
||||||
|
if (CalcEFieldExcitation()==false) return -1;
|
||||||
|
CalcPEC();
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
void Operator::ApplyElectricBC(bool* dirs)
|
||||||
|
{
|
||||||
|
if (dirs==NULL) return;
|
||||||
|
unsigned int pos[3];
|
||||||
|
unsigned int ipos;
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
int nP = (n+1)%3;
|
||||||
|
int nPP = (n+2)%3;
|
||||||
|
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
|
||||||
|
{
|
||||||
|
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
|
||||||
|
{
|
||||||
|
pos[n]=0;
|
||||||
|
vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
pos[n]=numLines[n]-1;
|
||||||
|
vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void Operator::ApplyMagneticBC(bool* dirs)
|
||||||
|
{
|
||||||
|
if (dirs==NULL) return;
|
||||||
|
unsigned int pos[3];
|
||||||
|
unsigned int ipos;
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
int nP = (n+1)%3;
|
||||||
|
int nPP = (n+2)%3;
|
||||||
|
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
|
||||||
|
{
|
||||||
|
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
|
||||||
|
{
|
||||||
|
pos[n]=0;
|
||||||
|
ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
|
||||||
|
|
||||||
|
pos[n]=numLines[n]-2;
|
||||||
|
ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
|
||||||
|
{
|
||||||
|
double coord[3];
|
||||||
|
double shiftCoord[3];
|
||||||
|
int nP = (n+1)%3;
|
||||||
|
int nPP = (n+2)%3;
|
||||||
|
coord[0] = discLines[0][pos[0]];
|
||||||
|
coord[1] = discLines[1][pos[1]];
|
||||||
|
coord[2] = discLines[2][pos[2]];
|
||||||
|
unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
double delta=MainOp->GetIndexDelta(n,pos[n]);
|
||||||
|
double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
|
||||||
|
double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
|
||||||
|
double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
|
||||||
|
double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
|
||||||
|
double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
|
||||||
|
|
||||||
|
//******************************* epsilon,kappa averaging *****************************//
|
||||||
|
//shift up-right
|
||||||
|
shiftCoord[n] = coord[n]+delta*0.5;
|
||||||
|
shiftCoord[nP] = coord[nP]+deltaP*0.25;
|
||||||
|
shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
|
||||||
|
CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[0] = 1*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] = 0;
|
||||||
|
}
|
||||||
|
//shift up-left
|
||||||
|
shiftCoord[n] = coord[n]+delta*0.5;
|
||||||
|
shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
|
||||||
|
shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
|
||||||
|
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[0] += 1*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
//shift down-right
|
||||||
|
shiftCoord[n] = coord[n]+delta*0.5;
|
||||||
|
shiftCoord[nP] = coord[nP]+deltaP*0.25;
|
||||||
|
shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
|
||||||
|
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[0] += 1*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
//shift down-left
|
||||||
|
shiftCoord[n] = coord[n]+delta*0.5;
|
||||||
|
shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
|
||||||
|
shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
|
||||||
|
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[0] += 1*fabs(deltaP*deltaPP);
|
||||||
|
inEC[1] += 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
inEC[0]*=gridDelta/fabs(delta)/4*__EPS0__;
|
||||||
|
inEC[1]*=gridDelta/fabs(delta)/4;
|
||||||
|
|
||||||
|
//******************************* mu,sigma averaging *****************************//
|
||||||
|
//shift down
|
||||||
|
shiftCoord[n] = coord[n]-delta_M*0.25;
|
||||||
|
shiftCoord[nP] = coord[nP]+deltaP*0.5;
|
||||||
|
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
|
||||||
|
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[2] = fabs(delta_M) / mat->GetMueWeighted(n,shiftCoord);
|
||||||
|
if (mat->GetSigma(n))
|
||||||
|
inEC[3] = fabs(delta_M) / mat->GetSigmaWeighted(n,shiftCoord);
|
||||||
|
else
|
||||||
|
inEC[3] = 0;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[2] = fabs(delta_M);
|
||||||
|
inEC[3] = 0;
|
||||||
|
}
|
||||||
|
//shift up
|
||||||
|
shiftCoord[n] = coord[n]+delta*0.25;
|
||||||
|
shiftCoord[nP] = coord[nP]+deltaP*0.5;
|
||||||
|
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
|
||||||
|
prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropMaterial* mat = prop->ToMaterial();
|
||||||
|
inEC[2] += mat->GetMue(n)*fabs(delta);
|
||||||
|
if (mat->GetSigmaWeighted(n,shiftCoord))
|
||||||
|
inEC[3] += fabs(delta)/mat->GetSigmaWeighted(n,shiftCoord);
|
||||||
|
else
|
||||||
|
inEC[3] = 0;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
inEC[2] += 1*fabs(delta);
|
||||||
|
inEC[3] = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2 * __MUE0__ / inEC[2];
|
||||||
|
if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2 / inEC[3];
|
||||||
|
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool Operator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
|
||||||
|
{
|
||||||
|
int nP = (n+1)%3;
|
||||||
|
int nPP = (n+2)%3;
|
||||||
|
|
||||||
|
unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
double delta=MainOp->GetIndexDelta(n,pos[n]);
|
||||||
|
double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
|
||||||
|
double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
|
||||||
|
|
||||||
|
double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
|
||||||
|
double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
|
||||||
|
double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
|
||||||
|
|
||||||
|
this->Calc_ECPos(n,pos,inMat);
|
||||||
|
|
||||||
|
inMat[0] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
|
||||||
|
inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
|
||||||
|
|
||||||
|
inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
|
||||||
|
inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
|
||||||
|
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
bool Operator::Calc_EC()
|
||||||
|
{
|
||||||
|
if (CSX==NULL) {cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false;}
|
||||||
|
|
||||||
|
unsigned int ipos;
|
||||||
|
unsigned int pos[3];
|
||||||
|
double inEC[4];
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
//init x-cell-array
|
||||||
|
delete[] EC_C[n];
|
||||||
|
delete[] EC_G[n];
|
||||||
|
delete[] EC_L[n];
|
||||||
|
delete[] EC_R[n];
|
||||||
|
EC_C[n] = new double[MainOp->GetSize()];
|
||||||
|
EC_G[n] = new double[MainOp->GetSize()];
|
||||||
|
EC_L[n] = new double[MainOp->GetSize()];
|
||||||
|
EC_R[n] = new double[MainOp->GetSize()];
|
||||||
|
for (unsigned int i=0;i<MainOp->GetSize();i++) //init all
|
||||||
|
{
|
||||||
|
EC_C[n][i]=0;
|
||||||
|
EC_G[n][i]=0;
|
||||||
|
EC_L[n][i]=0;
|
||||||
|
EC_R[n][i]=0;
|
||||||
|
}
|
||||||
|
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
||||||
|
{
|
||||||
|
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
||||||
|
{
|
||||||
|
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
||||||
|
{
|
||||||
|
Calc_ECPos(n,pos,inEC);
|
||||||
|
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
EC_C[n][ipos]=inEC[0];
|
||||||
|
EC_G[n][ipos]=inEC[1];
|
||||||
|
EC_L[n][ipos]=inEC[2];
|
||||||
|
EC_R[n][ipos]=inEC[3];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
double Operator::CalcTimestep()
|
||||||
|
{
|
||||||
|
dT=1e200;
|
||||||
|
double newT;
|
||||||
|
unsigned int pos[3];
|
||||||
|
unsigned int ipos;
|
||||||
|
unsigned int ipos_PM;
|
||||||
|
unsigned int ipos_PPM;
|
||||||
|
MainOp->SetReflection2Cell();
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
int nP = (n+1)%3;
|
||||||
|
int nPP = (n+2)%3;
|
||||||
|
|
||||||
|
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
||||||
|
{
|
||||||
|
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
||||||
|
{
|
||||||
|
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
||||||
|
{
|
||||||
|
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
ipos_PM = MainOp->Shift(nP,-1);
|
||||||
|
MainOp->ResetShift();
|
||||||
|
ipos_PPM= MainOp->Shift(nPP,-1);
|
||||||
|
MainOp->ResetShift();
|
||||||
|
newT = 2/sqrt( ( 4/EC_L[nP][ipos] + 4/EC_L[nP][ipos_PPM] + 4/EC_L[nPP][ipos] + 4/EC_L[nPP][ipos_PM]) / EC_C[n][ipos] );
|
||||||
|
if (newT<dT) dT=newT;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// cerr << "Timestep: " << dT << endl;
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool Operator::CalcEFieldExcitation()
|
||||||
|
{
|
||||||
|
if (dT==0) return false;
|
||||||
|
vector<unsigned int> vIndex[3];
|
||||||
|
vector<FDTD_FLOAT> vExcit[3];
|
||||||
|
vector<unsigned int> vDelay;
|
||||||
|
unsigned int ipos;
|
||||||
|
unsigned int pos[3];
|
||||||
|
double coord[3];
|
||||||
|
|
||||||
|
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
||||||
|
{
|
||||||
|
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
||||||
|
{
|
||||||
|
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
||||||
|
{
|
||||||
|
coord[0] = discLines[0][pos[0]];
|
||||||
|
coord[1] = discLines[1][pos[1]];
|
||||||
|
coord[2] = discLines[2][pos[2]];
|
||||||
|
CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::ELECTRODE);
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
CSPropElectrode* elec = prop->ToElectrode();
|
||||||
|
if ((elec->GetExcitType()==0) || (elec->GetExcitType()==1)) //soft or hard E-Field excite!
|
||||||
|
{
|
||||||
|
vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
vIndex[n].push_back(pos[n]);
|
||||||
|
double delta=MainOp->GetIndexDelta(n,pos[n])*gridDelta;
|
||||||
|
if (elec->GetActiveDir(n))
|
||||||
|
vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta);
|
||||||
|
else
|
||||||
|
vExcit[n].push_back(0);
|
||||||
|
if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite
|
||||||
|
{
|
||||||
|
vv[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
vi[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
vv[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
vi[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
E_Ex_Count = vIndex[0].size();
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
delete[] E_Ex_index[n];
|
||||||
|
E_Ex_index[n] = new unsigned int[E_Ex_Count];
|
||||||
|
for (unsigned int i=0;i<E_Ex_Count;++i)
|
||||||
|
E_Ex_index[n][i]=vIndex[n].at(i);
|
||||||
|
}
|
||||||
|
delete[] E_Ex_delay;
|
||||||
|
E_Ex_delay = new unsigned int[E_Ex_Count];
|
||||||
|
for (unsigned int i=0;i<E_Ex_Count;++i)
|
||||||
|
E_Ex_delay[i]=vDelay.at(i);
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
delete[] E_Ex_amp[n];
|
||||||
|
E_Ex_amp[n] = new FDTD_FLOAT[E_Ex_Count];
|
||||||
|
for (unsigned int i=0;i<E_Ex_Count;++i)
|
||||||
|
E_Ex_amp[n][i]=vExcit[n].at(i);
|
||||||
|
}
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool Operator::CalcPEC()
|
||||||
|
{
|
||||||
|
unsigned int pos[3];
|
||||||
|
double coord[3];
|
||||||
|
double delta;
|
||||||
|
|
||||||
|
for (int n=0;n<3;++n)
|
||||||
|
{
|
||||||
|
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
|
||||||
|
{
|
||||||
|
coord[2] = discLines[2][pos[2]];
|
||||||
|
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
|
||||||
|
{
|
||||||
|
coord[1] = discLines[1][pos[1]];
|
||||||
|
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
|
||||||
|
{
|
||||||
|
coord[0] = discLines[0][pos[0]];
|
||||||
|
MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||||
|
delta=MainOp->GetIndexDelta(n,pos[n]);
|
||||||
|
coord[n]= discLines[n][pos[n]] + delta*0.5;
|
||||||
|
CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
|
||||||
|
if (prop)
|
||||||
|
{
|
||||||
|
if (prop->GetType()==CSProperties::METAL) //set to PEC
|
||||||
|
{
|
||||||
|
vv[n][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
vi[n][pos[0]][pos[1]][pos[2]] = 0;
|
||||||
|
// cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; " << pos[1] << " ; " << pos[2] << endl;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|||||||
@ -19,12 +19,12 @@ public:
|
|||||||
|
|
||||||
virtual void SetGeometryCSX(ContinuousStructure* geo);
|
virtual void SetGeometryCSX(ContinuousStructure* geo);
|
||||||
|
|
||||||
virtual int CalcECOperator() {return -1;};
|
virtual int CalcECOperator();
|
||||||
|
|
||||||
virtual void CalcGaussianPulsExcitation(double f0, double fc);
|
virtual void CalcGaussianPulsExcitation(double f0, double fc);
|
||||||
|
|
||||||
virtual void ApplyElectricBC(bool* dirs) {}; //applied by default to all boundaries
|
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
|
||||||
virtual void ApplyMagneticBC(bool* dirs) {};
|
virtual void ApplyMagneticBC(bool* dirs);
|
||||||
|
|
||||||
double GetTimestep() {return dT;};
|
double GetTimestep() {return dT;};
|
||||||
unsigned int GetNyquistNum(double fmax);
|
unsigned int GetNyquistNum(double fmax);
|
||||||
@ -40,12 +40,17 @@ public:
|
|||||||
|
|
||||||
protected:
|
protected:
|
||||||
virtual void Init();
|
virtual void Init();
|
||||||
|
virtual void InitOperator();
|
||||||
|
|
||||||
ContinuousStructure* CSX;
|
ContinuousStructure* CSX;
|
||||||
double gridDelta;
|
double gridDelta;
|
||||||
|
|
||||||
double* discLines[3];
|
double* discLines[3];
|
||||||
unsigned int numLines[3];
|
unsigned int numLines[3];
|
||||||
|
|
||||||
|
AdrOp* MainOp;
|
||||||
|
AdrOp* DualOp;
|
||||||
|
|
||||||
//EC operator
|
//EC operator
|
||||||
FDTD_FLOAT**** vv; //calc new voltage from old voltage
|
FDTD_FLOAT**** vv; //calc new voltage from old voltage
|
||||||
FDTD_FLOAT**** vi; //calc new voltage from old current
|
FDTD_FLOAT**** vi; //calc new voltage from old current
|
||||||
@ -58,15 +63,26 @@ protected:
|
|||||||
|
|
||||||
//E-Field Excitation
|
//E-Field Excitation
|
||||||
//! Calc the electric field excitation.
|
//! Calc the electric field excitation.
|
||||||
virtual bool CalcEFieldExcitation() {return false;};
|
virtual bool CalcEFieldExcitation();
|
||||||
unsigned int E_Ex_Count;
|
unsigned int E_Ex_Count;
|
||||||
unsigned int* E_Ex_index[3];
|
unsigned int* E_Ex_index[3];
|
||||||
FDTD_FLOAT* E_Ex_amp[3]; //represented as edge-voltages!!
|
FDTD_FLOAT* E_Ex_amp[3]; //represented as edge-voltages!!
|
||||||
unsigned int* E_Ex_delay;
|
unsigned int* E_Ex_delay;
|
||||||
|
|
||||||
|
virtual bool CalcPEC();
|
||||||
|
|
||||||
//Calc timestep only internal use
|
//Calc timestep only internal use
|
||||||
virtual double CalcTimestep() {return 0;};
|
virtual double CalcTimestep();
|
||||||
double dT; //FDTD timestep!
|
double dT; //FDTD timestep!
|
||||||
|
|
||||||
|
//EC elements, internal only!
|
||||||
|
bool Calc_EC();
|
||||||
|
bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
|
||||||
|
bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
|
||||||
|
double* EC_C[3];
|
||||||
|
double* EC_G[3];
|
||||||
|
double* EC_L[3];
|
||||||
|
double* EC_R[3];
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif // OPERATOR_H
|
#endif // OPERATOR_H
|
||||||
|
|||||||
4
main.cpp
4
main.cpp
@ -3,7 +3,7 @@
|
|||||||
#include <sstream>
|
#include <sstream>
|
||||||
#include <time.h>
|
#include <time.h>
|
||||||
#include "tools/array_ops.h"
|
#include "tools/array_ops.h"
|
||||||
#include "FDTD/cartoperator.h"
|
#include "FDTD/operator.h"
|
||||||
#include "FDTD/engine.h"
|
#include "FDTD/engine.h"
|
||||||
#include "FDTD/processvoltage.h"
|
#include "FDTD/processvoltage.h"
|
||||||
#include "FDTD/processcurrent.h"
|
#include "FDTD/processcurrent.h"
|
||||||
@ -30,7 +30,7 @@ int main(int argc, char *argv[])
|
|||||||
|
|
||||||
//*************** setup operator ************//
|
//*************** setup operator ************//
|
||||||
cerr << "Create Operator..." << endl;
|
cerr << "Create Operator..." << endl;
|
||||||
CartOperator cop;
|
Operator cop;
|
||||||
cop.SetGeometryCSX(&CSX);
|
cop.SetGeometryCSX(&CSX);
|
||||||
|
|
||||||
cop.CalcECOperator();
|
cop.CalcECOperator();
|
||||||
|
|||||||
@ -8,7 +8,7 @@ CONFIG -= app_bundle
|
|||||||
TEMPLATE = app
|
TEMPLATE = app
|
||||||
OBJECTS_DIR = obj
|
OBJECTS_DIR = obj
|
||||||
INCLUDEPATH += ../CSXCAD \
|
INCLUDEPATH += ../CSXCAD \
|
||||||
../fparser
|
../fparser
|
||||||
LIBS += -L../CSXCAD \
|
LIBS += -L../CSXCAD \
|
||||||
-lCSXCAD \
|
-lCSXCAD \
|
||||||
-L../fparser \
|
-L../fparser \
|
||||||
@ -16,7 +16,6 @@ LIBS += -L../CSXCAD \
|
|||||||
-L../tinyxml \
|
-L../tinyxml \
|
||||||
-ltinyxml
|
-ltinyxml
|
||||||
SOURCES += main.cpp \
|
SOURCES += main.cpp \
|
||||||
FDTD/cartoperator.cpp \
|
|
||||||
tools/ErrorMsg.cpp \
|
tools/ErrorMsg.cpp \
|
||||||
tools/AdrOp.cpp \
|
tools/AdrOp.cpp \
|
||||||
FDTD/engine.cpp \
|
FDTD/engine.cpp \
|
||||||
@ -27,8 +26,7 @@ SOURCES += main.cpp \
|
|||||||
FDTD/processfields.cpp \
|
FDTD/processfields.cpp \
|
||||||
FDTD/processfields_td.cpp \
|
FDTD/processfields_td.cpp \
|
||||||
FDTD/processcurrent.cpp
|
FDTD/processcurrent.cpp
|
||||||
HEADERS += FDTD/cartoperator.h \
|
HEADERS += tools/ErrorMsg.h \
|
||||||
tools/ErrorMsg.h \
|
|
||||||
tools/AdrOp.h \
|
tools/AdrOp.h \
|
||||||
tools/constants.h \
|
tools/constants.h \
|
||||||
FDTD/engine.h \
|
FDTD/engine.h \
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user