add field processing for electric and magnetic flux densities

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
pull/32/head
Thorsten Liebig 2017-05-28 12:01:04 +02:00
parent 92939becd0
commit 6133dea5b0
6 changed files with 130 additions and 15 deletions

View File

@ -56,6 +56,10 @@ public:
virtual double* GetJField(const unsigned int* pos, double* out) const =0;
//! Get the total current density field by rot(H) at \p pos. \sa SetInterpolationType
virtual double* GetRotHField(const unsigned int* pos, double* out) const =0;
//! Get the (interpolated) electric flux density field at \p pos. \sa SetInterpolationType
virtual double* GetDField(const unsigned int* pos, double* out) const =0;
//! Get the (interpolated) magnetic flux density field at \p pos. \sa SetInterpolationType
virtual double* GetBField(const unsigned int* pos, double* out) const =0;
//! Calculate the electric field integral along a given line
virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const =0;

View File

@ -71,6 +71,10 @@ string ProcessFields::GetFieldNameByType(DumpType type)
return "J-Field";
case ROTH_FIELD_DUMP:
return "RotH-Field";
case D_FIELD_DUMP:
return "D-Field";
case B_FIELD_DUMP:
return "B-Field";
case SAR_LOCAL_DUMP:
return "SAR-local";
case SAR_1G_DUMP:
@ -95,6 +99,30 @@ bool ProcessFields::NeedConductivity() const
return false;
}
bool ProcessFields::NeedPermittivity() const
{
switch (m_DumpType)
{
case D_FIELD_DUMP:
return true;
default:
return false;
}
return false;
}
bool ProcessFields::NeedPermeability() const
{
switch (m_DumpType)
{
case B_FIELD_DUMP:
return true;
default:
return false;
}
return false;
}
void ProcessFields::InitProcess()
{
if (Enabled==false) return;
@ -333,6 +361,44 @@ FDTD_FLOAT**** ProcessFields::CalcField()
}
}
return field;
case D_FIELD_DUMP:
for (unsigned int i=0; i<numLines[0]; ++i)
{
pos[0]=posLines[0][i];
for (unsigned int j=0; j<numLines[1]; ++j)
{
pos[1]=posLines[1][j];
for (unsigned int k=0; k<numLines[2]; ++k)
{
pos[2]=posLines[2][k];
m_Eng_Interface->GetDField(pos,out);
field[0][i][j][k] = out[0];
field[1][i][j][k] = out[1];
field[2][i][j][k] = out[2];
}
}
}
return field;
case B_FIELD_DUMP:
for (unsigned int i=0; i<numLines[0]; ++i)
{
pos[0]=posLines[0][i];
for (unsigned int j=0; j<numLines[1]; ++j)
{
pos[1]=posLines[1][j];
for (unsigned int k=0; k<numLines[2]; ++k)
{
pos[2]=posLines[2][k];
m_Eng_Interface->GetBField(pos,out);
field[0][i][j][k] = out[0];
field[1][i][j][k] = out[1];
field[2][i][j][k] = out[2];
}
}
}
return field;
default:
cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
return field;

View File

@ -40,7 +40,7 @@ public:
Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP),
(conduction) electric current density (kappa*E) (J_FIELD_DUMP) and total current density (rotH)
*/
enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29};
enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, D_FIELD_DUMP=4, B_FIELD_DUMP=5, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29};
virtual std::string GetProcessingName() const {return "common field processing";}
@ -75,6 +75,8 @@ public:
static std::string GetFieldNameByType(DumpType type);
virtual bool NeedConductivity() const;
virtual bool NeedPermittivity() const;
virtual bool NeedPermeability() const;
protected:
DumpType m_DumpType;

View File

@ -47,6 +47,11 @@ double* Engine_Interface_FDTD::GetJField(const unsigned int* pos, double* out) c
return GetRawInterpolatedField(pos, out, 1);
}
double* Engine_Interface_FDTD::GetDField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedField(pos, out, 3);
}
double* Engine_Interface_FDTD::GetRotHField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedField(pos, out, 2);
@ -116,6 +121,27 @@ double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos,
}
double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedDualField(pos, out, 0);
}
double* Engine_Interface_FDTD::GetBField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedDualField(pos, out, 1);
}
double Engine_Interface_FDTD::GetRawDualField(unsigned int n, const unsigned int* pos, int type) const
{
double value = m_Eng->GetCurr(n,pos[0],pos[1],pos[2]);
double delta = m_Op->GetEdgeLength(n,pos,true);
if ((type==0) && (delta))
return value/delta;
if ((type==1) && (m_Op->m_mueR) && (delta))
return value*m_Op->m_mueR[n][pos[0]][pos[1]][pos[2]]/delta;
return 0.0;
}
double* Engine_Interface_FDTD::GetRawInterpolatedDualField(const unsigned int* pos, double* out, int type) const
{
unsigned int iPos[] = {pos[0],pos[1],pos[2]};
int nP,nPP;
@ -124,9 +150,9 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
{
default:
case NO_INTERPOLATION:
out[0] = m_Eng->GetCurr(0,pos) / m_Op->GetEdgeLength(0,pos,true);
out[1] = m_Eng->GetCurr(1,pos) / m_Op->GetEdgeLength(1,pos,true);
out[2] = m_Eng->GetCurr(2,pos) / m_Op->GetEdgeLength(2,pos,true);
out[0] = GetRawDualField(0, pos, type);
out[1] = GetRawDualField(1, pos, type);
out[2] = GetRawDualField(2, pos, type);
break;
case NODE_INTERPOLATE:
for (int n=0; n<3; ++n)
@ -138,13 +164,13 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
out[n] = 0;
continue;
}
out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
out[n] = GetRawDualField(n, iPos, type);
--iPos[nP];
out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
out[n]+= GetRawDualField(n, iPos, type);
--iPos[nPP];
out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
out[n]+= GetRawDualField(n, iPos, type);
++iPos[nP];
out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
out[n]+= GetRawDualField(n, iPos, type);
++iPos[nPP];
out[n]/=4;
}
@ -153,7 +179,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
for (int n=0; n<3; ++n)
{
delta = m_Op->GetEdgeLength(n,iPos,true);
out[n] = m_Eng->GetCurr(n,iPos);
out[n] = GetRawDualField(n, iPos, type);
if ((pos[n]>=m_Op->GetNumberOfLines(n,true)-1))
{
out[n] = 0; //magnetic field on the outer boundaries is always zero
@ -162,7 +188,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
++iPos[n];
double deltaUp = m_Op->GetEdgeLength(n,iPos,true);
double deltaRel = delta / (delta+deltaUp);
out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetCurr(n,iPos)/deltaUp*deltaRel;
out[n] = out[n]*(1.0-deltaRel) + (double)GetRawDualField(n, iPos, type)*deltaRel;
--iPos[n];
}
break;
@ -201,6 +227,8 @@ double Engine_Interface_FDTD::GetRawField(unsigned int n, const unsigned int* po
return value/delta;
if ((type==1) && (m_Op->m_kappa) && (delta))
return value*m_Op->m_kappa[n][pos[0]][pos[1]][pos[2]]/delta;
if ((type==3) && (m_Op->m_epsR) && (delta))
return value*m_Op->m_epsR[n][pos[0]][pos[1]][pos[2]]/delta;
if (type==2) //calc rot(H)
{
int nP = (n+1)%3;

View File

@ -44,6 +44,8 @@ public:
virtual double* GetHField(const unsigned int* pos, double* out) const;
virtual double* GetJField(const unsigned int* pos, double* out) const;
virtual double* GetRotHField(const unsigned int* pos, double* out) const;
virtual double* GetDField(const unsigned int* pos, double* out) const;
virtual double* GetBField(const unsigned int* pos, double* out) const;
virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const;
@ -56,10 +58,15 @@ protected:
Operator* m_Op;
Engine* m_Eng;
//! Internal method to get an interpolated field of a given type. (0: E, 1: J, 2: rotH)
//! Internal method to get an interpolated field of a given type. (0: E, 1: J, 2: rotH, 3: D)
virtual double* GetRawInterpolatedField(const unsigned int* pos, double* out, int type) const;
//! Internal method to get a raw field of a given type. (0: E, 1: J, 2: rotH)
//! Internal method to get a raw field of a given type. (0: E, 1: J, 2: rotH, 3: D)
virtual double GetRawField(unsigned int n, const unsigned int* pos, int type) const;
//! Internal method to get an interpolated dual field of a given type. (0: H, 1: B)
virtual double* GetRawInterpolatedDualField(const unsigned int* pos, double* out, int type) const;
//! Internal method to get a raw dual field of a given type. (0: H, 1: B)
virtual double GetRawDualField(unsigned int n, const unsigned int* pos, int type) const;
};
#endif // ENGINE_INTERFACE_FDTD_H

View File

@ -472,9 +472,9 @@ bool openEMS::SetupProcessing()
CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox();
if (db)
{
if ((db->GetDumpType()>=0) && (db->GetDumpType()<=3))
if ((db->GetDumpType()>=0) && (db->GetDumpType()<=5))
ProcField = new ProcessFieldsTD(NewEngineInterface(db->GetMultiGridLevel()));
else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13))
else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=15))
ProcField = new ProcessFieldsFD(NewEngineInterface(db->GetMultiGridLevel()));
else if ( ((db->GetDumpType()>=20) && (db->GetDumpType()<=22)) || (db->GetDumpType()==29) )
{
@ -518,8 +518,12 @@ bool openEMS::SetupProcessing()
ProcField->SetDumpType(ProcessFields::SAR_RAW_DATA);
//SetupMaterialStorages() has previewed storage needs... refresh here to prevent cleanup!!!
if ( ProcField->NeedConductivity() && Enable_Dumps )
if ( ProcField->NeedPermittivity() && Enable_Dumps)
FDTD_Op->SetMaterialStoreFlags(0,true);
if ( ProcField->NeedConductivity() && Enable_Dumps)
FDTD_Op->SetMaterialStoreFlags(1,true);
if ( ProcField->NeedPermeability() && Enable_Dumps)
FDTD_Op->SetMaterialStoreFlags(2,true);
ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());
ProcField->SetFileType((ProcessFields::FileType)db->GetFileType());
@ -567,6 +571,10 @@ bool openEMS::SetupMaterialStorages()
(db->GetDumpType()==20) || (db->GetDumpType()==21) || (db->GetDumpType()==22)) && // SAR dump types
Enable_Dumps )
FDTD_Op->SetMaterialStoreFlags(1,true); //tell operator to store kappa material data
if ( ((db->GetDumpType()==4) || (db->GetDumpType()==14)) || Enable_Dumps) // electric flux density storage
FDTD_Op->SetMaterialStoreFlags(0,true); //tell operator to store epsR material data
if ( ((db->GetDumpType()==5) || (db->GetDumpType()==15)) || Enable_Dumps) // magnetic flux density storage
FDTD_Op->SetMaterialStoreFlags(2,true); //tell operator to store mueR material data
}
return true;
}