diff --git a/Common/engine_interface_base.h b/Common/engine_interface_base.h
index 6ce466b..15be488 100644
--- a/Common/engine_interface_base.h
+++ b/Common/engine_interface_base.h
@@ -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;
diff --git a/Common/processfields.cpp b/Common/processfields.cpp
index f52fc14..b486c32 100644
--- a/Common/processfields.cpp
+++ b/Common/processfields.cpp
@@ -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;
diff --git a/Common/processfields.h b/Common/processfields.h
index 3d4085b..d954484 100644
--- a/Common/processfields.h
+++ b/Common/processfields.h
@@ -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;
diff --git a/FDTD/engine_interface_fdtd.cpp b/FDTD/engine_interface_fdtd.cpp
index b8fd9de..b5155fb 100644
--- a/FDTD/engine_interface_fdtd.cpp
+++ b/FDTD/engine_interface_fdtd.cpp
@@ -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;
diff --git a/FDTD/engine_interface_fdtd.h b/FDTD/engine_interface_fdtd.h
index 4e824ac..83281be 100644
--- a/FDTD/engine_interface_fdtd.h
+++ b/FDTD/engine_interface_fdtd.h
@@ -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
diff --git a/openems.cpp b/openems.cpp
index 057b573..a3c15d2 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -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;
 }