mur_abc_extension: now allows material distribution

pull/1/head
Thorsten Liebig 2010-04-28 21:45:05 +02:00
parent 03d8935cb8
commit 138624d2d8
5 changed files with 76 additions and 14 deletions

View File

@ -31,7 +31,8 @@ Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Ex
m_LineNr = m_Op_mur->m_LineNr;
m_LineNr_Shift = m_Op_mur->m_LineNr_Shift;
m_Mur_Coeff = m_Op_mur->m_Mur_Coeff;
m_Mur_Coeff_nyP = m_Op_mur->m_Mur_Coeff_nyP;
m_Mur_Coeff_nyPP = m_Op_mur->m_Mur_Coeff_nyPP;
m_volt_nyP = Create2DArray(m_numLines);
m_volt_nyPP = Create2DArray(m_numLines);
@ -66,7 +67,6 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates()
{
if (IsActive()==false) return;
if (m_Eng==NULL) return;
if (m_Mur_Coeff==0) return;
unsigned int pos[] = {0,0,0};
unsigned int pos_shift[] = {0,0,0};
pos[m_ny] = m_LineNr;
@ -78,8 +78,8 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates()
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
{
pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyP,pos);
m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos);
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos);
m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos);
}
}
}
@ -88,7 +88,6 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates()
{
if (IsActive()==false) return;
if (m_Eng==NULL) return;
if (m_Mur_Coeff==0) return;
unsigned int pos[] = {0,0,0};
unsigned int pos_shift[] = {0,0,0};
pos[m_ny] = m_LineNr;
@ -100,8 +99,8 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates()
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
{
pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Mur_Coeff * m_Eng->GetVolt(m_nyP,pos_shift);
m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos_shift);
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos_shift);
m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos_shift);
}
}
}
@ -110,7 +109,6 @@ void Engine_Ext_Mur_ABC::Apply2Voltages()
{
if (IsActive()==false) return;
if (m_Eng==NULL) return;
if (m_Mur_Coeff==0) return;
unsigned int pos[] = {0,0,0};
pos[m_ny] = m_LineNr;

View File

@ -46,7 +46,8 @@ protected:
int m_LineNr_Shift;
unsigned int m_numLines[2];
FDTD_FLOAT m_Mur_Coeff;
FDTD_FLOAT** m_Mur_Coeff_nyP;
FDTD_FLOAT** m_Mur_Coeff_nyPP;
FDTD_FLOAT** m_volt_nyP; //n+1 direction
FDTD_FLOAT** m_volt_nyPP; //n+2 direction
};

View File

@ -35,6 +35,7 @@ public:
virtual ~Operator();
virtual bool SetGeometryCSX(ContinuousStructure* geo);
virtual ContinuousStructure* GetGeometryCSX() {return CSX;}
virtual int CalcECOperator();

View File

@ -22,23 +22,35 @@
Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op)
{
m_Mur_Coeff = NULL;
m_ny = -1;
m_nyP = -1;
m_nyPP = -1;
m_LineNr = 0;
m_LineNr_Shift = 0;
m_Mur_Coeff = 0;
m_Mur_Coeff_nyP = NULL;
m_Mur_Coeff_nyPP = NULL;
m_numLines[0]=0;
m_numLines[1]=0;
}
Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC()
{
Delete2DArray(m_Mur_Coeff_nyP,m_numLines);
m_Mur_Coeff_nyP = NULL;
Delete2DArray(m_Mur_Coeff_nyPP,m_numLines);
m_Mur_Coeff_nyPP = NULL;
}
void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
{
if ((ny<0) || (ny>2))
return;
Delete2DArray(m_Mur_Coeff_nyP,m_numLines);
Delete2DArray(m_Mur_Coeff_nyPP,m_numLines);
m_ny = ny;
m_nyP = (ny+1)%3;
m_nyPP = (ny+2)%3;
@ -55,6 +67,10 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
m_numLines[0] = m_Op->GetNumberOfLines(m_nyP);
m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP);
m_Mur_Coeff_nyP = Create2DArray(m_numLines);
m_Mur_Coeff_nyPP = Create2DArray(m_numLines);
}
bool Operator_Ext_Mur_ABC::BuildExtension()
@ -68,8 +84,53 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
unsigned int pos[] = {0,0,0};
pos[m_ny] = m_LineNr;
double delta = fabs(m_Op->GetMeshDelta(m_ny,pos));
m_Mur_Coeff = (__C0__ * dT - delta) / (__C0__ * dT + delta);
// cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_ny << " @ " << m_LineNr << " = "<< m_Mur_Coeff << endl;
double coord[] = {0,0,0};
coord[0] = m_Op->GetDiscLine(0,pos[0]);
coord[1] = m_Op->GetDiscLine(1,pos[1]);
coord[2] = m_Op->GetDiscLine(2,pos[2]);
double eps,mue;
double c0t;
if (m_LineNr==0)
coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) + delta/2 / m_Op->GetGridDelta();
else
coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta();
for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP])
{
coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]);
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
{
coord[m_nyPP] = m_Op->GetDiscLine(m_nyPP,pos[m_nyPP]);
CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
if (prop)
{
CSPropMaterial* mat = prop->ToMaterial();
//nP
eps = mat->GetEpsilonWeighted(m_nyP,coord);
mue = mat->GetMueWeighted(m_nyP,coord);
c0t = __C0__ * dT / sqrt(eps*mue);
m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
//nPP
eps = mat->GetEpsilonWeighted(m_nyPP,coord);
mue = mat->GetMueWeighted(m_nyPP,coord);
c0t = __C0__ * dT / sqrt(eps*mue);
m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
}
else
{
c0t = __C0__ * dT;
m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]];
}
// cerr << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << " : " << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << endl;
}
}
// cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_ny << " @ " << m_LineNr << endl;
return true;
}

View File

@ -43,7 +43,8 @@ protected:
unsigned int m_numLines[2];
FDTD_FLOAT m_Mur_Coeff;
FDTD_FLOAT** m_Mur_Coeff_nyP;
FDTD_FLOAT** m_Mur_Coeff_nyPP;
};
#endif // OPERATOR_EXT_MUR_ABC_H