lumpedRLC: replace all "uint" to "unsigned int".

The original contributor used "uint" in all code, which is a
non-standard language extension that doesn't exist in ISO C
or ISO C++, and causes build failures on macOS as reported by
PR #144. Replace all "uint" with "unsigned int" for standard
conformance to maximize portability.

Signed-off-by: Yifeng Li <tomli@tomli.me>

CMakeLists.txt

@@ -7,12 +7,12 @@ ELSE()
 ENDIF()
 
 PROJECT(openEMS CXX C)
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.0)
 
 # default
 set(LIB_VERSION_MAJOR 0)
 set(LIB_VERSION_MINOR 0)
-set(LIB_VERSION_PATCH 35)
+set(LIB_VERSION_PATCH 36)
 set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
 
 set(VERSION "v${LIB_VERSION_STRING}")
@@ -23,7 +23,7 @@ ENDIF()
 
 #ADD_DEFINITIONS( -D__SSE2__ )
 
-set(VERSION "v0.0.35")
+set(VERSION "v0.0.36")
 
 # add git revision
 IF(EXISTS ${PROJECT_SOURCE_DIR}/.git )
@@ -61,7 +61,6 @@ if(UNIX AND ENABLE_RPATH)
   set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
 endif()
 
-
 if (WITH_MPI)
 	ADD_DEFINITIONS(-DMPI_SUPPORT)
 	# Require MPI for this project:
@@ -138,14 +137,7 @@ INCLUDE_DIRECTORIES (${Boost_INCLUDE_DIRS})
 find_package(VTK COMPONENTS vtkIOXML vtkIOGeometry vtkIOLegacy vtkIOPLY NO_MODULE REQUIRED)
 
 message(STATUS "Found package VTK. Using version " ${VTK_VERSION})
-if("${VTK_MAJOR_VERSION}" GREATER 5)
+set( vtk_LIBS ${VTK_LIBRARIES} )
-    set( vtk_LIBS ${VTK_LIBRARIES} )
-else()
-    set( vtk_LIBS
-        vtkCommon
-        vtkIO
-    )    
-endif()
 message(STATUS "vtk libraries " ${vtk_LIBS})
 
 include(${VTK_USE_FILE})
@@ -159,7 +151,7 @@ elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "AMD64")
   set(ARCH "x86_64")
 elseif(CMAKE_SYSTEM_PROCESSOR MATCHES "^(ppc64.*|PPC64.*|powerpc64.*)")
   set(ARCH "ppc64")
-elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "aarch64")
+elseif(${CMAKE_SYSTEM_PROCESSOR} MATCHES "^(aarch64|arm64)")
   set(ARCH "aarch64")
 elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "unknown")
   set(ARCH "unknown")
@@ -170,13 +162,13 @@ endif()
 
 if(${ARCH} STREQUAL "x86_64")
   message(STATUS "Detected 64-bit x86 target")
-  #set(CMAKE_CXX_FLAGS "-msse -march=native")
+  #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse -march=native")
 elseif(${ARCH} STREQUAL "ppc64")
   message(STATUS "Detected 64-bit POWER target")
-  set(CMAKE_CXX_FLAGS "-DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
 elseif(${ARCH} STREQUAL "aarch64")
   message(STATUS "Detected 64-bit ARM target")
-  set(CMAKE_CXX_FLAGS "-DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
 elseif(${ARCH} STREQUAL "unsupported")
   message(FATAL_ERROR "Unsupported target architecture!  Try porting openEMS to your architecture...")
 else()

FDTD/engine.cpp

@@ -123,16 +123,37 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
 				shift[2]=pos[2];
 				//do the updates here
 				//for x
-				volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]];
+				volt[0][pos[0]][pos[1]][pos[2]] *=
-				volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]);
+				    Op->vv[0][pos[0]][pos[1]][pos[2]];
+				volt[0][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[0][pos[0]][pos[1]][pos[2]] * (
+				        curr[2][pos[0]][pos[1]         ][pos[2]         ] -
+				        curr[2][pos[0]][pos[1]-shift[1]][pos[2]         ] -
+				        curr[1][pos[0]][pos[1]         ][pos[2]         ] +
+				        curr[1][pos[0]][pos[1]         ][pos[2]-shift[2]]
+				    );
 
 				//for y
-				volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]];
+				volt[1][pos[0]][pos[1]][pos[2]] *=
-				volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]);
+				    Op->vv[1][pos[0]][pos[1]][pos[2]];
+				volt[1][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[1][pos[0]][pos[1]][pos[2]] * (
+				        curr[0][pos[0]         ][pos[1]][pos[2]         ] -
+				        curr[0][pos[0]         ][pos[1]][pos[2]-shift[2]] -
+				        curr[2][pos[0]         ][pos[1]][pos[2]         ] +
+				        curr[2][pos[0]-shift[0]][pos[1]][pos[2]         ]
+				    );
 
 				//for z
-				volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]];
+				volt[2][pos[0]][pos[1]][pos[2]] *=
-				volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]);
+				    Op->vv[2][pos[0]][pos[1]][pos[2]];
+				volt[2][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[2][pos[0]][pos[1]][pos[2]] * (
+				        curr[1][pos[0]         ][pos[1]         ][pos[2]] -
+				        curr[1][pos[0]-shift[0]][pos[1]         ][pos[2]] -
+				        curr[0][pos[0]         ][pos[1]         ][pos[2]] +
+				        curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]]
+				    );
 			}
 		}
 		++pos[0];
@@ -151,16 +172,37 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
 			{
 				//do the updates here
 				//for x
-				curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]];
+				curr[0][pos[0]][pos[1]][pos[2]] *=
-				curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]);
+				    Op->ii[0][pos[0]][pos[1]][pos[2]];
+				curr[0][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[0][pos[0]][pos[1]][pos[2]] * (
+				        volt[2][pos[0]][pos[1]  ][pos[2]  ] -
+				        volt[2][pos[0]][pos[1]+1][pos[2]  ] -
+				        volt[1][pos[0]][pos[1]  ][pos[2]  ] +
+				        volt[1][pos[0]][pos[1]  ][pos[2]+1]
+				    );
 
 				//for y
-				curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]];
+				curr[1][pos[0]][pos[1]][pos[2]] *=
-				curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]);
+				    Op->ii[1][pos[0]][pos[1]][pos[2]];
+				curr[1][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[1][pos[0]][pos[1]][pos[2]] * (
+				        volt[0][pos[0]  ][pos[1]][pos[2]  ] -
+				        volt[0][pos[0]  ][pos[1]][pos[2]+1] -
+				        volt[2][pos[0]  ][pos[1]][pos[2]  ] +
+				        volt[2][pos[0]+1][pos[1]][pos[2]  ]
+				    );
 
 				//for z
-				curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]];
+				curr[2][pos[0]][pos[1]][pos[2]] *=
-				curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]);
+				    Op->ii[2][pos[0]][pos[1]][pos[2]];
+				curr[2][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[2][pos[0]][pos[1]][pos[2]] * (
+				        volt[1][pos[0]  ][pos[1]  ][pos[2]] -
+				        volt[1][pos[0]+1][pos[1]  ][pos[2]] -
+				        volt[0][pos[0]  ][pos[1]  ][pos[2]] +
+				        volt[0][pos[0]  ][pos[1]+1][pos[2]]
+				    );
 			}
 		}
 		++pos[0];

FDTD/engine.h

@@ -47,6 +47,8 @@ public:
 
 	virtual unsigned int GetNumberOfTimesteps() {return numTS;}
 
+	virtual void NextInterval(float curr_speed) {};
+
 	//this access functions muss be overloaded by any new engine using a different storage model
 	inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z )		const { return volt[n][x][y][z]; }
 	inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] )							const { return volt[n][pos[0]][pos[1]][pos[2]]; }

FDTD/engine_multithread.cpp

@@ -55,6 +55,12 @@ Engine_Multithread::Engine_Multithread(const Operator_Multithread* op) : ENGINE_
 	m_IterateBarrier = 0;
 	m_startBarrier = 0;
 	m_stopBarrier = 0;
+	m_thread_group = 0;
+	m_max_numThreads = boost::thread::hardware_concurrency();
+	m_numThreads = 0;
+	m_last_speed = 0;
+	m_opt_speed = false;
+	m_stopThreads = true;
 
 #ifdef ENABLE_DEBUG_TIME
 	m_MPI_Barrier = 0;
@@ -92,27 +98,85 @@ void Engine_Multithread::setNumThreads( unsigned int numThreads )
 void Engine_Multithread::Init()
 {
 	m_stopThreads = true;
+	m_opt_speed = false;
 	ENGINE_MULTITHREAD_BASE::Init();
 
 	// initialize threads
 	m_stopThreads = false;
 	if (m_numThreads == 0)
-		m_numThreads = boost::thread::hardware_concurrency();
+	{
+		m_opt_speed = true;
+		m_numThreads = 1;
+	}
+	else if (m_numThreads > m_max_numThreads)
+		m_numThreads = m_max_numThreads;
+
+#ifdef MPI_SUPPORT
+	m_MPI_Barrier = 0;
+#endif
+	this->changeNumThreads(m_numThreads);
+}
+
+void Engine_Multithread::Reset()
+{
+	if (m_thread_group!=0) // prevent multiple invocations
+	{
+		ClearExtensions(); //prevent extensions from interfering with thread reset...
+
+		// stop the threads
+		//NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;
+		m_thread_group->interrupt_all();
+		m_thread_group->join_all(); // wait for termination
+		delete m_IterateBarrier;
+		m_IterateBarrier = 0;
+		delete m_startBarrier;
+		m_startBarrier = 0;
+		delete m_stopBarrier;
+		m_stopBarrier = 0;
+		delete m_thread_group;
+		m_thread_group = 0;
+	}
+
+	ENGINE_MULTITHREAD_BASE::Reset();
+}
+
+void Engine_Multithread::changeNumThreads(unsigned int numThreads)
+{
+	if (m_thread_group!=0)
+	{
+		m_thread_group->interrupt_all();
+		m_thread_group->join_all(); // wait for termination
+		delete m_thread_group;
+		m_thread_group = 0;
+		//m_stopThreads = false;
+	}
+
+	m_numThreads = numThreads;
+
+	if (g_settings.GetVerboseLevel()>0)
+		cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
 
 	vector<unsigned int> m_Start_Lines;
 	vector<unsigned int> m_Stop_Lines;
 	m_Op_MT->CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines );
 
-	if (g_settings.GetVerboseLevel()>0)
+	if (m_IterateBarrier!=0)
-		cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
+		delete m_IterateBarrier;
+	// make sure all threads are waiting
 	m_IterateBarrier = new boost::barrier(m_numThreads); // numThread workers
 
+	if (m_startBarrier!=0)
+		delete m_startBarrier;
 	m_startBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
+
+	if (m_stopBarrier!=0)
+		delete m_stopBarrier;
 	m_stopBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
 #ifdef MPI_SUPPORT
 	m_MPI_Barrier = 0;
 #endif
 
+	m_thread_group = new boost::thread_group();
 	for (unsigned int n=0; n<m_numThreads; n++)
 	{
 		unsigned int start = m_Start_Lines.at(n);
@@ -130,50 +194,44 @@ void Engine_Multithread::Init()
 				cout << stop-start+1 << ";";
 //		NS_Engine_Multithread::DBG().cout() << "###DEBUG## Thread " << n << ": start=" << start << " stop=" << stop  << " stop_h=" << stop_h << std::endl;
 		boost::thread *t = new boost::thread( NS_Engine_Multithread::thread(this,start,stop,stop_h,n) );
-		m_thread_group.add_thread( t );
+		m_thread_group->add_thread( t );
 	}
 
 	for (size_t n=0; n<m_Eng_exts.size(); ++n)
 		m_Eng_exts.at(n)->SetNumberOfThreads(m_numThreads);
 }
 
-void Engine_Multithread::Reset()
-{
-	if (!m_stopThreads) // prevent multiple invocations
-	{
-		ClearExtensions(); //prevent extensions from interfering with thread reset...
-
-		// stop the threads
-		//NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;
-		m_iterTS = 1;
-		m_startBarrier->wait(); // start the threads
-		m_stopThreads = true;
-		m_stopBarrier->wait(); // wait for the threads to finish
-		m_thread_group.join_all(); // wait for termination
-		delete m_IterateBarrier;
-		m_IterateBarrier = 0;
-		delete m_startBarrier;
-		m_startBarrier = 0;
-		delete m_stopBarrier;
-		m_stopBarrier = 0;
-	}
-
-	ENGINE_MULTITHREAD_BASE::Reset();
-}
-
 bool Engine_Multithread::IterateTS(unsigned int iterTS)
 {
 	m_iterTS = iterTS;
 
-	//cout << "bool Engine_Multithread::IterateTS(): starting threads ...";
+	//cerr << "bool Engine_Multithread::IterateTS(): starting threads ...";
 	m_startBarrier->wait(); // start the threads
 
-	//cout << "... threads started";
+	//cerr << "... threads started" << endl;
 
 	m_stopBarrier->wait(); // wait for the threads to finish <iterTS> time steps
+
 	return true;
 }
 
+void Engine_Multithread::NextInterval(float curr_speed)
+{
+	ENGINE_MULTITHREAD_BASE::NextInterval(curr_speed);
+	if (!m_opt_speed) return;
+	if (curr_speed<m_last_speed)
+	{
+		this->changeNumThreads(m_numThreads-1);
+		cout << "Multithreaded Engine: Best performance found using " << m_numThreads << " threads." << std::endl;
+		m_opt_speed = false;
+	}
+	else if (m_numThreads<m_max_numThreads)
+	{
+		m_last_speed = curr_speed;
+		this->changeNumThreads(m_numThreads+1);
+	}
+}
+
 void Engine_Multithread::DoPreVoltageUpdates(int threadID)
 {
 	//execute extensions in reverse order -> highest priority gets access to the voltages last
@@ -269,6 +327,12 @@ void thread::operator()()
 		m_enginePtr->m_startBarrier->wait();
 		//cout << "Thread " << boost::this_thread::get_id() << " waiting... started." << endl;
 
+		if (m_enginePtr->m_stopThreads)
+		{
+			//DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") stop!." << endl;
+			return;
+		}
+
 		DEBUG_TIME( Timer timer1 );
 
 		for (unsigned int iter=0; iter<m_enginePtr->m_iterTS; ++iter)

FDTD/engine_multithread.h

@@ -90,6 +90,7 @@ public:
 	virtual void setNumThreads( unsigned int numThreads );
 	virtual void Init();
 	virtual void Reset();
+	virtual void NextInterval(float curr_speed);
 
 	//! Iterate \a iterTS number of timesteps
 	virtual bool IterateTS(unsigned int iterTS);
@@ -104,13 +105,17 @@ public:
 
 protected:
 	Engine_Multithread(const Operator_Multithread* op);
+	void changeNumThreads(unsigned int numThreads);
 	const Operator_Multithread* m_Op_MT;
-	boost::thread_group m_thread_group;
+	boost::thread_group *m_thread_group;
 	boost::barrier *m_startBarrier, *m_stopBarrier;
 	boost::barrier *m_IterateBarrier;
 	volatile unsigned int m_iterTS;
 	unsigned int m_numThreads; //!< number of worker threads
+	unsigned int m_max_numThreads; //!< max. number of worker threads
 	volatile bool m_stopThreads;
+	bool m_opt_speed;
+	float m_last_speed;
 
 #ifdef MPI_SUPPORT
 	/*! Workaround needed for subgridding scheme... (see Engine_CylinderMultiGrid)

FDTD/engine_sse.cpp

@@ -91,16 +91,37 @@ void Engine_sse::UpdateVoltages(unsigned int startX, unsigned int numX)
 			for (pos[2]=1; pos[2]<numVectors; ++pos[2])
 			{
 				// x-polarization
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[0][pos[0]][pos[1]][pos[2]].v * ( f4_curr[2][pos[0]][pos[1]][pos[2]].v - f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v - f4_curr[1][pos[0]][pos[1]][pos[2]].v + f4_curr[1][pos[0]][pos[1]][pos[2]-1].v );
+				    Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[0][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[2][pos[0]][pos[1]         ][pos[2]].v -
+				        f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v -
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]].v +
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]-1].v
+				    );
 
 				// y-polarization
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[1][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[1][pos[0]][pos[1]][pos[2]].v * ( f4_curr[0][pos[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]-1].v - f4_curr[2][pos[0]][pos[1]][pos[2]].v + f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]].v);
+				    Op->f4_vv[1][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[1][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]  ].v -
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]-1].v -
+				        f4_curr[2][pos[0]         ][pos[1]][pos[2]  ].v +
+				        f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]  ].v
+				    );
 
 				// z-polarization
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[2][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[2][pos[0]][pos[1]][pos[2]].v * ( f4_curr[1][pos[0]][pos[1]][pos[2]].v - f4_curr[1][pos[0]-shift[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]].v + f4_curr[0][pos[0]][pos[1]-shift[1]][pos[2]].v);
+				    Op->f4_vv[2][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[2][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[1][pos[0]         ][pos[1]         ][pos[2]].v -
+				        f4_curr[1][pos[0]-shift[0]][pos[1]         ][pos[2]].v -
+				        f4_curr[0][pos[0]         ][pos[1]         ][pos[2]].v +
+				        f4_curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = 0
@@ -109,20 +130,41 @@ void Engine_sse::UpdateVoltages(unsigned int startX, unsigned int numX)
 			temp.f[1] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[2];
-			f4_volt[0][pos[0]][pos[1]][0].v *= Op->f4_vv[0][pos[0]][pos[1]][0].v;
+			f4_volt[0][pos[0]][pos[1]][0].v *=
-			f4_volt[0][pos[0]][pos[1]][0].v += Op->f4_vi[0][pos[0]][pos[1]][0].v * ( f4_curr[2][pos[0]][pos[1]][0].v - f4_curr[2][pos[0]][pos[1]-shift[1]][0].v - f4_curr[1][pos[0]][pos[1]][0].v + temp.v );
+			    Op->f4_vv[0][pos[0]][pos[1]][0].v;
+			f4_volt[0][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[0][pos[0]][pos[1]][0].v * (
+			        f4_curr[2][pos[0]][pos[1]         ][0].v -
+			        f4_curr[2][pos[0]][pos[1]-shift[1]][0].v -
+			        f4_curr[1][pos[0]][pos[1]         ][0].v +
+			        temp.v
+			    );
 
 			// y-polarization
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[2];
-			f4_volt[1][pos[0]][pos[1]][0].v *= Op->f4_vv[1][pos[0]][pos[1]][0].v;
+			f4_volt[1][pos[0]][pos[1]][0].v *=
-			f4_volt[1][pos[0]][pos[1]][0].v += Op->f4_vi[1][pos[0]][pos[1]][0].v * ( f4_curr[0][pos[0]][pos[1]][0].v - temp.v - f4_curr[2][pos[0]][pos[1]][0].v + f4_curr[2][pos[0]-shift[0]][pos[1]][0].v);
+			    Op->f4_vv[1][pos[0]][pos[1]][0].v;
+			f4_volt[1][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[1][pos[0]][pos[1]][0].v * (
+			        f4_curr[0][pos[0]         ][pos[1]][0].v -
+			        temp.v -
+			        f4_curr[2][pos[0]         ][pos[1]][0].v +
+			        f4_curr[2][pos[0]-shift[0]][pos[1]][0].v
+			    );
 
 			// z-polarization
-			f4_volt[2][pos[0]][pos[1]][0].v *= Op->f4_vv[2][pos[0]][pos[1]][0].v;
+			f4_volt[2][pos[0]][pos[1]][0].v *=
-			f4_volt[2][pos[0]][pos[1]][0].v += Op->f4_vi[2][pos[0]][pos[1]][0].v * ( f4_curr[1][pos[0]][pos[1]][0].v - f4_curr[1][pos[0]-shift[0]][pos[1]][0].v - f4_curr[0][pos[0]][pos[1]][0].v + f4_curr[0][pos[0]][pos[1]-shift[1]][0].v);
+			    Op->f4_vv[2][pos[0]][pos[1]][0].v;
+			f4_volt[2][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[2][pos[0]][pos[1]][0].v * (
+			        f4_curr[1][pos[0]         ][pos[1]         ][0].v -
+			        f4_curr[1][pos[0]-shift[0]][pos[1]         ][0].v -
+			        f4_curr[0][pos[0]         ][pos[1]         ][0].v +
+			        f4_curr[0][pos[0]         ][pos[1]-shift[1]][0].v
+			    );
 		}
 		++pos[0];
 	}
@@ -141,16 +183,37 @@ void Engine_sse::UpdateCurrents(unsigned int startX, unsigned int numX)
 			for (pos[2]=0; pos[2]<numVectors-1; ++pos[2])
 			{
 				// x-pol
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[0][pos[0]][pos[1]][pos[2]].v * ( f4_volt[2][pos[0]][pos[1]][pos[2]].v - f4_volt[2][pos[0]][pos[1]+1][pos[2]].v - f4_volt[1][pos[0]][pos[1]][pos[2]].v + f4_volt[1][pos[0]][pos[1]][pos[2]+1].v);
+				    Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[0][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[2][pos[0]][pos[1]  ][pos[2]  ].v -
+				        f4_volt[2][pos[0]][pos[1]+1][pos[2]  ].v -
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]  ].v +
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]+1].v
+				    );
 
 				// y-pol
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[1][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[1][pos[0]][pos[1]][pos[2]].v * ( f4_volt[0][pos[0]][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]+1].v - f4_volt[2][pos[0]][pos[1]][pos[2]].v + f4_volt[2][pos[0]+1][pos[1]][pos[2]].v);
+				    Op->f4_ii[1][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[1][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]  ].v -
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]+1].v -
+				        f4_volt[2][pos[0]  ][pos[1]][pos[2]  ].v +
+				        f4_volt[2][pos[0]+1][pos[1]][pos[2]  ].v
+				    );
 
 				// z-pol
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[2][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[2][pos[0]][pos[1]][pos[2]].v * ( f4_volt[1][pos[0]][pos[1]][pos[2]].v - f4_volt[1][pos[0]+1][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]].v + f4_volt[0][pos[0]][pos[1]+1][pos[2]].v);
+				    Op->f4_ii[2][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[2][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[1][pos[0]  ][pos[1]  ][pos[2]].v -
+				        f4_volt[1][pos[0]+1][pos[1]  ][pos[2]].v -
+				        f4_volt[0][pos[0]  ][pos[1]  ][pos[2]].v +
+				        f4_volt[0][pos[0]  ][pos[1]+1][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = numVectors-1
@@ -159,20 +222,41 @@ void Engine_sse::UpdateCurrents(unsigned int startX, unsigned int numX)
 			temp.f[1] = f4_volt[1][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[1][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[0][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[0][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[2][pos[0]][pos[1]][numVectors-1].v - f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v - f4_volt[1][pos[0]][pos[1]][numVectors-1].v + temp.v);
+			    Op->f4_ii[0][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[0][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[2][pos[0]][pos[1]  ][numVectors-1].v -
+			        f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v -
+			        f4_volt[1][pos[0]][pos[1]  ][numVectors-1].v +
+			        temp.v
+			    );
 
 			// y-pol
 			temp.f[0] = f4_volt[0][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[0][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[0][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[1][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[1][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[0][pos[0]][pos[1]][numVectors-1].v - temp.v - f4_volt[2][pos[0]][pos[1]][numVectors-1].v + f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v);
+			    Op->f4_ii[1][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[1][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[0][pos[0]  ][pos[1]][numVectors-1].v -
+			        temp.v -
+			        f4_volt[2][pos[0]  ][pos[1]][numVectors-1].v +
+			        f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v
+			    );
 
 			// z-pol
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[2][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[2][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[1][pos[0]][pos[1]][numVectors-1].v - f4_volt[1][pos[0]+1][pos[1]][numVectors-1].v - f4_volt[0][pos[0]][pos[1]][numVectors-1].v + f4_volt[0][pos[0]][pos[1]+1][numVectors-1].v);
+			    Op->f4_ii[2][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[2][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[1][pos[0]  ][pos[1]  ][numVectors-1].v -
+			        f4_volt[1][pos[0]+1][pos[1]  ][numVectors-1].v -
+			        f4_volt[0][pos[0]  ][pos[1]  ][numVectors-1].v +
+			        f4_volt[0][pos[0]  ][pos[1]+1][numVectors-1].v
+			    );
 		}
 		++pos[0];
 	}

FDTD/engine_sse_compressed.cpp

@@ -55,47 +55,93 @@ void Engine_SSE_Compressed::UpdateVoltages(unsigned int startX, unsigned int num
 			{
 				index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
 				// x-polarization
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[0][index].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[0][index].v * ( f4_curr[2][pos[0]][pos[1]][pos[2]].v - f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v - f4_curr[1][pos[0]][pos[1]][pos[2]].v + f4_curr[1][pos[0]][pos[1]][pos[2]-1].v );
+				    Op->f4_vv_Compressed[0][index].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[0][index].v * (
+				        f4_curr[2][pos[0]][pos[1]         ][pos[2]  ].v -
+				        f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]  ].v -
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]  ].v +
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]-1].v
+				    );
 
 				// y-polarization
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[1][index].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[1][index].v * ( f4_curr[0][pos[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]-1].v - f4_curr[2][pos[0]][pos[1]][pos[2]].v + f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]].v);
+				    Op->f4_vv_Compressed[1][index].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[1][index].v * (
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]  ].v -
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]-1].v -
+				        f4_curr[2][pos[0]         ][pos[1]][pos[2]  ].v +
+				        f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]  ].v
+				    );
 
 				// z-polarization
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[2][index].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v *=
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[2][index].v * ( f4_curr[1][pos[0]][pos[1]][pos[2]].v - f4_curr[1][pos[0]-shift[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]].v + f4_curr[0][pos[0]][pos[1]-shift[1]][pos[2]].v);
+				    Op->f4_vv_Compressed[2][index].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[2][index].v * (
+				        f4_curr[1][pos[0]         ][pos[1]]         [pos[2]].v -
+				        f4_curr[1][pos[0]-shift[0]][pos[1]]         [pos[2]].v -
+				        f4_curr[0][pos[0]         ][pos[1]]         [pos[2]].v +
+				        f4_curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = 0
 			// x-polarization
 			index = Op->m_Op_index[pos[0]][pos[1]][0];
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_slli_si128( (__m128i)f4_curr[1][pos[0]][pos[1]][numVectors-1].v, 4 );
+			temp.v = (__m128)_mm_slli_si128(
+			             (__m128i)f4_curr[1][pos[0]][pos[1]][numVectors-1].v, 4
+			         );
 #else
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[2];
 #endif
-			f4_volt[0][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[0][index].v;
+			f4_volt[0][pos[0]][pos[1]][0].v *=
-			f4_volt[0][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[0][index].v * ( f4_curr[2][pos[0]][pos[1]][0].v - f4_curr[2][pos[0]][pos[1]-shift[1]][0].v - f4_curr[1][pos[0]][pos[1]][0].v + temp.v );
+			    Op->f4_vv_Compressed[0][index].v;
+			f4_volt[0][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[0][index].v * (
+			        f4_curr[2][pos[0]][pos[1]         ][0].v -
+			        f4_curr[2][pos[0]][pos[1]-shift[1]][0].v -
+			        f4_curr[1][pos[0]][pos[1]         ][0].v +
+			        temp.v
+			    );
 
 			// y-polarization
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_slli_si128( (__m128i)f4_curr[0][pos[0]][pos[1]][numVectors-1].v, 4 );
+			temp.v = (__m128)_mm_slli_si128(
+			             (__m128i)f4_curr[0][pos[0]][pos[1]][numVectors-1].v, 4
+			         );
 #else
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[2];
 #endif
-			f4_volt[1][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[1][index].v;
+			f4_volt[1][pos[0]][pos[1]][0].v *=
-			f4_volt[1][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[1][index].v * ( f4_curr[0][pos[0]][pos[1]][0].v - temp.v - f4_curr[2][pos[0]][pos[1]][0].v + f4_curr[2][pos[0]-shift[0]][pos[1]][0].v);
+			    Op->f4_vv_Compressed[1][index].v;
+			f4_volt[1][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[1][index].v * (
+			        f4_curr[0][pos[0]         ][pos[1]][0].v -
+			        temp.v -
+			        f4_curr[2][pos[0]         ][pos[1]][0].v +
+			        f4_curr[2][pos[0]-shift[0]][pos[1]][0].v
+			    );
 
 			// z-polarization
-			f4_volt[2][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[2][index].v;
+			f4_volt[2][pos[0]][pos[1]][0].v *=
-			f4_volt[2][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[2][index].v * ( f4_curr[1][pos[0]][pos[1]][0].v - f4_curr[1][pos[0]-shift[0]][pos[1]][0].v - f4_curr[0][pos[0]][pos[1]][0].v + f4_curr[0][pos[0]][pos[1]-shift[1]][0].v);
+			    Op->f4_vv_Compressed[2][index].v;
+			f4_volt[2][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[2][index].v * (
+			        f4_curr[1][pos[0]         ][pos[1]         ][0].v -
+			        f4_curr[1][pos[0]-shift[0]][pos[1]         ][0].v -
+			        f4_curr[0][pos[0]         ][pos[1]         ][0].v +
+			        f4_curr[0][pos[0]         ][pos[1]-shift[1]][0].v
+			    );
 		}
 		++pos[0];
 	}
@@ -116,47 +162,93 @@ void Engine_SSE_Compressed::UpdateCurrents(unsigned int startX, unsigned int num
 			{
 				index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
 				// x-pol
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[0][index].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[0][index].v * ( f4_volt[2][pos[0]][pos[1]][pos[2]].v - f4_volt[2][pos[0]][pos[1]+1][pos[2]].v - f4_volt[1][pos[0]][pos[1]][pos[2]].v + f4_volt[1][pos[0]][pos[1]][pos[2]+1].v);
+				    Op->f4_ii_Compressed[0][index].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[0][index].v * (
+				        f4_volt[2][pos[0]][pos[1]  ][pos[2]  ].v -
+				        f4_volt[2][pos[0]][pos[1]+1][pos[2]  ].v -
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]  ].v +
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]+1].v
+				    );
 
 				// y-pol
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[1][index].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[1][index].v * ( f4_volt[0][pos[0]][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]+1].v - f4_volt[2][pos[0]][pos[1]][pos[2]].v + f4_volt[2][pos[0]+1][pos[1]][pos[2]].v);
+				    Op->f4_ii_Compressed[1][index].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[1][index].v * (
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]  ].v -
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]+1].v -
+				        f4_volt[2][pos[0]  ][pos[1]][pos[2]  ].v +
+				        f4_volt[2][pos[0]+1][pos[1]][pos[2]  ].v
+				    );
 
 				// z-pol
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[2][index].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v *=
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[2][index].v * ( f4_volt[1][pos[0]][pos[1]][pos[2]].v - f4_volt[1][pos[0]+1][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]].v + f4_volt[0][pos[0]][pos[1]+1][pos[2]].v);
+				    Op->f4_ii_Compressed[2][index].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[2][index].v * (
+				        f4_volt[1][pos[0]  ][pos[1]  ][pos[2]].v -
+				        f4_volt[1][pos[0]+1][pos[1]  ][pos[2]].v -
+				        f4_volt[0][pos[0]  ][pos[1]  ][pos[2]].v +
+				        f4_volt[0][pos[0]  ][pos[1]+1][pos[2]].v
+				    );
 			}
 
 			index = Op->m_Op_index[pos[0]][pos[1]][numVectors-1];
 			// for pos[2] = numVectors-1
 			// x-pol
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_srli_si128( (__m128i)f4_volt[1][pos[0]][pos[1]][0].v, 4 );
+			temp.v = (__m128)_mm_srli_si128(
+			             (__m128i)f4_volt[1][pos[0]][pos[1]][0].v, 4
+			         );
 #else
 			temp.f[0] = f4_volt[1][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[1][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[1][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
 #endif
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[0][index].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[0][index].v * ( f4_volt[2][pos[0]][pos[1]][numVectors-1].v - f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v - f4_volt[1][pos[0]][pos[1]][numVectors-1].v + temp.v);
+			    Op->f4_ii_Compressed[0][index].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[0][index].v * (
+			        f4_volt[2][pos[0]][pos[1]  ][numVectors-1].v -
+			        f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v -
+			        f4_volt[1][pos[0]][pos[1]  ][numVectors-1].v +
+			        temp.v
+			    );
 
 			// y-pol
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_srli_si128( (__m128i)f4_volt[0][pos[0]][pos[1]][0].v, 4 );
+			temp.v = (__m128)_mm_srli_si128(
+			             (__m128i)f4_volt[0][pos[0]][pos[1]][0].v, 4
+			         );
 #else
 			temp.f[0] = f4_volt[0][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[0][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[0][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
 #endif
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[1][index].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[1][index].v * ( f4_volt[0][pos[0]][pos[1]][numVectors-1].v - temp.v - f4_volt[2][pos[0]][pos[1]][numVectors-1].v + f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v);
+			    Op->f4_ii_Compressed[1][index].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[1][index].v * (
+			        f4_volt[0][pos[0]  ][pos[1]][numVectors-1].v -
+			        temp.v -
+			        f4_volt[2][pos[0]  ][pos[1]][numVectors-1].v +
+			        f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v
+			    );
 
 			// z-pol
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[2][index].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *=
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[2][index].v * ( f4_volt[1][pos[0]][pos[1]][numVectors-1].v - f4_volt[1][pos[0]+1][pos[1]][numVectors-1].v - f4_volt[0][pos[0]][pos[1]][numVectors-1].v + f4_volt[0][pos[0]][pos[1]+1][numVectors-1].v);
+			    Op->f4_ii_Compressed[2][index].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[2][index].v * (
+			        f4_volt[1][pos[0]  ][pos[1]  ][numVectors-1].v -
+			        f4_volt[1][pos[0]+1][pos[1]  ][numVectors-1].v -
+			        f4_volt[0][pos[0]  ][pos[1]  ][numVectors-1].v +
+			        f4_volt[0][pos[0]  ][pos[1]+1][numVectors-1].v
+			    );
 		}
 		++pos[0];
 	}

FDTD/extensions/CMakeLists.txt

@@ -24,6 +24,8 @@ set(SOURCES
   ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_steadystate.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_lumpedRLC.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_lumpedRLC.cpp
   PARENT_SCOPE
 )
 

FDTD/extensions/engine_ext_lumpedRLC.cpp

@@ -0,0 +1,195 @@
+/*
+*	Additional
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "engine_ext_lumpedRLC.h"
+#include "operator_ext_lumpedRLC.h"
+
+#include "FDTD/engine_sse.h"
+
+Engine_Ext_LumpedRLC::Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC* op_ext_RLC) : Engine_Extension(op_ext_RLC)
+{
+	// Local pointer of the operator.
+	m_Op_Ext_RLC = op_ext_RLC;
+
+	v_Vdn		= new FDTD_FLOAT*[3];
+	v_Jn		= new FDTD_FLOAT*[3];
+
+	// No additional allocations are required if there are no actual lumped elements.
+	if (!(m_Op_Ext_RLC->RLC_count))
+		return;
+
+	// Initialize ADE containers for currents and voltages
+	v_Il 		= new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+
+	for (unsigned int posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx)
+		v_Il[posIdx] 	= 0.0;
+
+	for (unsigned int k = 0 ; k < 3 ; k++)
+	{
+		v_Vdn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+		v_Jn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+
+		for (unsigned int posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx)
+		{
+			v_Jn[k][posIdx] = 0.0;
+			v_Vdn[k][posIdx] = 0.0;;
+		}
+	}
+
+}
+
+Engine_Ext_LumpedRLC::~Engine_Ext_LumpedRLC()
+{
+	// Only delete if values were allocated in the first place
+	if (m_Op_Ext_RLC->RLC_count)
+	{
+		delete[] v_Il;
+
+		for (unsigned int k = 0 ; k < 3 ; k++)
+		{
+			delete[] v_Vdn[k];
+			delete[] v_Jn[k];
+		}
+	}
+
+	delete[] v_Vdn;
+	delete[] v_Jn;
+
+	v_Il	= NULL;
+
+	v_Vdn	= NULL;
+	v_Jn	= NULL;
+
+	m_Op_Ext_RLC = NULL;
+
+
+}
+
+void Engine_Ext_LumpedRLC::DoPreVoltageUpdates()
+{
+	unsigned int **pos = m_Op_Ext_RLC->v_RLC_pos;
+	int *dir = m_Op_Ext_RLC->v_RLC_dir;
+
+	// Iterate Vd containers
+	FDTD_FLOAT	*v_temp;
+	v_temp = v_Vdn[2];
+	v_Vdn[2] = v_Vdn[1];
+	v_Vdn[1] = v_Vdn[0];
+	v_Vdn[0] = v_temp;
+
+	// In pre-process, only update the parallel inductor current:
+	for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+		v_Il[pIdx] += (m_Op_Ext_RLC->v_RLC_i2v[pIdx])*(m_Op_Ext_RLC->v_RLC_ilv[pIdx])*v_Vdn[1][pIdx];
+
+	return;
+}
+
+void Engine_Ext_LumpedRLC::Apply2Voltages()
+{
+	unsigned int **pos = m_Op_Ext_RLC->v_RLC_pos;
+	int *dir = m_Op_Ext_RLC->v_RLC_dir;
+
+	// Iterate J containers
+	FDTD_FLOAT	*v_temp;
+	v_temp = v_Jn[2];
+	v_Jn[2] = v_Jn[1];
+	v_Jn[1] = v_Jn[0];
+	v_Jn[0] = v_temp;
+
+
+	// Read engine calculated node voltage
+	switch (m_Eng->GetType())
+	{
+		case Engine::BASIC:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = m_Eng->Engine::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);
+
+			break;
+		}
+		case Engine::SSE:
+		{
+			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = eng_sse->Engine_sse::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);
+
+			break;
+		}
+		default:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = m_Eng->GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);;
+
+			break;
+		}
+	}
+
+	// Post process: Calculate node voltage with respect to the lumped RLC auxilliary quantity, J
+	for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+	{
+		// Calculate updated node voltage, with series and parallel additions
+		v_Vdn[0][pIdx] = 	(m_Op_Ext_RLC->v_RLC_vvd[pIdx])*(
+							v_Vdn[0][pIdx] - v_Il[pIdx]						// Addition for Parallel inductor
+							+
+							(m_Op_Ext_RLC->v_RLC_vv2[pIdx])*v_Vdn[2][pIdx]	// Vd[n-2] addition
+							+
+							(m_Op_Ext_RLC->v_RLC_vj1[pIdx])*v_Jn[1][pIdx]	// J[n-1] addition
+							+
+							(m_Op_Ext_RLC->v_RLC_vj2[pIdx])*v_Jn[2][pIdx]);	// J[n-2] addition
+
+		// Update J[0]
+		v_Jn[0][pIdx] =	(m_Op_Ext_RLC->v_RLC_ib0[pIdx])*(v_Vdn[0][pIdx] - v_Vdn[2][pIdx])
+						-
+						((m_Op_Ext_RLC->v_RLC_b1[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[1][pIdx]
+						-
+						((m_Op_Ext_RLC->v_RLC_b2[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[2][pIdx];
+	}
+
+
+	// Update node voltage
+	switch (m_Eng->GetType())
+	{
+		case Engine::BASIC:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				m_Eng->Engine::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+		case Engine::SSE:
+		{
+			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				eng_sse->Engine_sse::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+		default:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				m_Eng->SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+	}
+
+
+	return;
+}
+
+

FDTD/extensions/engine_ext_lumpedRLC.h

@@ -0,0 +1,54 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ENGINE_EXT_LUMPEDRLC_H
+#define ENGINE_EXT_LUMPEDRLC_H
+
+#include "engine_extension.h"
+#include "FDTD/engine.h"
+#include "FDTD/operator.h"
+
+class Operator_Ext_LumpedRLC;
+
+class Engine_Ext_LumpedRLC : public Engine_Extension
+{
+	friend class Operator_Ext_LumpedRLC;
+	friend class Operator;
+	friend class ContinuousStructure;
+
+public:
+
+	Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC *op_ext_RLC);
+	virtual ~Engine_Ext_LumpedRLC();
+
+	virtual void DoPreVoltageUpdates();
+	virtual void Apply2Voltages();
+
+protected:
+	Operator_Ext_LumpedRLC* m_Op_Ext_RLC;
+
+	// Auxilliary containers
+
+	// Array setup: volt_C_ADE[mesh_pos]
+	FDTD_FLOAT *v_Il;		// Container for current on inductor- Parallel RLC
+
+	FDTD_FLOAT **v_Vdn;		// Container for nodal vd at [n],[n-1],[n-2]
+	FDTD_FLOAT **v_Jn;		// Container for nodal J at [n],[n-1],[n-2]
+
+};
+
+#endif // ENGINE_EXT_LORENTZMATERIAL_H

FDTD/extensions/operator_ext_lumpedRLC.cpp

@@ -0,0 +1,514 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "../operator.h"
+#include "tools/array_ops.h"
+#include "tools/constants.h"
+//#include "cond_sheet_parameter.h"
+#include "tools/AdrOp.h"
+
+#include "operator_ext_lumpedRLC.h"
+#include "engine_ext_lumpedRLC.h"
+
+#include "CSPrimBox.h"
+#include "CSProperties.h"
+#include "CSPropLumpedElement.h"
+
+#define COPY_V2A(V,A) std::copy(V.begin(),V.end(),A)
+
+Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op) : Operator_Extension(op)
+{
+	// Parallel circuit coefficients
+	v_RLC_ilv = NULL;
+	v_RLC_i2v = NULL;
+
+	// Series circuit coefficients
+	v_RLC_vv2	= NULL;	// Coefficient for [n-2] time of Vd update in Vd equation
+	v_RLC_vj1	= NULL;	// Coefficient for [n-1] time of J update in Vd equation
+	v_RLC_vj2	= NULL;	// Coefficient for [n-2] time of J update in Vd equation
+	v_RLC_vvd	= NULL;	// Coefficient to multiply all Vd in the Vd update equation
+	v_RLC_ib0	= NULL;	// Inverse of beta_0
+	v_RLC_b1	= NULL;	// beta_1
+	v_RLC_b2	= NULL;	// beta_2
+
+	// Additional containers
+	v_RLC_dir = NULL;
+	v_RLC_pos = NULL;
+
+	RLC_count = 0;
+}
+
+Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext) : Operator_Extension(op,op_ext)
+{
+	// Parallel circuit coefficients
+	v_RLC_ilv = NULL;
+	v_RLC_i2v = NULL;
+
+	// Series circuit coefficients
+	v_RLC_vv2	= NULL;	// Coefficient for [n-2] time of Vd update in Vd equation
+	v_RLC_vj1	= NULL;	// Coefficient for [n-1] time of J update in Vd equation
+	v_RLC_vj2	= NULL;	// Coefficient for [n-2] time of J update in Vd equation
+	v_RLC_vvd	= NULL;	// Coefficient to multiply all Vd in the Vd update equation
+	v_RLC_ib0	= NULL;	// Inverse of beta_0
+	v_RLC_b1	= NULL;	// beta_1
+	v_RLC_b2	= NULL;	// beta_2
+
+	// Additional containers
+	v_RLC_dir = NULL;
+	v_RLC_pos = NULL;
+
+	RLC_count = 0;
+}
+
+Operator_Ext_LumpedRLC::~Operator_Ext_LumpedRLC()
+{
+	if (this->RLC_count)
+	{
+		// Parallel circuit coefficients
+		delete[] v_RLC_ilv;
+		delete[] v_RLC_i2v;
+
+		// Series circuit coefficients
+		delete[] v_RLC_vv2;
+		delete[] v_RLC_vj1;
+		delete[] v_RLC_vj2;
+		delete[] v_RLC_vvd;
+		delete[] v_RLC_ib0;
+		delete[] v_RLC_b1;
+		delete[] v_RLC_b2;
+
+		// Additional containers
+		delete[] v_RLC_dir;
+
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; dIdx++)
+			delete[] v_RLC_pos[dIdx];
+
+		delete[] v_RLC_pos;
+
+	}
+}
+
+Operator_Extension* Operator_Ext_LumpedRLC::Clone(Operator* op)
+{
+	if (dynamic_cast<Operator_Ext_LumpedRLC*>(this)==NULL)
+		return NULL;
+	return new Operator_Ext_LumpedRLC(op, this);
+}
+
+bool Operator_Ext_LumpedRLC::BuildExtension()
+{
+	double 	dT = m_Op->GetTimestep();
+
+	double	fMax = 	m_Op->GetExcitationSignal()->GetCenterFreq()
+					+
+					m_Op->GetExcitationSignal()->GetCutOffFreq();
+
+	unsigned int 	pos[] = {0,0,0};
+
+	vector<CSProperties*> cs_props;
+
+	int 			dir;
+	CSPropLumpedElement::LEtype lumpedType;
+
+	vector<unsigned int> 	v_pos[3];
+
+	vector<int>		v_dir;
+
+	vector<double>  v_ilv;
+	vector<double>	v_i2v;
+
+	vector<double>	v_vv2;
+	vector<double>	v_vj1;
+	vector<double>	v_vj2;
+	vector<double>	v_vvd;
+	vector<double>	v_ib0;
+	vector<double>	v_b1;
+	vector<double>	v_b2;
+
+	// Lumped RLC parameters
+	double R, L, C;
+
+	// clear all vectors to initialize them
+	for (unsigned int dIdx = 0 ; dIdx < 3 ; dIdx++)
+		v_pos[dIdx].clear();
+
+	v_dir.clear();
+
+	v_ilv.clear();
+	v_i2v.clear();
+
+	v_vv2.clear();
+	v_vj1.clear();
+	v_vj2.clear();
+	v_vvd.clear();
+	v_ib0.clear();
+	v_b1.clear();
+	v_b2.clear();
+
+	// Obtain from CSX (continuous structure) all the lumped RLC properties
+	// Properties are material properties, not the objects themselves
+	cs_props = m_Op->CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
+
+	// Iterate through various properties. In theory, there should be a property set per-
+	// primitive, as each lumped element should have it's own unique properties.
+	for(size_t n = 0 ; n < cs_props.size() ; ++n)
+	{
+		// Cast current property to lumped RLC property continuous structure properties
+		CSPropLumpedElement* cs_RLC_props = dynamic_cast<CSPropLumpedElement*>(cs_props.at(n));
+		if (cs_RLC_props==NULL)
+			return false; //sanity check: this should never happen!
+
+		// Store direction and type
+		dir = cs_RLC_props->GetDirection();
+		lumpedType = cs_RLC_props->GetLEtype();
+
+		//		if (lumpedType == LEtype::INVALID
+		if (lumpedType == CSPropLumpedElement::INVALID)
+		{
+			cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: RLCtype is invalid! considering as parallel. "
+					<< " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+			lumpedType = CSPropLumpedElement::PARALLEL;
+		}
+
+		// Extract R, L and C from property class
+		C = cs_RLC_props->GetCapacity();
+		if (C < 0)
+			C = NAN;
+		R = cs_RLC_props->GetResistance();
+		if (R < 0)
+			R = NAN;
+		L = cs_RLC_props->GetInductance();
+		if (L < 0)
+			L = NAN;
+
+		// Check that this is a lumped RLC
+		if (!(this->IsLElumpedRLC(cs_RLC_props)))
+			continue;
+
+		if ((dir < 0) || (dir > 2))
+		{
+			cerr << "Operator_Ext_LumpedRLC::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. "
+					<< " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+			continue;
+		}
+
+		// Initialize other two direction containers
+		int dir_p1 = (dir + 1) % 3;
+		int dir_p2 = (dir + 2) % 3;
+
+		// Now iterate through primitive(s). I still think there should be only one per-
+		// material definition, but maybe I'm wrong...
+		vector<CSPrimitives*> cs_RLC_prims = cs_RLC_props->GetAllPrimitives();
+
+		for (size_t boxIdx = 0 ; boxIdx < cs_RLC_prims.size() ; ++boxIdx)
+		{
+			CSPrimBox* cBox = dynamic_cast<CSPrimBox*>(cs_RLC_prims.at(boxIdx));
+
+			if (cBox)
+			{
+
+				// Get box start and stop positions
+				unsigned int 	uiStart[3],
+								uiStop[3];
+
+
+				// snap to the native coordinate system
+				int Snap_Dimension =
+						m_Op->SnapBox2Mesh(
+								cBox->GetStartCoord()->GetCoords(m_Op->m_MeshType),	// Start Coord
+								cBox->GetStopCoord()->GetCoords(m_Op->m_MeshType),	// Stop Coord
+								uiStart,	// Start Index
+								uiStop,		// Stop Index
+								false,		// Dual (doublet) Grid?
+								true);		// Full mesh?
+
+				// Verify that snapped dimension is correct
+				if (Snap_Dimension<=0)
+				{
+					if (Snap_Dimension>=-1)
+						cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC snapping failed! Dimension is: " << Snap_Dimension << " skipping. "
+								<< " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+					// Snap_Dimension == -2 means outside the simulation domain --> no special warning, but box probably marked as unused!
+					continue;
+				}
+
+				// Verify that in the direction of the current propagation, the size isn't zero.
+				if (uiStart[dir]==uiStop[dir])
+				{
+					cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC with zero (snapped) length is invalid! skipping. "
+							<< " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+					continue;
+				}
+
+				// Calculate number of cells per-direction
+				unsigned int 	Ncells_0 = uiStop[dir] - uiStart[dir],
+								Ncells_1 = uiStop[dir_p1] - uiStart[dir_p1] + 1,
+								Ncells_2 = uiStop[dir_p2] - uiStart[dir_p2] + 1;
+
+				// All cells in directions 1 and 2 are considered parallel connection
+				unsigned int Npar = Ncells_1*Ncells_2;
+
+				// Separate elements such that individual elements can be calculated.
+				double	dL = L*Npar/Ncells_0,
+						dR = R*Npar/Ncells_0,
+						dG = (1/R)*Ncells_0/Npar,
+						dC = C*Ncells_0/Npar;
+
+				double	ib0	= 2.0*dT*dC/(4.0*dL*dC + 2.0*dT*dR*dC + dT*dT),
+						b1	= (dT*dT - 4.0*dL*dC)/(dT*dC),
+						b2	= (4.0*dL*dC - 2.0*dT*dR*dC + dT*dT)/(2.0*dT*dC);
+
+				// Special case: If this is a parallel resonant circuit, and there is no
+				// parallel resistor, use zero conductivity. May be risky when low-loss
+				// simulations are involved
+				if (lumpedType == CSPropLumpedElement::PARALLEL)
+					if (R == 0.0)
+						dG = 0.0;
+
+				int iPos = 0;
+
+				double Zmin,Zcd_min;
+
+				// In the various positions, update the capacitors and "inverse" resistors
+				for (pos[dir] = uiStart[dir] ; pos[dir] < uiStop[dir] ; ++pos[dir])
+				{
+					for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1])
+					{
+						for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2])
+						{
+							iPos = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]);
+
+
+							// Separate to two different cases. Parallel and series
+							switch (lumpedType)
+							{
+								case CSPropLumpedElement::PARALLEL:
+									// Update capacitor either way.
+									if (dC > 0)
+										m_Op->EC_C[dir][iPos] = dC;
+									else
+										// This case takes the "natural" capacitor into account.
+										dC = m_Op->EC_C[dir][iPos];
+
+									v_i2v.push_back((dT/dC)/(1.0 + dT*dG/(2.0*dC)));
+
+									// Update conductivity
+									if (R >= 0)
+										m_Op->EC_G[dir][iPos] = dG;
+
+									// Update coefficients with respect to the parallel inductance
+									if (L > 0)
+										v_ilv.push_back(dT/dL);
+									else
+										v_ilv.push_back(0.0);
+
+									// Take into account the case that the "natural" capacitor is too small
+									// with respect to the inductor or the resistor, and add a warning.
+									if (dC == 0)
+									{
+										double Cd = m_Op->EC_C[dir][iPos];
+										Zmin = max(dR,2*PI*fMax*dL);
+										Zcd_min = 1.0/(2.0*PI*fMax*Cd);
+
+										// Check if the "parasitic" capcitance is not small enough
+										if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin)
+										{
+											Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT);
+											m_Op->EC_C[dir][iPos] = Cd;
+										}
+									}
+
+									v_vv2.push_back(0.0);
+									v_vj1.push_back(0.0);
+									v_vj2.push_back(0.0);
+									v_vvd.push_back(1.0);
+									v_ib0.push_back(0.0);
+									v_b1.push_back(0.0);
+									v_b2.push_back(0.0);
+
+									// Update with discrete component values of
+									m_Op->Calc_ECOperatorPos(dir,pos);
+
+									v_dir.push_back(dir);
+
+									break;
+
+								case CSPropLumpedElement::SERIES:
+									m_Op->EC_G[dir][iPos] = 0.0;
+
+									// is a series inductor, modeled separately.
+									FDTD_FLOAT Cd = m_Op->EC_C[dir][iPos];
+
+									// Calculate minimum impedance, at maximum frequency
+									Zmin = sqrt(pow(dR,2) + pow(2*PI*fMax*dL - 1.0/(dC*2*PI*fMax),2));
+									Zcd_min = 1.0/(2.0*PI*fMax*Cd);
+
+									// Check if the "parasitic" capcitance is not small enough
+									if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin)
+									{
+										Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT);
+										m_Op->EC_C[dir][iPos] = Cd;
+									}
+
+									// No contribution from parallel inductor
+									v_ilv.push_back(0.0);
+									v_i2v.push_back(0.0);
+
+									// Contributions from series resistor and inductor
+									v_vv2.push_back(0.5*dT*ib0/Cd);
+									v_vj1.push_back(0.5*dT*(b1*ib0 - 1.0)/Cd);
+									v_vj2.push_back(0.5*dT*b2*ib0/Cd);
+									v_vvd.push_back(1.0/(1.0 + 0.5*dT*ib0/Cd));
+									v_ib0.push_back(ib0);
+									v_b1.push_back(b1);
+									v_b2.push_back(b2);
+
+									m_Op->Calc_ECOperatorPos(dir,pos);
+
+									v_dir.push_back(dir);
+
+									break;
+							}
+
+							// Store position and direction
+							for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+								v_pos[dIdx].push_back(pos[dIdx]);
+
+						}
+					}
+				}
+
+				// Build metallic caps
+				if (cs_RLC_props->GetCaps())
+					for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1])
+					{
+						for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2])
+						{
+							pos[dir]=uiStart[dir];
+							if (pos[dir_p1]<uiStop[dir_p1])
+							{
+								m_Op->SetVV(dir_p1,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p1]);
+							}
+
+							if (pos[dir_p2]<uiStop[dir_p2])
+							{
+								m_Op->SetVV(dir_p2,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p2]);
+							}
+
+							pos[dir]=uiStop[dir];
+							if (pos[dir_p1]<uiStop[dir_p1])
+							{
+								m_Op->SetVV(dir_p1,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p1]);
+							}
+
+							if (pos[dir_p2]<uiStop[dir_p2])
+							{
+								m_Op->SetVV(dir_p2,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p2]);
+							}
+						}
+					}
+
+
+				// Mark as used
+				cBox->SetPrimitiveUsed(true);
+			}
+		}
+
+	}
+
+	// Start data storage
+	RLC_count = v_dir.size();
+
+	// values
+	if (RLC_count)
+	{
+		// Allocate space to all variables
+		v_RLC_dir 	= new int[RLC_count];
+
+		// Parallel circuit coefficients
+		v_RLC_ilv 	= new FDTD_FLOAT[RLC_count];
+		v_RLC_i2v 	= new FDTD_FLOAT[RLC_count];
+
+		// Series circuit coefficients
+		v_RLC_vv2 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vj1 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vj2 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vvd = new FDTD_FLOAT[RLC_count];
+		v_RLC_ib0 = new FDTD_FLOAT[RLC_count];
+		v_RLC_b1 = new FDTD_FLOAT[RLC_count];
+		v_RLC_b2 = new FDTD_FLOAT[RLC_count];
+
+		v_RLC_pos = new unsigned int*[3];
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+			v_RLC_pos[dIdx] = new unsigned int[RLC_count];
+
+		// Copy all vectors to arrays
+		COPY_V2A(v_dir, v_RLC_dir);
+
+		COPY_V2A(v_ilv, v_RLC_ilv);
+		COPY_V2A(v_i2v, v_RLC_i2v);
+
+		COPY_V2A(v_vv2,v_RLC_vv2);
+		COPY_V2A(v_vj1,v_RLC_vj1);
+		COPY_V2A(v_vj2,v_RLC_vj2);
+		COPY_V2A(v_vvd,v_RLC_vvd);
+		COPY_V2A(v_ib0,v_RLC_ib0);
+		COPY_V2A(v_b1,v_RLC_b1);
+		COPY_V2A(v_b2,v_RLC_b2);
+
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+			COPY_V2A(v_pos[dIdx],v_RLC_pos[dIdx]);
+	}
+
+	return true;
+}
+
+Engine_Extension* Operator_Ext_LumpedRLC::CreateEngineExtention()
+{
+	Engine_Ext_LumpedRLC* eng_ext_RLC = new Engine_Ext_LumpedRLC(this);
+	return eng_ext_RLC;
+}
+
+void Operator_Ext_LumpedRLC::ShowStat(ostream &ostr)  const
+{
+	Operator_Extension::ShowStat(ostr);
+	string On_Off[2] = {"Off", "On"};
+
+	ostr << "Active cells\t\t: " << RLC_count << endl;
+}
+
+bool Operator_Ext_LumpedRLC::IsLElumpedRLC(const CSPropLumpedElement* const p_prop)
+{
+	CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype();
+
+	double L = p_prop->GetInductance();
+
+	bool isParallelRLC = (lumpedType == CSPropLumpedElement::PARALLEL) && (L > 0.0);
+	bool isSeriesRLC = lumpedType == CSPropLumpedElement::SERIES;
+
+	// This needs to be something that isn't a parallel RC circuit to add data to this extension.
+	return isParallelRLC || isSeriesRLC;
+}
+

FDTD/extensions/operator_ext_lumpedRLC.h

@@ -0,0 +1,88 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPERATOR_EXT_LUMPEDRLC_H_
+#define OPERATOR_EXT_LUMPEDRLC_H_
+
+#include "vector"
+
+#include "FDTD/operator.h"
+#include "operator_extension.h"
+#include "operator_ext_cylinder.h"
+
+#include "engine_ext_lumpedRLC.h"
+
+#define LUMPED_RLC_Z_FACT 20.0
+
+class Operator_Ext_LumpedRLC : public Operator_Extension
+{
+	friend class Engine_Ext_LumpedRLC;
+public:
+	Operator_Ext_LumpedRLC(Operator* op);
+
+	virtual ~Operator_Ext_LumpedRLC();
+
+	virtual Operator_Extension* Clone(Operator* op);
+
+	virtual bool BuildExtension();
+
+	virtual Engine_Extension* CreateEngineExtention();
+
+	virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;}
+	virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;}
+	virtual bool IsMPISave() const {return true;}
+
+	virtual string GetExtensionName() const {return string("Series\\Parallel Lumped RLC load");}
+
+	virtual void ShowStat(ostream &ostr) const;
+
+	virtual bool IsLElumpedRLC(const CSPropLumpedElement* const p_prop);
+
+protected:
+	//! Copy constructor
+	Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext);
+
+	// ADE update coefficients, array setup: coeff[mesh_pos_index]
+
+	// Parallel circuit coefficients
+	FDTD_FLOAT		*v_RLC_ilv;
+	FDTD_FLOAT		*v_RLC_i2v;
+
+	// Series circuit coefficients
+	FDTD_FLOAT		*v_RLC_vv2;		// Coefficient for [n-2] time of Vd update in Vd equation
+	FDTD_FLOAT		*v_RLC_vj1;		// Coefficient for [n-1] time of J update in Vd equation
+	FDTD_FLOAT		*v_RLC_vj2;		// Coefficient for [n-2] time of J update in Vd equation
+	FDTD_FLOAT		*v_RLC_vvd;		// Coefficient to multiply all Vd in the Vd update equation
+	FDTD_FLOAT		*v_RLC_ib0;		// Inverse of beta_0
+	FDTD_FLOAT		*v_RLC_b1;		// beta_1
+	FDTD_FLOAT		*v_RLC_b2;		// beta_2
+
+	// Additional containers
+	int				*v_RLC_dir;
+	unsigned int			**v_RLC_pos;
+
+	// Vector length indicator
+	unsigned int			RLC_count;
+
+
+
+
+};
+
+
+
+#endif /* OPERATOR_EXT_LUMPEDRLC_H_ */

FDTD/extensions/operator_ext_upml.cpp

@@ -404,7 +404,9 @@ bool Operator_Ext_UPML::BuildExtension()
 					CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i);
 					nP = (n+1)%3;
 					nPP = (n+2)%3;
-					if ((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0)
+					// if eff_Mat[1] > 1e3 assume a metal and disable PML to continue a signal layer
+					// sigma > 1000 S/m is "very lossy" already...
+					if (((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0) && (eff_Mat[1]<1e3))
 					{
 						//check if pos is on PEC
 						if ( (m_Op->GetVV(n,pos[0],pos[1],pos[2]) + m_Op->GetVI(n,pos[0],pos[1],pos[2])) != 0 )

FDTD/operator.cpp

@@ -379,7 +379,6 @@ int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned in
 	return ret;
 }
 
-
 Grid_Path Operator::FindPath(double start[], double stop[])
 {
 	Grid_Path path;
@@ -790,7 +789,7 @@ void Operator::DumpMaterial2File(string filename)
 	delete vtk_Writer;
 }
 
- bool Operator::SetupCSXGrid(CSRectGrid* grid)
+bool Operator::SetupCSXGrid(CSRectGrid* grid)
  {
 	 for (int n=0; n<3; ++n)
 	 {
@@ -1536,11 +1535,15 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, v
 bool Operator::Calc_LumpedElements()
 {
 	vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
+
 	for (size_t i=0;i<props.size();++i)
 	{
+
 		CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
+
 		if (PLE==NULL)
 			return false; //sanity check: this should never happen!
+
 		vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
 		for (size_t bn=0;bn<prims.size();++bn)
 		{
@@ -1555,6 +1558,11 @@ bool Operator::Calc_LumpedElements()
 				if (R<0)
 					R = NAN;
 
+				// If this is not a parallel RC, skip this.
+				if (!(this->IsLEparRC(PLE)))
+					continue;
+
+
 				if ((std::isnan(R)) && (std::isnan(C)))
 				{
 					cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
@@ -1703,6 +1711,18 @@ bool Operator::Calc_LumpedElements()
 	return true;
 }
 
+bool Operator::IsLEparRC(const CSPropLumpedElement* const p_prop)
+{
+	CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype();
+
+	double L = p_prop->GetInductance();
+
+	bool IsParallelRC = (lumpedType == CSPropLumpedElement::PARALLEL) && !(L > 0.0);
+
+	// This needs to be something that isn't a parallel RC circuit to add data to this extension.
+	return IsParallelRC;
+}
+
 void Operator::Init_EC()
 {
 	for (int n=0; n<3; ++n)

FDTD/operator.h

@@ -33,12 +33,16 @@ class Operator : public Operator_Base
 {
 	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_ConductingSheet; //we need to find a way around this... friend class Operator_Extension only would be nice
+	friend class Operator_Ext_ConductingSheet; 	// 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_Excitation;
 	friend class Operator_Ext_UPML;
 	friend class Operator_Ext_Cylinder;
+	friend class Operator_Ext_LumpedRLC;		// Gadi: I now know why the two previous remarks are here.
+
+	// So apparaently I have to use functionality from operator
+	// in my "lumpedRLC" class. This is ugly...
 public:
 	enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
 
@@ -244,6 +248,9 @@ protected:
 	//! Calculate and setup lumped elements
 	virtual bool Calc_LumpedElements();
 
+	//! Condition to determine if this is a lumped RC, to invoke Calc_LumpedElements
+	virtual bool IsLEparRC(const CSPropLumpedElement* const p_prop);
+
 	//! Store the size of the applied boundary conditions
 	int m_BC_Size[6];
 

FDTD/operator_cylinder.cpp

@@ -46,7 +46,7 @@ Operator_Cylinder::~Operator_Cylinder()
 Engine* Operator_Cylinder::CreateEngine()
 {
 	//! create a special cylindrical-engine
-	m_Engine = Engine_Cylinder::New(this, m_numThreads);
+	m_Engine = Engine_Cylinder::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 

FDTD/operator_cylindermultigrid.cpp

@@ -51,7 +51,7 @@ Operator_CylinderMultiGrid* Operator_CylinderMultiGrid::New(vector<double> Split
 
 Engine* Operator_CylinderMultiGrid::CreateEngine()
 {
-	m_Engine = Engine_CylinderMultiGrid::New(this,m_numThreads);
+	m_Engine = Engine_CylinderMultiGrid::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 

FDTD/operator_multithread.cpp

@@ -36,11 +36,12 @@ Operator_Multithread::~Operator_Multithread()
 void Operator_Multithread::setNumThreads( unsigned int numThreads )
 {
 	m_numThreads = numThreads;
+	m_orig_numThreads = numThreads;
 }
 
 Engine* Operator_Multithread::CreateEngine()
 {
-	m_Engine = Engine_Multithread::New(this,m_numThreads);
+	m_Engine = Engine_Multithread::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 
@@ -106,7 +107,7 @@ void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<u
 
 int Operator_Multithread::CalcECOperator( DebugFlags debugFlags )
 {
-	if (m_numThreads == 0)
+	if ((m_numThreads == 0) || (m_numThreads > boost::thread::hardware_concurrency()))
 		m_numThreads = boost::thread::hardware_concurrency();
 
 	vector<unsigned int> m_Start_Lines;

FDTD/operator_multithread.h

@@ -64,6 +64,7 @@ protected:
 
 	boost::thread_group m_thread_group;
 	unsigned int m_numThreads; // number of worker threads
+	unsigned int m_orig_numThreads;
 
 	//! Calculate the start/stop lines for the multithreading operator and engine.
 	/*!

matlab/setup.m

@@ -18,15 +18,22 @@ if isOctave()
     disp('compiling oct files')
     fflush(stdout);
     if isunix
-        [res, fn_so] = unix('find /usr/lib     -name libhdf5.so');
+        dylib_extension = 'so';
-        [res, fn_h]  = unix('find /usr/include -name hdf5.h | sort -r | head -1');
+        search_path = "/usr";
+        if ismac
+            dylib_extension = 'dylib';
+            [res, search_path] = unix('brew --prefix');
+            search_path = strtrim(search_path);
+        end
+        [res, fn_so] = unix(['find ' search_path '/lib     -name libhdf5.' dylib_extension]);
+        [res, fn_h]  = unix(['find ' search_path '/include -name hdf5.h | grep -v opencv | sort -r | head -1']);
         if length(fn_so)>0 && length(fn_h)>0
             [hdf5lib_dir, hdf5lib_fn, ext] = fileparts(fn_so);
             disp(["HDF5 library path found at: " hdf5lib_dir])
 
             [hdf5inc_dir, hdf5inc_fn, ext] = fileparts(fn_h);
             disp(["HDF5 include path found at: " hdf5inc_dir])
-            mkoctfile("h5readatt_octave.cc", ["-L" hdf5lib_dir], ["-I" hdf5inc_dir], "-L hdf5")
+            mkoctfile("h5readatt_octave.cc", ["-L" hdf5lib_dir], ["-I" hdf5inc_dir], "-lhdf5")
         else
             mkoctfile -lhdf5 h5readatt_octave.cc
         end

nf2ff/CMakeLists.txt

@@ -33,6 +33,7 @@ set(HEADERS
 
 add_library( nf2ff SHARED ${SOURCES})
 set_target_properties(nf2ff PROPERTIES VERSION ${LIB_VERSION_STRING} SOVERSION ${LIB_VERSION_MAJOR})
+set_target_properties(nf2ff PROPERTIES CXX_STANDARD 11)
 if (WIN32)
     target_compile_definitions(nf2ff PRIVATE -DBUILD_NF2FF_LIB )
 endif (WIN32)

nf2ff/main.cpp

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2023 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -25,7 +25,7 @@ int main(int argc, char *argv[])
 {
 	cout << " ---------------------------------------------------------------------- " << endl;
 	cout << " | nf2ff, near-field to far-field transformation for openEMS " << endl;
-	cout << " | (C) 2012-2014 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license" << endl;
+	cout << " | (C) 2012-2023 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license" << endl;
 	cout << " ---------------------------------------------------------------------- " << endl;
 
 	if (argc<=1)

openems.cpp

@@ -20,6 +20,7 @@
 #include <iostream>
 #include <fstream>
 #include "tools/array_ops.h"
+#include "tools/signal.h"
 #include "tools/useful.h"
 #include "FDTD/operator_cylinder.h"
 #include "FDTD/operator_cylindermultigrid.h"
@@ -30,6 +31,7 @@
 #include "FDTD/extensions/operator_ext_mur_abc.h"
 #include "FDTD/extensions/operator_ext_upml.h"
 #include "FDTD/extensions/operator_ext_lorentzmaterial.h"
+#include "FDTD/extensions/operator_ext_lumpedRLC.h"
 #include "FDTD/extensions/operator_ext_conductingsheet.h"
 #include "FDTD/extensions/operator_ext_steadystate.h"
 #include "FDTD/extensions/engine_ext_steadystate.h"
@@ -245,6 +247,13 @@ bool openEMS::parseCommandLineArgument( const char *argv )
 	return false;
 }
 
+void openEMS::SetNumberOfThreads(int val)
+{
+	if ((val<0) || (val>boost::thread::hardware_concurrency()))
+		val = boost::thread::hardware_concurrency();
+	m_engine_numThreads = val;
+}
+
 string openEMS::GetExtLibsInfo(string prefix)
 {
 	stringstream str;
@@ -285,8 +294,8 @@ void openEMS::WelcomeScreen()
 #endif
 
 	cout << " ---------------------------------------------------------------------- " << endl;
-	cout << " | openEMS " << bits << " -- version " GIT_VERSION << endl;
+	cout << " | openEMS " << bits << " -- version " << GIT_VERSION << endl;
-	cout << " | (C) 2010-2018 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license"   << endl;
+	cout << " | (C) 2010-2023 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license"   << endl;
 	cout << " ---------------------------------------------------------------------- " << endl;
 	cout << openEMS::GetExtLibsInfo("\t") << endl;
 }
@@ -873,6 +882,12 @@ void openEMS::SetStepExcite(double f_max)
 	m_Exc->SetupStepExcite(f_max);
 }
 
+void openEMS::SetCustomExcite(std::string str, double f0, double fmax)
+{
+	this->InitExcitation();
+	m_Exc->SetupCustomExcite(str, f0, fmax);
+}
+
 Excitation* openEMS::InitExcitation()
 {
 	delete m_Exc;
@@ -891,14 +906,18 @@ int openEMS::SetupFDTD()
 	timeval startTime;
 	gettimeofday(&startTime,NULL);
 
+	Signal::SetupHandlerForSIGINT(SIGNAL_EXIT_FORCE);
+
 	if (m_CSX==NULL)
 	{
 		cerr << "openEMS::SetupFDTD: Error: CSXCAD is not set!" << endl;
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 3;
 	}
 	if (m_CSX==NULL)
 	{
 		cerr << "openEMS::SetupFDTD: Error: CSXCAD is not set!" << endl;
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 3;
 	}
 	std::string ec = m_CSX->Update();
@@ -919,7 +938,10 @@ int openEMS::SetupFDTD()
 
 	//*************** setup operator ************//
 	if (SetupOperator()==false)
+	{
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 2;
+	}
 
 	// default material averaging is quarter cell averaging
 	FDTD_Op->SetQuarterCellMaterialAvg();
@@ -934,6 +956,7 @@ int openEMS::SetupFDTD()
 	if (m_Exc==NULL)
 	{
 		cerr << "openEMS::SetupFDTD: Error, excitation is not defined! Abort!" << endl;
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 3;
 	}
 
@@ -942,7 +965,11 @@ int openEMS::SetupFDTD()
 	if (!CylinderCoords)
 		FDTD_Op->AddExtension(new Operator_Ext_TFSF(FDTD_Op));
 
-	if (FDTD_Op->SetGeometryCSX(m_CSX)==false) return(2);
+	if (FDTD_Op->SetGeometryCSX(m_CSX)==false)
+	{
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
+		return(2);
+	}
 
 	SetupBoundaryConditions();
 
@@ -988,6 +1015,9 @@ int openEMS::SetupFDTD()
 		FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op));
 	if (m_CSX->GetQtyPropertyType(CSProperties::CONDUCTINGSHEET)>0)
 		FDTD_Op->AddExtension(new Operator_Ext_ConductingSheet(FDTD_Op, m_Exc->GetMaxFreq()));
+	if (m_CSX->GetQtyPropertyType(CSProperties::LUMPED_ELEMENT)>0)
+		FDTD_Op->AddExtension(new Operator_Ext_LumpedRLC(FDTD_Op));
+
 
 	//check all properties to request material storage during operator creation...
 	SetupMaterialStorages();
@@ -1054,6 +1084,7 @@ int openEMS::SetupFDTD()
 	if (m_no_simulation)
 	{
 		// simulation was disabled (to generate debug output only)
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 1;
 	}
 
@@ -1068,7 +1099,10 @@ int openEMS::SetupFDTD()
 
 	//setup all processing classes
 	if (SetupProcessing()==false)
+	{
+		Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 		return 2;
+	}
 
 	// Cleanup all unused material storages...
 	FDTD_Op->CleanupMaterialStorage();
@@ -1082,6 +1116,7 @@ int openEMS::SetupFDTD()
 		PA->DumpBoxes2File("box_dump_");
 	}
 
+	Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 	return 0;
 }
 
@@ -1107,12 +1142,19 @@ bool openEMS::CheckAbortCond()
 	if (m_Abort) //abort was set externally
 		return true;
 
+	//check whether SIGINT is received
+	if (Signal::ReceivedSIGINT())
+	{
+		cerr << "openEMS::CheckAbortCond(): Received SIGINT, aborting simulation gracefully..." << endl;
+		return true;
+	}
+
 	//check whether the file "ABORT" exist in current working directory
 	ifstream ifile("ABORT");
 	if (ifile)
 	{
 		ifile.close();
-		cerr << "openEMS::CheckAbortCond(): Found file \"ABORT\", aborting simulation..." << endl;
+		cerr << "openEMS::CheckAbortCond(): Found file \"ABORT\", aborting simulation gracefully..." << endl;
 		return true;
 	}
 
@@ -1123,6 +1165,8 @@ void openEMS::RunFDTD()
 {
 	cout << "Running FDTD engine... this may take a while... grab a cup of coffee?!?" << endl;
 
+	Signal::SetupHandlerForSIGINT(SIGNAL_EXIT_GRACEFUL);
+
 	//special handling of a field processing, needed to realize the end criteria...
 	ProcessFields* ProcField = new ProcessFields(NewEngineInterface());
 	PA->AddProcessing(ProcField);
@@ -1203,6 +1247,7 @@ void openEMS::RunFDTD()
 
 			if (m_DumpStats)
 				DumpRunStatistics(__OPENEMS_RUN_STAT_FILE__, t_run, currTS, speed, currE);
+			FDTD_Eng->NextInterval(speed);
 		}
 	}
 	if ((change>endCrit) && (FDTD_Op->GetExcitationSignal()->GetExciteType()==0))
@@ -1220,6 +1265,8 @@ void openEMS::RunFDTD()
 
 	//*************** postproc ************//
 	PA->PostProcess();
+
+	Signal::SetupHandlerForSIGINT(SIGNAL_ORIGINAL);
 }
 
 bool openEMS::DumpStatistics(const string& filename, double time)

openems.h

@@ -76,7 +76,7 @@ public:
 	void SetTimeStepFactor(double val) {m_TS_fac=val;}
 	void SetMaxTime(double val) {m_maxTime=val;}
 
-	void SetNumberOfThreads(unsigned int val) {m_engine_numThreads = val;}
+	void SetNumberOfThreads(int val);
 
 	void DebugMaterial() {DebugMat=true;}
 	void DebugOperator() {DebugOp=true;}
@@ -105,6 +105,7 @@ public:
 	void SetSinusExcite(double f0);
 	void SetDiracExcite(double f_max);
 	void SetStepExcite(double f_max);
+	void SetCustomExcite(std::string str, double f0, double fmax);
 
 	Excitation* InitExcitation();
 

python/Tutorials/Bent_Patch_Antenna.py

@@ -3,10 +3,10 @@
  Bent Patch Antenna Tutorial
 
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.33+
+  - openEMS v0.0.35+
 
- (C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
+ (c) 2016-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 """
 
@@ -131,11 +131,12 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 
 if not post_proc_only:
-    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+    FDTD.Run(Sim_Path, cleanup=True)
 
 ### Postprocessing & plotting
 f = np.linspace(max(1e9,f0-fc),f0+fc,401)
@@ -174,7 +175,7 @@ else:
 
     figure(figsize=(15, 7))
     ax = subplot(121, polar=True)
-    E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + nf2ff_res_phi0.Dmax
+    E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + 10.0*np.log10(nf2ff_res_phi0.Dmax)
     ax.plot(np.deg2rad(theta), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xz-plane')
     ax.grid(True)
     ax.set_xlabel('theta (deg)')
@@ -186,11 +187,11 @@ else:
     nf2ff_res_theta90 = nf2ff.CalcNF2FF(Sim_Path, f_res, 90, phi, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xy.h5')
 
     ax = subplot(122, polar=True)
-    E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + nf2ff_res_theta90.Dmax
+    E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + 10.0*np.log10(nf2ff_res_theta90.Dmax)
     ax.plot(np.deg2rad(phi), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xy-plane')
     ax.grid(True)
     ax.set_xlabel('phi (deg)')
-    suptitle('Bent Patch Anteanna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14)
+    suptitle('Bent Patch Antenna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14)
     ax.legend(loc=3)
 
     print( 'radiated power: Prad = {:.2e} Watt'.format(nf2ff_res_theta90.Prad[0]))

python/Tutorials/CRLH_Extraction.py

@@ -2,14 +2,12 @@
 """
  Tutorials / CRLH_Extraction
 
- Description at:
- http://openems.de/index.php/Tutorial:_CRLH_Extraction
-
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.34+
+  - openEMS v0.0.35+
+
+ (c) 2016-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
- (C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
 """
 
 
@@ -20,6 +18,7 @@ from pylab import *
 from CSXCAD  import ContinuousStructure
 from openEMS import openEMS
 from openEMS.physical_constants import *
+from openEMS.automesh import mesh_hint_from_box
 
 ### Class to represent single CRLH unit cells
 class CRLH_Cells:
@@ -44,36 +43,29 @@ class CRLH_Cells:
         self.edge_resolution = res
 
     def createCell(self, translate = [0,0,0]):
-        def append_mesh(mesh1, mesh2):
+        mesh = [None,None,None]
-            for n in range(3):
+        third_res = self.edge_resolution/3
-                if mesh1[n] is None:
-                    mesh1[n] = mesh2[n]
-                elif mesh2[n] is None:
-                    continue
-                else:
-                    mesh1[n] += mesh2[n]
-            return mesh1
         translate = array(translate)
         start = [-self.LL/2 , -self.LW/2, self.Top] + translate
         stop  = [-self.GLT/2,  self.LW/2, self.Top] + translate
         box = self.props['metal_top'].AddBox(start, stop, priority=10)
-        mesh = box.GetGridHint('x', metal_edge_res=self.edge_resolution, down_dir=False)
+        mesh = mesh_hint_from_box(box, 'x', metal_edge_res=self.edge_resolution, down_dir=False, mesh=mesh)
-        append_mesh(mesh, box.GetGridHint('y', metal_edge_res=self.edge_resolution) )
+        mesh = mesh_hint_from_box(box, 'y', metal_edge_res=self.edge_resolution, mesh=mesh)
 
         start = [+self.LL/2 , -self.LW/2, self.Top] + translate
         stop  = [+self.GLT/2,  self.LW/2, self.Top] + translate
         box = self.props['metal_top'].AddBox(start, stop, priority=10)
-        append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution, up_dir=False) )
+        mesh = mesh_hint_from_box(box, 'x', metal_edge_res=self.edge_resolution, up_dir=False, mesh=mesh)
 
         start = [-(self.LL-self.GLB)/2, -self.LW/2, self.Bot] + translate
         stop  = [+(self.LL-self.GLB)/2,  self.LW/2, self.Bot] + translate
         box = self.props['metal_bot'].AddBox(start, stop, priority=10)
-        append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution) )
+        mesh = mesh_hint_from_box(box, 'x', metal_edge_res=self.edge_resolution, mesh=mesh)
 
         start = [-self.SW/2, -self.LW/2-self.SL, self.Bot] + translate
         stop  = [+self.SW/2,  self.LW/2+self.SL, self.Bot] + translate
         box = self.props['metal_bot'].AddBox(start, stop, priority=10)
-        append_mesh(mesh, box.GetGridHint('xy', metal_edge_res=self.edge_resolution) )
+        mesh = mesh_hint_from_box(box, 'xy', metal_edge_res=self.edge_resolution, mesh=mesh)
 
         start = [0, -self.LW/2-self.SL+self.SW/2, 0       ] + translate
         stop  = [0, -self.LW/2-self.SL+self.SW/2, self.Bot] + translate
@@ -170,10 +162,11 @@ if __name__ == '__main__':
         if not os.path.exists(Sim_Path):
             os.mkdir(Sim_Path)
         CSX.Write2XML(CSX_file)
-        os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+        from CSXCAD import AppCSXCAD_BIN
+        os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
     if not post_proc_only:
-        FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+        FDTD.Run(Sim_Path, cleanup=True)
 
     ### Post-Processing
     f = linspace( f_start, f_stop, 1601 )

python/Tutorials/Helical_Antenna.py

@@ -3,10 +3,10 @@
  Helical Antenna Tutorial
 
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.33+
+  - openEMS v0.0.35+
 
- (C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
+ (c) 2015-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 """
 
@@ -123,10 +123,11 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 if not post_proc_only:
-    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+    FDTD.Run(Sim_Path, cleanup=True)
 
 ### Postprocessing & plotting
 freq = linspace( f0-fc, f0+fc, 501 )

python/Tutorials/MSL_NotchFilter.py

@@ -2,14 +2,11 @@
 """
  Microstrip Notch Filter Tutorial
 
- Description at:
- http://openems.de/doc/openEMS/Tutorials.html#microstrip-notch-filter
-
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.34+
+  - openEMS v0.0.35+
 
- (C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
+ (c) 2016-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 """
 
@@ -98,11 +95,12 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 
 if not post_proc_only:
-    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+    FDTD.Run(Sim_Path, cleanup=True)
 
 ### Post-processing and plotting
 f = linspace( 1e6, f_max, 1601 )

python/Tutorials/RCS_Sphere.py

@@ -3,10 +3,10 @@
  Tutorials / radar cross section of a metal sphere
 
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.34+
+  - openEMS v0.0.35+
 
- (C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
+ (c) 2016-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 """
 
 ### Import Libraries
@@ -79,11 +79,12 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 
 if not post_proc_only:
-    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+    FDTD.Run(Sim_Path, cleanup=True)
 
 ### Postprocessing & plotting
 # get Gaussian pulse strength at frequency f0

python/Tutorials/Rect_Waveguide.py

@@ -2,14 +2,11 @@
 """
  Rectangular Waveguide Tutorial
 
- Description at:
- http://openems.de/doc/openEMS/Tutorials.html#rectangular-waveguide
-
  Tested with
-  - python 3.4
+  - python 3.10
-  - openEMS v0.0.34+
+  - openEMS v0.0.35+
 
- (C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
+ (c) 2015-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 """
 
@@ -88,10 +85,11 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 if not post_proc_only:
-    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
+    FDTD.Run(Sim_Path, cleanup=True)
 
 ### Postprocessing & plotting
 freq = linspace(f_start,f_stop,201)

python/Tutorials/Simple_Patch_Antenna.py

@@ -1,8 +1,13 @@
 # -*- coding: utf-8 -*-
 """
-Created on Fri Dec 18 20:56:53 2015
+ Simple Patch Antenna Tutorial
+
+ Tested with
+  - python 3.10
+  - openEMS v0.0.34+
+
+ (c) 2015-2023 Thorsten Liebig <thorsten.liebig@gmx.de>
 
-@author: thorsten
 """
 
 ### Import Libraries
@@ -102,7 +107,8 @@ if 0:  # debugging only
     if not os.path.exists(Sim_Path):
         os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
-    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
+    from CSXCAD import AppCSXCAD_BIN
+    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))
 
 if not post_proc_only:
     FDTD.Run(Sim_Path, verbose=3, cleanup=True)
@@ -130,7 +136,7 @@ else:
     nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3])
 
     figure()
-    E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + nf2ff_res.Dmax[0]
+    E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10.0*np.log10(nf2ff_res.Dmax[0])
     plot(theta, np.squeeze(E_norm[:,0]), 'k-', linewidth=2, label='xz-plane')
     plot(theta, np.squeeze(E_norm[:,1]), 'r--', linewidth=2, label='yz-plane')
     grid()

python/openEMS/_nf2ff.pxd

@@ -20,7 +20,6 @@ from libcpp.string cimport string
 from libcpp.vector cimport vector
 from libcpp.complex cimport complex
 from libcpp cimport bool
-cimport cython.numeric
 
 cdef extern from "openEMS/nf2ff.h":
     cdef cppclass cpp_nf2ff "nf2ff":

python/openEMS/automesh.py

@@ -10,11 +10,25 @@ import numpy as np
 
 from CSXCAD import CSPrimitives
 from CSXCAD.Utilities import CheckNyDir, GetMultiDirs
+from openEMS.physical_constants import C0
+
+def mesh_combine(mesh1, mesh2, sort=True):
+    mesh = [None, None, None]
+    for ny in range(3):
+        if mesh1[ny] is None and mesh2[ny] is None:
+            continue
+        elif mesh1[ny] is None:
+            mesh[ny] = mesh2[ny]
+        elif mesh2[ny] is None:
+            mesh[ny] = mesh1[ny]
+        else:
+            mesh[ny] = list(sorted(mesh1[ny] + mesh2[ny]))
+    return mesh
 
 def mesh_hint_from_primitive(primitive, dirs, **kw):
-    if primitive.GetType() is CSPrimitives.POINT:
+    if primitive.GetType() == CSPrimitives.POINT:
         return mesh_hint_from_point(primitive, dirs, **kw)
-    if primitive.GetType() is CSPrimitives.BOX:
+    if primitive.GetType() == CSPrimitives.BOX:
         return mesh_hint_from_box(primitive, dirs, **kw)
     else:
         return None
@@ -25,12 +39,15 @@ def mesh_hint_from_point(point, dirs, **kw):
     Get a grid hint for the coordinates of the point.
 
     :param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all'
+    :param mesh: combine mesh hint to existing mesh
     :returns: (3,) list of mesh hints
     """
     hint = [None, None, None]
     coord = point.GetCoord()
     for ny in GetMultiDirs(dirs):
         hint[ny] = [coord[ny],]
+    if 'mesh' in kw:
+        return mesh_combine(hint, kw['mesh'])
     return hint
 
 def mesh_hint_from_box(box, dirs, **kw):
@@ -41,6 +58,9 @@ def mesh_hint_from_box(box, dirs, **kw):
 
     :param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all'
     :param metal_edge_res: float -- 2D flat edge resolution
+    :param up_dir: bool -- Enable upper edge
+    :param down_dir: bool -- Enable lower edge
+    :param mesh: combine mesh hint to existing mesh
     :returns: (3,) list of mesh hints
     """
     metal_edge_res = kw.get('metal_edge_res', None)
@@ -73,5 +93,23 @@ def mesh_hint_from_box(box, dirs, **kw):
                 hint[ny].append(stop[ny])
         else:
             hint[ny].append(start[ny])
+    if 'mesh' in kw:
+        return mesh_combine(hint, kw['mesh'])
     return hint
 
+def mesh_estimate_cfl_timestep(mesh):
+    """ mesh_estimate_cfl_timestep(mesh)
+
+    Estimate the maximum CFL time step of the given mesh needed to ensure numerical stability,
+    assuming propagation in pure vacuum.
+
+    :returns: the maximum CFL time step, in seconds
+    """
+    invMinDiff = [None, None, None]
+    for ny in range(3):
+        invMinDiff[ny] = np.min(np.diff(mesh.GetLines(ny))) ** -2
+
+    delta_t = mesh.GetDeltaUnit() / (C0 * np.sqrt(np.sum(invMinDiff)))
+
+    return delta_t
+

python/openEMS/openEMS.pxd

@@ -51,12 +51,15 @@ cdef extern from "openEMS/openems.h":
         void SetSinusExcite(double f0)
         void SetDiracExcite(double f_max)
         void SetStepExcite(double f_max)
+        void SetCustomExcite(string _str, double f0, double fmax)
 
         void SetAbort(bool val)
 
         void SetVerboseLevel(int level)
-        void DebugPEC()      nogil
         void DebugMaterial() nogil
+        void DebugPEC()      nogil
+        void DebugOperator() nogil
+        void DebugBox()      nogil
         void DebugCSX()      nogil
 
         int SetupFDTD() nogil

python/openEMS/openEMS.pyx

@@ -259,6 +259,17 @@ cdef class openEMS:
         """
         self.thisptr.SetStepExcite(f_max)
 
+    def SetCustomExcite(self, _str, f0, fmax):
+        """ SetCustomExcite(_str, f0, fmax)
+
+        Set a custom function as excitation signal. The custom function is supplied as a string which gets
+        parsed using `Function Parser for C++ <http://warp.povusers.org/FunctionParser/fparser.html>`_.
+
+        :param _str: str -- Custom function as string literal
+        :param f0: -- Base frequency.
+        :param fmax: -- Maximum frequency.
+        """
+        self.thisptr.SetCustomExcite(_str, f0, fmax)
 
     def SetBoundaryCond(self, BC):
         """ SetBoundaryCond(BC)
@@ -270,7 +281,7 @@ cdef class openEMS:
         * 0 or 'PEC' : perfect electric conductor (default)
         * 1 or 'PMC' : perfect magnetic conductor, useful for symmetries
         * 2 or 'MUR' : simple MUR absorbing boundary conditions
-        * 3 or 'PML-8' : PML absorbing boundary conditions
+        * 3 or 'PML_8' : PML absorbing boundary conditions
 
         :param BC: (8,) array or list -- see options above
         """
@@ -447,7 +458,8 @@ cdef class openEMS:
                     continue
                 grid.AddLine(n, hint[n])
 
-    def Run(self, sim_path, cleanup=False, setup_only=False, debug_pec=False, verbose=None, **kw):
+    def Run(self, sim_path, cleanup=False,setup_only=False, debug_material=False, debug_pec=False,
+            debug_operator=False, debug_boxes=False, debug_csx=False, verbose=None, **kw):
         """ Run(sim_path, cleanup=False, setup_only=False, verbose=None)
 
         Run the openEMS FDTD simulation.
@@ -461,19 +473,31 @@ cdef class openEMS:
         :param numThreads: int -- set the number of threads (default 0 --> max)
         """
         if cleanup and os.path.exists(sim_path):
-            shutil.rmtree(sim_path)
+            shutil.rmtree(sim_path, ignore_errors=True)
             os.mkdir(sim_path)
         if not os.path.exists(sim_path):
             os.mkdir(sim_path)
         os.chdir(sim_path)
         if verbose is not None:
             self.thisptr.SetVerboseLevel(verbose)
+        if debug_material:
+            with nogil:
+                self.thisptr.DebugMaterial()
         if debug_pec:
             with nogil:
                 self.thisptr.DebugPEC()
+        if debug_operator:
+            with nogil:
+                self.thisptr.DebugOperator()
+        if debug_boxes:
+            with nogil:
+                self.thisptr.DebugBox()
+        if debug_csx:
+            with nogil:
+                self.thisptr.DebugCSX()
         if 'numThreads' in kw:
             self.thisptr.SetNumberOfThreads(int(kw['numThreads']))
-        assert os.getcwd() == sim_path
+        assert os.getcwd() == os.path.realpath(sim_path)
         _openEMS.WelcomeScreen()
         cdef int EC
         with nogil:

python/openEMS/ports.py

@@ -64,8 +64,8 @@ class Port(object):
         self.CSX      = CSX
         self.number   = port_nr
         self.excite   = excite
-        self.start    = np.array(start, np.float)
+        self.start    = np.array(start, np.double)
-        self.stop     = np.array(stop, np.float)
+        self.stop     = np.array(stop, np.double)
         self.Z_ref    = None
         self.U_filenames = kw.get('U_filenames', [])
         self.I_filenames = kw.get('I_filenames', [])
@@ -126,7 +126,7 @@ class Port(object):
         self.uf_ref = self.uf_tot - self.uf_inc
         self.if_ref = self.if_inc - self.if_tot
 
-        if type(self.Z_ref) == float:
+        if type(self.Z_ref) in [int, float]:
             self.ut_inc = 0.5 * ( self.ut_tot + self.it_tot * self.Z_ref )
             self.it_inc = 0.5 * ( self.it_tot + self.ut_tot / self.Z_ref )
             self.ut_ref = self.ut_tot - self.ut_inc
@@ -243,7 +243,7 @@ class MSLPort(Port):
         if meas_pos_idx>=len(prop_lines)-1:
             meas_pos_idx=len(prop_lines)-2
         self.measplane_shift = np.abs(self.start[self.prop_ny]-prop_lines[meas_pos_idx])
-        prope_idx = np.array([meas_pos_idx-1, meas_pos_idx, meas_pos_idx+1], np.int)
+        prope_idx = np.array([meas_pos_idx-1, meas_pos_idx, meas_pos_idx+1], int)
         if self.direction<0:
             prope_idx = np.flipud(prope_idx)
         u_prope_pos = prop_lines[prope_idx]
@@ -300,9 +300,11 @@ class MSLPort(Port):
     def ReadUIData(self, sim_path, freq, signal_type ='pulse'):
         self.u_data = UI_data(self.U_filenames, sim_path, freq, signal_type )
         self.uf_tot = self.u_data.ui_f_val[1]
+        self.ut_tot = self.u_data.ui_val[1]
 
         self.i_data = UI_data(self.I_filenames, sim_path, freq, signal_type )
         self.if_tot = 0.5*(self.i_data.ui_f_val[0]+self.i_data.ui_f_val[1])
+        self.it_tot = 0.5*(self.i_data.ui_val[0]+self.i_data.ui_val[1])
 
         unit = self.CSX.GetGrid().GetDeltaUnit()
         Et = self.u_data.ui_f_val[1]

python/setup.py

@@ -22,7 +22,7 @@ extensions = [
 
 setup(
   name="openEMS",
-  version = '0.0.33',
+  version = '0.0.36',
   description = "Python interface for the openEMS FDTD library",
   classifiers = [
       'Development Status :: 3 - Alpha',
@@ -31,6 +31,12 @@ setup(
       'Intended Audience :: Science/Research',
       'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
       'Programming Language :: Python',
+      'Programming Language :: Python :: 3',
+      'Programming Language :: Python :: 3.9',
+      'Programming Language :: Python :: 3.10',
+      'Programming Language :: Python :: 3.11',
+      'Programming Language :: Python :: 3.12',
+      'Programming Language :: Python :: Implementation :: CPython',
       'Topic :: Scientific/Engineering',
       'Topic :: Software Development :: Libraries :: Python Modules',
       'Operating System :: POSIX :: Linux',
@@ -43,5 +49,11 @@ setup(
   url = 'https://openEMS.de',
   packages=["openEMS", ],
   package_data={'openEMS': ['*.pxd']},
-  ext_modules = cythonize(extensions, language_level = "3")
+  python_requires='>=3.9',
+  install_requires=[
+    'cython>=3.0.6',  # Apache License 2.0 (https://github.com/cython/cython/blob/master/LICENSE.txt)
+    'h5py>=3.10.0',   # BSD 3-Clause (https://github.com/h5py/h5py/blob/master/LICENSE)
+    'numpy>=1.26.2',  # BSD 3-Clause (https://github.com/numpy/numpy/blob/main/LICENSE.txt)
+  ],
+  ext_modules = cythonize(extensions, language_level = 3)
  )

tools/CMakeLists.txt

@@ -4,6 +4,7 @@ set(SOURCES
   ${CMAKE_CURRENT_SOURCE_DIR}/AdrOp.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/ErrorMsg.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/array_ops.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/signal.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/global.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/hdf5_file_reader.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/hdf5_file_writer.cpp

tools/signal.cpp

@@ -0,0 +1,263 @@
+/*
+*	Copyright (C) 2023 Yifeng Li <tomli@tomli.me>
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cstdio>
+#include <cstring>
+#include <cstdlib>
+#include "signal.h"
+
+void Signal::SetupHandlerForSIGINT(int type)
+{
+	m_sigintAbort = 0;
+
+#ifndef WIN32
+	UnixSetupHandlerForSIGINT(type);
+#else
+	Win32SetupHandlerForConsoleCtrl(type);
+#endif
+}
+
+#ifndef WIN32
+void Signal::UnixSetupHandlerForSIGINT(int type)
+{
+	if (type == SIGNAL_ORIGINAL && m_sigHandlerOriginal)
+	{
+		// If we're acting as a shared library and a program (such as
+		// the Python interpreter) calls us, changing the SIGINT handler
+		// unilaterally may overwrite the original handler and affect
+		// the functionality of the original program. Thus, we save the
+		// original handler and restore it after the end of SetupFDTD()
+		// or RunFDTD() to minimize the disruption.
+		auto retval = std::signal(SIGINT, m_sigHandlerOriginal);
+		if (retval == SIG_ERR)
+		{
+			fprintf(stderr, "Signal::UnixSetupHandlerForSIGINT(): "
+					"Failed to restore signal handler!\n");
+		}
+		m_sigHandlerOriginal = NULL;
+	}
+	else if (type == SIGNAL_EXIT_GRACEFUL)
+	{
+		m_sigHandlerOriginal = std::signal(SIGINT, UnixGracefulExitHandler);
+		if (m_sigHandlerOriginal == SIG_ERR)
+		{
+			fprintf(stderr, "Signal::UnixSetupHandlerForSIGINT(): "
+					"Failed to set UnixGracefulExitHandler!\n");
+			m_sigHandlerOriginal = NULL;
+		}
+	}
+	else if (type == SIGNAL_EXIT_FORCE)
+	{
+		m_sigHandlerOriginal = std::signal(SIGINT, UnixForceExitHandler);
+		if (m_sigHandlerOriginal == SIG_ERR)
+		{
+			fprintf(stderr, "Signal::UnixSetupHandlerForSIGINT(): "
+					"Failed to set UnixForceExitHandler!\n");
+			m_sigHandlerOriginal = NULL;
+		}
+	}
+}
+
+void Signal::UnixGracefulExitHandler(int signal)
+{
+	m_sigintAbort = 1;
+
+	// C standard only guarantees that a sig_atomic_t variable is safe
+	// to read or write, but it's not necessarily safe to increment by
+	// one, and also not safe to set one sig_atomic_t depending on the
+	// result of another sig_atomic_t.
+	//
+	// Thus, we switch the signal handler itself instead of recording
+	// the number of times SIGINT is raised.
+	auto retval = std::signal(SIGINT, UnixForceExitHandler);
+	if (retval == SIG_ERR)
+	{
+		SafeStderrWrite("\nSignal::UnixGracefulExitHandler(): "
+				"Failed to set UnixForceExitHandler!");
+	}
+	else
+	{
+		SafeStderrWrite("\nSignal::UnixGracefulExitHandler(): "
+				"Gracefully aborting simulation "
+				"now, this may take a few seconds...\n"
+				"Signal::UnixGracefulExitHandler(): "
+				"To force-exit, send Ctrl-C again, "
+				"but simulation results may be lost.\n");
+	}
+}
+
+void Signal::UnixForceExitHandler(int signal)
+{
+	SafeStderrWrite("\nSignal::UnixForceExitHandler(): "
+			"Force-exit simulation process now!\n");
+
+	// By convention, if a program is (uncleanly) aborted due to
+	// an external signal, preferably it should return 128 + signal.
+	// For SIGINT, it's 130.
+	std::_Exit(128 + signal);
+}
+
+#else
+
+void Signal::Win32SetupHandlerForConsoleCtrl(int type)
+{
+	if (type == SIGNAL_ORIGINAL || m_sigHandlerRegistered)
+	{
+		// On Windows, SetConsoleCtrlHandler appends a new ConsoleCtrlHandler
+		// in addition to the existing handlers. Thus, we need to record
+		// the ConsoleCtrlHandler installed by us (instead of getting the
+		// pre-existing handlers on Unix). Then, before we install a new
+		// signal handler, we need to use the argument "Add == FALSE" to
+		// remove the handler we previously installed.
+		//
+		// We also need to do the same in case that we're restoring the
+		// ConsoleCtrlHandler to the original state (note how on Unix, the
+		// if expression uses "AND", but on Windows, the if expression uses
+		// "OR".
+		BOOL success = SetConsoleCtrlHandler(m_sigHandlerRegistered, FALSE);
+		m_sigHandlerRegistered = NULL;
+
+		if (!success)
+		{
+			fprintf(stderr, "Signal::Win32SetupHandlerForConsoleCtrl(): "
+					"Failed to unregister ConsoleCtrlHandler!\n");
+			return;
+		}
+	}
+
+	// Assume m_sigHandlerRegistered has already been unregistered.
+	if (type == SIGNAL_EXIT_GRACEFUL)
+	{
+		m_sigHandlerRegistered = (PHANDLER_ROUTINE) Win32GracefulExitHandler;
+		BOOL success = SetConsoleCtrlHandler(m_sigHandlerRegistered, TRUE);
+
+		if (!success)
+		{
+			fprintf(stderr, "Signal::Win32SetupHandlerForConsoleCtrl(): "
+					"Failed to register Win32GracefulExitHandler!\n");
+		}
+	}
+	else if (type == SIGNAL_EXIT_FORCE)
+	{
+		m_sigHandlerRegistered = (PHANDLER_ROUTINE) Win32ForceExitHandler;
+		BOOL success = SetConsoleCtrlHandler(m_sigHandlerRegistered, TRUE);
+
+		if (!success)
+		{
+			fprintf(stderr, "Signal::Win32SetupHandlerForConsoleCtrl(): "
+					"Failed to register Win32ForceExitHandler!\n");
+		}
+	}
+}
+
+BOOL Signal::Win32GracefulExitHandler(DWORD fdwCtrlType)
+{
+	m_sigintAbort = 1;
+
+	// unregister the current handler
+	BOOL success = SetConsoleCtrlHandler(m_sigHandlerRegistered, FALSE);
+	if (!success)
+	{
+		SafeStderrWrite("Signal::Win32GracefulExitHandler(): "
+				"Failed to unregister Win32GracefulExitHandler!\n");
+		return true;
+	}
+
+	// install a new handler
+	m_sigHandlerRegistered = (PHANDLER_ROUTINE) Win32ForceExitHandler;
+	success = SetConsoleCtrlHandler(m_sigHandlerRegistered, TRUE);
+	if (!success)
+	{
+		SafeStderrWrite("Signal::Win32GracefulExitHandler(): "
+				"Failed to register Win32ForceExitHandler!\n");
+	}
+	else
+	{
+		SafeStderrWrite("\nSignal::Win32GracefulExitHandler(): "
+				"Gracefully aborting simulation "
+				"now, this may take a few seconds...\n"
+				"Signal::Win32GracefulExitHandler(): "
+				"To force-exit, send Ctrl-C again, "
+				"but simulation results may be lost.\n");
+	}
+
+	return true;
+}
+
+BOOL Signal::Win32ForceExitHandler(DWORD fdwCtrlType)
+{
+	SafeStderrWrite("\nSignal::Win32ForceExitHandler(): "
+			"Force-exit simulation process now!\n");
+
+	// On Windows, the exit code for SIGINT is always 3.
+	std::_Exit(3);
+
+	// unreachable
+	return true;
+}
+#endif
+
+bool Signal::ReceivedSIGINT(void)
+{
+	if (m_sigintAbort)
+		return true;
+	else
+		return false;
+}
+
+void Signal::SafeStderrWrite(const char *buf)
+{
+#ifdef WIN32
+	// On Windows, using any kind of system calls in a ANSI C signal
+	// handler is prohibited, in this case, this function should return
+	// immediately without doing anything. But, when the official way
+	// SetConsoleCtrlHandler() is used (instead of using ANSI C signals),
+	// there's no such restriction.
+	fprintf(stderr, "%s", buf);
+	fflush(stderr);
+	return;
+#else
+	// On Unix, in a signal handler, it's unsafe to use normal I/O
+	// functions such as iostream, puts(), printf(), fprintf(). The
+	// only safe option is the system call write().
+	size_t buf_len = strlen(buf);
+	ssize_t bytes = 0;
+
+	while (buf_len > 0)
+	{
+		bytes = write(STDERR_FILENO, buf, buf_len);
+		if (bytes < 0)
+		{
+			// write failure, nothing we can do.
+			return;
+		}
+
+		if ((size_t) bytes > buf_len)
+		{
+			// Assertion: This should never happen. bytes is
+			// always less or equal to buf_len, and buf_len
+			// will never underflow under any circumstances
+			// (unless the write system call is broken).
+			return;
+		}
+
+		buf += bytes;		    // advance buffer position
+		buf_len -= (size_t) bytes;  // decrement limiter
+	}
+	return;  // write completed.
+#endif
+}

tools/signal.h

@@ -0,0 +1,59 @@
+/*
+*	Copyright (C) 2023 Yifeng Li <tomli@tomli.me>
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef SIGNAL_H
+#define SIGNAL_H
+
+#include <csignal>
+
+#ifndef WIN32
+#include <unistd.h>
+#else
+#include <windows.h>
+#endif
+
+enum
+{
+	SIGNAL_ORIGINAL,
+	SIGNAL_EXIT_GRACEFUL,
+	SIGNAL_EXIT_FORCE,
+};
+
+class Signal
+{
+public:
+	static void SetupHandlerForSIGINT(int type);
+	static bool ReceivedSIGINT(void);
+
+private:
+	inline static volatile std::sig_atomic_t m_sigintAbort = 0;
+	static void SafeStderrWrite(const char *buf);
+
+#ifndef WIN32
+	inline static void (*m_sigHandlerOriginal)(int) = NULL;
+	static void UnixSetupHandlerForSIGINT(int type);
+	static void UnixGracefulExitHandler(int signal);
+	static void UnixForceExitHandler(int signal);
+#else
+	inline static PHANDLER_ROUTINE m_sigHandlerRegistered = NULL;
+	static void Win32SetupHandlerForConsoleCtrl(int type);
+	static BOOL Win32GracefulExitHandler(DWORD fdwCtrlType);
+	static BOOL Win32ForceExitHandler(DWORD fdwCtrlType);
+#endif
+};
+
+#endif