zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

new dump type: electric current density (J = kappa * E)

pull/1/head
Thorsten Liebig 2011-01-07 16:12:07 +01:00
parent ea496b6129
commit 2c3ebe5a7d
7 changed files with 83 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Common/engine_interface_base.h
View File

@ -50,6 +50,8 @@ public:
virtual double* GetEField(const unsigned int* pos, double* out) const =0; virtual double* GetEField(const unsigned int* pos, double* out) const =0;
//! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType //! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType
virtual double* GetHField(const unsigned int* pos, double* out) const =0; virtual double* GetHField(const unsigned int* pos, double* out) const =0;
//! Get the (interpolated) electric current density field at \p pos. \sa SetInterpolationType
virtual double* GetJField(const unsigned int* pos, double* out) const =0;
//! Calculate the electric field integral along a given line //! Calculate the electric field integral along a given line
virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const =0; virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const =0;

29
Common/processfields.cpp
View File

@ -63,6 +63,8 @@ string ProcessFields::GetFieldNameByType(DumpType type)
return "E-Field"; return "E-Field";
case H_FIELD_DUMP: case H_FIELD_DUMP:
return "H-Field"; return "H-Field";
case J_FIELD_DUMP:
return "J-Field";
} }
return "unknown field"; return "unknown field";
} }
@ -539,8 +541,9 @@ FDTD_FLOAT**** ProcessFields::CalcField()
double out[3]; double out[3];
//create array //create array
FDTD_FLOAT**** field = Create_N_3DArray<FDTD_FLOAT>(numLines); FDTD_FLOAT**** field = Create_N_3DArray<FDTD_FLOAT>(numLines);
if (m_DumpType==E_FIELD_DUMP) switch (m_DumpType)
{ {
case E_FIELD_DUMP:
for (unsigned int i=0; i<numLines[0]; ++i) for (unsigned int i=0; i<numLines[0]; ++i)
{ {
pos[0]=posLines[0][i]; pos[0]=posLines[0][i];
@ -559,10 +562,7 @@ FDTD_FLOAT**** ProcessFields::CalcField()
} }
} }
return field; return field;
} case H_FIELD_DUMP:
if (m_DumpType==H_FIELD_DUMP)
{
for (unsigned int i=0; i<numLines[0]; ++i) for (unsigned int i=0; i<numLines[0]; ++i)
{ {
pos[0]=posLines[0][i]; pos[0]=posLines[0][i];
@ -581,6 +581,25 @@ FDTD_FLOAT**** ProcessFields::CalcField()
} }
} }
return field; return field;
case J_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->GetJField(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;
} }
cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl; cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;

8
Common/processfields.h
View File

@ -29,8 +29,14 @@ public:
ProcessFields(Engine_Interface_Base* eng_if); ProcessFields(Engine_Interface_Base* eng_if);
virtual ~ProcessFields(); virtual ~ProcessFields();
//! File type definition.
enum FileType { VTK_FILETYPE, HDF5_FILETYPE}; enum FileType { VTK_FILETYPE, HDF5_FILETYPE};
enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP};
//! Dump type definitions.
/*!
Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP) and (conduction) electric current density (kappa*E) (J_FIELD_DUMP).
*/
enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP};
virtual void InitProcess(); virtual void InitProcess();

36
FDTD/engine_interface_fdtd.cpp
View File

@ -28,6 +28,16 @@ Engine_Interface_FDTD::~Engine_Interface_FDTD()
} }
double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) const double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedField(pos, out, 0);
}
double* Engine_Interface_FDTD::GetJField(const unsigned int* pos, double* out) const
{
return GetRawInterpolatedField(pos, out, 1);
}
double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos, double* out, int type) const
{ {
unsigned int iPos[] = {pos[0],pos[1],pos[2]}; unsigned int iPos[] = {pos[0],pos[1],pos[2]};
int nP,nPP; int nP,nPP;
@ -40,7 +50,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
{ {
delta = m_Op->GetEdgeLength(n,pos,false); delta = m_Op->GetEdgeLength(n,pos,false);
if (delta) if (delta)
out[n] = m_Eng->GetVolt(n,pos) / delta; out[n] = GetRawField(n,pos,type) / delta;
else else
out[n] = 0.0; out[n] = 0.0;
} }
@ -49,7 +59,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
delta = m_Op->GetEdgeLength(n,iPos); delta = m_Op->GetEdgeLength(n,iPos);
out[n] = m_Eng->GetVolt(n,iPos); out[n] = GetRawField(n,iPos);
if (delta==0) if (delta==0)
{ {
out[n]=0; out[n]=0;
@ -63,7 +73,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
--iPos[n]; --iPos[n];
double deltaDown = m_Op->GetEdgeLength(n,iPos); double deltaDown = m_Op->GetEdgeLength(n,iPos);
double deltaRel = delta / (delta+deltaDown); double deltaRel = delta / (delta+deltaDown);
out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetVolt(n,iPos)/deltaDown*deltaRel; out[n] = out[n]*(1.0-deltaRel)/delta + (double)GetRawField(n,iPos)/deltaDown*deltaRel;
++iPos[n]; ++iPos[n];
} }
break; break;
@ -79,19 +89,19 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
} }
delta = m_Op->GetEdgeLength(n,iPos); delta = m_Op->GetEdgeLength(n,iPos);
if (delta) if (delta)
out[n]=m_Eng->GetVolt(n,iPos)/delta; out[n]=GetRawField(n,iPos)/delta;
++iPos[nP]; ++iPos[nP];
delta = m_Op->GetEdgeLength(n,iPos); delta = m_Op->GetEdgeLength(n,iPos);
if (delta) if (delta)
out[n]+=m_Eng->GetVolt(n,iPos)/delta; out[n]+=GetRawField(n,iPos)/delta;
++iPos[nPP]; ++iPos[nPP];
delta = m_Op->GetEdgeLength(n,iPos); delta = m_Op->GetEdgeLength(n,iPos);
if (delta) if (delta)
out[n]+=m_Eng->GetVolt(n,iPos)/delta; out[n]+=GetRawField(n,iPos)/delta;
--iPos[nP]; --iPos[nP];
delta = m_Op->GetEdgeLength(n,iPos); delta = m_Op->GetEdgeLength(n,iPos);
if (delta) if (delta)
out[n]+=m_Eng->GetVolt(n,iPos)/delta; out[n]+=GetRawField(n,iPos)/delta;
--iPos[nPP]; --iPos[nPP];
out[n]/=4; out[n]/=4;
} }
@ -176,3 +186,15 @@ double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, con
} }
return result; return result;
} }
double Engine_Interface_FDTD::GetRawField(unsigned int n, const unsigned int* pos, int type) const
{
double value = m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
if (type==0)
return value;
if ((type==1) && (m_Op->m_kappa))
return value*m_Op->m_kappa[n][pos[0]][pos[1]][pos[2]];
return 0.0;
}

4
FDTD/engine_interface_fdtd.h
View File

@ -42,6 +42,7 @@ public:
virtual double* GetEField(const unsigned int* pos, double* out) const; virtual double* GetEField(const unsigned int* pos, double* out) const;
virtual double* GetHField(const unsigned int* pos, double* out) const; virtual double* GetHField(const unsigned int* pos, double* out) const;
virtual double* GetJField(const unsigned int* pos, double* out) const;
virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const; virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const;
@ -51,6 +52,9 @@ public:
protected: protected:
Operator* m_Op; Operator* m_Op;
Engine* m_Eng; Engine* m_Eng;
double* GetRawInterpolatedField(const unsigned int* pos, double* out, int type = 0) const;
double GetRawField(unsigned int n, const unsigned int* pos, int type = 0) const;
}; };
#endif // ENGINE_INTERFACE_FDTD_H #endif // ENGINE_INTERFACE_FDTD_H

1
FDTD/operator.h
View File

@ -32,6 +32,7 @@ class TiXmlElement;
class Operator : public Operator_Base class Operator : public Operator_Base
{ {
friend class Engine; friend class Engine;
friend class Engine_Interface_FDTD;
friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
friend class Operator_Ext_PML_SF_Plane; friend class Operator_Ext_PML_SF_Plane;
public: public:

20
openems.cpp
View File

@ -398,6 +398,9 @@ int openEMS::SetupFDTD(const char* file)
if (timestep) if (timestep)
FDTD_Op->SetTimestep(timestep); FDTD_Op->SetTimestep(timestep);
//save kappa material properties, maybe used for dump
FDTD_Op->SetMaterialStoreFlags(1,true);
Operator::DebugFlags debugFlags = Operator::None; Operator::DebugFlags debugFlags = Operator::None;
if (DebugMat) if (DebugMat)
debugFlags |= Operator::debugMaterial; debugFlags |= Operator::debugMaterial;
@ -407,6 +410,9 @@ int openEMS::SetupFDTD(const char* file)
debugFlags |= Operator::debugPEC; debugFlags |= Operator::debugPEC;
FDTD_Op->CalcECOperator( debugFlags ); FDTD_Op->CalcECOperator( debugFlags );
//reset flag for kappa material properties, if no dump-box resets it to true, it will be cleaned up...
FDTD_Op->SetMaterialStoreFlags(1,false);
unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep()); unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep());
if ((maxTime_TS>0) && (maxTime_TS<NrTS)) if ((maxTime_TS>0) && (maxTime_TS<NrTS))
NrTS = maxTime_TS; NrTS = maxTime_TS;
@ -529,9 +535,9 @@ int openEMS::SetupFDTD(const char* file)
CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox();
if (db) if (db)
{ {
if ((db->GetDumpType()>=0) && (db->GetDumpType()<2)) if ((db->GetDumpType()>=0) && (db->GetDumpType()<=2))
ProcField = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng)); ProcField = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
else if ((db->GetDumpType()>=10) && (db->GetDumpType()<12)) else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=12))
ProcField = new ProcessFieldsFD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng)); ProcField = new ProcessFieldsFD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
else else
cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl; cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl;
@ -545,13 +551,17 @@ int openEMS::SetupFDTD(const char* file)
//make dualMesh the default mesh for h-field dumps, maybe overwritten by interpolation type (node-interpolation) //make dualMesh the default mesh for h-field dumps, maybe overwritten by interpolation type (node-interpolation)
ProcField->SetDualMesh(true); ProcField->SetDualMesh(true);
} }
if ((db->GetDumpType()==10) || (db->GetDumpType()==11))
ProcField->AddFrequency(db->GetFDSamples());
if (db->GetDumpType()>=10) if (db->GetDumpType()>=10)
{
ProcField->AddFrequency(db->GetFDSamples());
ProcField->SetDumpType((ProcessFields::DumpType)(db->GetDumpType()-10)); ProcField->SetDumpType((ProcessFields::DumpType)(db->GetDumpType()-10));
}
else else
ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType()); ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12)) && Enable_Dumps )
FDTD_Op->SetMaterialStoreFlags(1,true); //keep kappa material data (prevent cleanup)
ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode()); ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());
ProcField->SetFileType((ProcessFields::FileType)db->GetFileType()); ProcField->SetFileType((ProcessFields::FileType)db->GetFileType());
if (CylinderCoords) if (CylinderCoords)
@ -573,6 +583,8 @@ int openEMS::SetupFDTD(const char* file)
} }
} }
FDTD_Op->CleanupMaterialStorage();
CSX.WarnUnusedPrimitves(cerr); CSX.WarnUnusedPrimitves(cerr);
// dump all boxes (voltage, current, fields, ...) // dump all boxes (voltage, current, fields, ...)