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>
set C++ version for boost thread
2024-09-16 19:49:23 +02:00 · 2024-03-11 21:48:06 +01:00 · 2023-12-16 15:48:58 +01:00 · 2023-12-16 15:41:31 +01:00 · 2023-12-16 11:05:02 +01:00 · 2023-12-16 11:01:36 +01:00
40 changed files with 1824 additions and 204 deletions
--- a/CMakeLists.txt
+++ b/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} )
-else()
-    set( vtk_LIBS
-        vtkCommon
-        vtkIO
-    )    
-endif()
+set( vtk_LIBS ${VTK_LIBRARIES} )
 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()
--- a/FDTD/engine.cpp
+++ b/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]] += 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]]);
+				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]] +=
+				    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]] += 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]]);
+				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]] +=
+				    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]] += 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]]);
+				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]] +=
+				    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]] += 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]);
+				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]] +=
+				    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]] += 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]]);
+				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]] +=
+				    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]] += 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]]);
+				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]] +=
+				    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];
--- a/FDTD/engine.h
+++ b/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]]; }
--- a/FDTD/engine_multithread.cpp
+++ b/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)
-		cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
+	if (m_IterateBarrier!=0)
+		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)
--- a/FDTD/engine_multithread.h
+++ b/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)
--- a/FDTD/engine_sse.cpp
+++ b/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 += 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 );
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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 );
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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];
 	}
--- a/FDTD/engine_sse_compressed.cpp
+++ b/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 += 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 );
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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 );
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+				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 +=
+				    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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 += 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);
+			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 +=
+			    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];
 	}
--- a/FDTD/extensions/CMakeLists.txt
+++ b/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
 )

--- a/FDTD/extensions/engine_ext_lumpedRLC.cpp
+++ b/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;
+}
+
+
--- a/FDTD/extensions/engine_ext_lumpedRLC.h
+++ b/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
--- a/FDTD/extensions/operator_ext_lumpedRLC.cpp
+++ b/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;
+}
+
--- a/FDTD/extensions/operator_ext_lumpedRLC.h
+++ b/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_ */
--- a/FDTD/extensions/operator_ext_upml.cpp
+++ b/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 )
--- a/FDTD/operator.cpp
+++ b/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)
--- a/FDTD/operator.h
+++ b/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_ConductingSheet; //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_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];

--- a/FDTD/operator_cylinder.cpp
+++ b/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;
 }

--- a/FDTD/operator_cylindermultigrid.cpp
+++ b/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;
 }

--- a/FDTD/operator_multithread.cpp
+++ b/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;
--- a/FDTD/operator_multithread.h
+++ b/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.
 	/*!
--- a/matlab/setup.m
+++ b/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');
-        [res, fn_h]  = unix('find /usr/include -name hdf5.h | sort -r | head -1');
+        dylib_extension = 'so';
+        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
--- a/nf2ff/CMakeLists.txt
+++ b/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)
--- a/nf2ff/main.cpp
+++ b/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)
--- a/openems.cpp
+++ b/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 << " | (C) 2010-2018 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license"   << endl;
+	cout << " | openEMS " << bits << " -- version " << GIT_VERSION << 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)
--- a/openems.h
+++ b/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();

--- a/python/Tutorials/Bent_Patch_Antenna.py
+++ b/python/Tutorials/Bent_Patch_Antenna.py
@ -3,10 +3,10 @@
 Bent Patch Antenna Tutorial

 Tested with
-  - python 3.4
-  - openEMS v0.0.33+
+  - python 3.10
+  - 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]))
--- a/python/Tutorials/CRLH_Extraction.py
+++ b/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
-  - openEMS v0.0.34+
+  - python 3.10
+  - 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):
-            for n in range(3):
-                if mesh1[n] is None:
-                    mesh1[n] = mesh2[n]
-                elif mesh2[n] is None:
-                    continue
-                else:
-                    mesh1[n] += mesh2[n]
-            return mesh1
+        mesh = [None,None,None]
+        third_res = self.edge_resolution/3
        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)
-        append_mesh(mesh, box.GetGridHint('y', metal_edge_res=self.edge_resolution) )
+        mesh = mesh_hint_from_box(box, 'x', metal_edge_res=self.edge_resolution, down_dir=False, mesh=mesh)
+        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 )
--- a/python/Tutorials/Helical_Antenna.py
+++ b/python/Tutorials/Helical_Antenna.py
@ -3,10 +3,10 @@
 Helical Antenna Tutorial

 Tested with
-  - python 3.4
-  - openEMS v0.0.33+
+  - python 3.10
+  - 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 )
--- a/python/Tutorials/MSL_NotchFilter.py
+++ b/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
-  - openEMS v0.0.34+
+  - python 3.10
+  - 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 )
--- a/python/Tutorials/RCS_Sphere.py
+++ b/python/Tutorials/RCS_Sphere.py
@ -3,10 +3,10 @@
 Tutorials / radar cross section of a metal sphere

 Tested with
-  - python 3.4
-  - openEMS v0.0.34+
+  - python 3.10
+  - 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
--- a/python/Tutorials/Rect_Waveguide.py
+++ b/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
-  - openEMS v0.0.34+
+  - python 3.10
+  - 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)
--- a/python/Tutorials/Simple_Patch_Antenna.py
+++ b/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()
--- a/python/openEMS/_nf2ff.pxd
+++ b/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":
--- a/python/openEMS/automesh.py
+++ b/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
+
--- a/python/openEMS/openEMS.pxd
+++ b/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
--- a/python/openEMS/openEMS.pyx
+++ b/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:
--- a/python/openEMS/ports.py
+++ b/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.stop     = np.array(stop, np.float)
+        self.start    = np.array(start, np.double)
+        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]
--- a/python/setup.py
+++ b/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)
 )
--- a/tools/CMakeLists.txt
+++ b/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
--- a/tools/signal.cpp
+++ b/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
+}
--- a/tools/signal.h
+++ b/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
Author	SHA1	Message	Date
Yifeng Li	7d7688a288	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>	2024-09-16 19:49:23 +02:00
aWZHY0yQH81uOYvH	1ccf094247	set C++ version for boost thread	2024-03-11 21:48:06 +01:00
Thorsten Liebig	6a13d81cf0	python: add missing definition for custom excite, #129 Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-12-16 15:48:58 +01:00
Tobias Oberstein	45c90c24e8	expose SetCustomExcite() in Python (ported from PR #129 )	2023-12-16 15:41:31 +01:00
Thorsten Liebig	e5db9de99b	python: fix language level setup, #130 Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-12-16 11:05:02 +01:00
Tobias Oberstein	a87a75efc3	remove import that is actually unneeded	2023-12-16 11:01:36 +01:00
Tobias Oberstein	e7620dcd72	also add cython (everything but setuptools)	2023-12-16 11:01:36 +01:00
Tobias Oberstein	5bc31d0255	add missing package python dependencies	2023-12-16 11:01:36 +01:00
Tobias Oberstein	eb4845975a	add tested python impls to pkg spec; fixes #127	2023-12-16 11:01:36 +01:00
Tobias Oberstein	831bd7c835	re-sync python package version with current openEMS release	2023-12-16 11:01:36 +01:00
Tobias Oberstein	bf30b7d94a	fix cython import	2023-12-16 11:01:36 +01:00
Yifeng Li	71dde7ea49	FDTD: reformat code of update equations. The original update equations in the FDTD engine have extremely long lines and are difficult to read and work with. This patch inserts line breaks, it aligns all array indexes by x/y/z coordinates to make it easy to visually compare. Signed-off-by: Yifeng Li <tomli@tomli.me>	2023-11-18 17:55:46 +01:00
Chris Madsen	5b8cf2f2ed	Dmax is in linear power units, so convert to dB power before adding to dB power.	2023-11-18 13:25:28 +01:00
gammaxy1	6212d1de68	PML_8 typo in docstring	2023-11-18 13:25:28 +01:00
Yifeng Li	840c9755d5	Handle SIGINT for openEMS and Python, with graceful exit support. Currently, openEMS doesn't have any special code to handle SIGINT (which is raised by pressing Control-C). By default, the program is terminated without saving data. This worked okay in the past, but now its limitations are becoming obvious. 1. When openEMS is used as a Python module, Control-C stops working because SIGINT is now managed by Python in order to generate KeyboardInterrupt exceptions, normally this isn't a problem, but if we are running an external C++ (Cython) function such as openEMS, the Python interpreter mainloop has no control until we return. As a result, SIGINT is received but never handled. In Cython, programs are expected to call PyErr_CheckSignals() in its blocking loop periodically to temporally transfer control back to Python to handle signals. But this introduces a dependency of Cython in the FDTD mainloop. 2. During a simulation, it's not possible to abort it gracefully by pressing Control-C, this is a limitation of openEMS itself, it's always a force exit. Currently the only supported method for graceful exit is creating a file called "ABORT" in the simulation directory. If we already need to implement a signal handler, adding a graceful exit at the same time would be a good idea. This commit installs SIGINT handlers during SetupFDTD() and RunFDTD(). 1. In RunFDTD(), if SIGINT is received once, a status flag is set, which is then checked in CheckAbortCond(), allowing a graceful exit with the same effect of an "ABORT" file. If SIGINT is received twice, openEMS force exit without saving data (just like the old default behavior). 2. In SetupFDTD(), if SIGINT is received, openEMS immediately force exit without saving data, identical to the old behavior. In a huge simulation, initializing and compressing operators may have a long time. so we want an early exit before RunFDTD(). 3. Before RunFDTD() and SetupFDTD() return, the original signal handler for SIGINT is restored. This is important since when we're acting as a shared library. When 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. For example, Python would never be able to raise KeyboardInterrupt again. Thus, we save the original handler and restore it later. Signed-off-by: Yifeng Li <tomli@tomli.me>	2023-11-18 12:32:44 +01:00
G. L	ee3f2b7d80	Lumped RLC parallel & series implementation (openEMS) (#121 )	2023-11-18 12:23:15 +01:00
Thorsten Liebig	5f36e7f3a2	version 0.0.36 Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-10-22 16:25:51 +02:00
Yifeng Li	782a7381bf	CMakeLists.txt: append instead of overwrite CMAKE_CXX_FLAGS Currently, on ARM and PPC, CMakeLists.txt uses: set(CMAKE_CXX_FLAGS "-DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS") but this overwrites the default value of CMAKE_CXX_FLAGS from CMake, including user-specified CXXFLAGS via environmental variable, making it impossible to change CXXFLAGS. This patch appends instead of overwrite CMAKE_CXX_FLAGS via: set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS") Signed-off-by: Yifeng Li <tomli@tomli.me>	2023-05-17 19:01:06 +02:00
Yifeng Li	486f3140cb	openEMS.pyx: check canonical path in assert, close #113 . Currently, running openEMS's example Python scripts on macOS always fails with the following error: $ python3 MSL_NotchFilter.py Traceback (most recent call last): File "/Users/gentoo/code/openEMS-Project/openEMS/python/Tutorials/MSL_NotchFilter.py", line 103, in <module> FDTD.Run(Sim_Path, cleanup=True) File "openEMS/openEMS.pyx", line 489, in openEMS.openEMS.openEMS.Run AssertionError This is caused by an oversight of an assertion in openEMS.pyx: os.chdir(sim_path) # ... assert os.getcwd() == sim_path The problem here is that "sim_path" is not a canonical path name, so the assertion would fail if the path we're switching into contains a symbolic link. This problem affects all operating systems, it's not limited to macOS. But on macOS, the problem is especially serious, since macOS's "/tmp" is a link to "/private/tmp" by default. Thus, it causes an AssertionError in all the included Python examples. Instead of doing "assert os.getcwd() == sim_path", we should write "assert os.getcwd() == os.path.realpath(sim_path)" to ensure that we're checking a canonical path. Signed-off-by: Yifeng Li <tomli@tomli.me>	2023-05-03 18:46:29 +02:00
Thorsten Liebig	c651cce61f	Merge remote-tracking branch 'drake/feature/estimate-cfl-timestep'	2023-03-19 11:45:14 +01:00
Thorsten Liebig	568cdbdfac	PML: try to fix pml working for a finite conductor waveguide sigma > 1000 S/m is considered a conductor (not ideal solution) Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-03-08 22:10:30 +01:00
Gonzalo José Carracedo Carballal	cf34998b01	Add missing import	2023-03-06 21:46:43 +01:00
Gonzalo José Carracedo Carballal	3eb4439959	Improve readability of mesh_estimate_cfl_timestep	2023-03-06 21:45:17 +01:00
Gonzalo José Carracedo Carballal	6440b408ac	Implement mesh_estimate_cfl_timestep	2023-03-06 20:36:50 +01:00
Thorsten Liebig	cb63ab01c4	python: fix automesh type detection Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-03-05 20:33:07 +01:00
Thorsten Liebig	704fad6dc4	python: cleanup should not crash if folder cannot be removed Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-02-22 18:57:52 +01:00
Thorsten Liebig	cbbae61c24	python: fix TD for MSL ports with set ref. impedance Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-02-22 18:57:17 +01:00
Yifeng Li	8e408307b8	python/ports.py: replace deprecated "np.int" with "int". Accroding to NumPy's development team, "for a long time, np.int has been an alias of the builtin int. This is repeatedly a cause of confusion for newcomers, and existed mainly for historic reasons." This and many other aliases have been deprecated since NumPy v1.20.0, and at this point they've been completely removed. Replace "np.int" with "int" allows ports.py to run again. Signed-off-by: Yifeng Li <tomli@tomli.me>	2023-02-22 18:41:43 +01:00
Thorsten Liebig	ecf0c160e0	nf2ff main: update year info Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-02-19 12:47:33 +01:00
Thorsten Liebig	b49bd2af80	MT engine: fix threads not cleaned up, #104 Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-02-11 10:43:28 +01:00
Thorsten Liebig	0342eefd27	python tutorials: use new/better automesh options for CRLH examples Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-07 20:52:36 +01:00
Thorsten Liebig	595c8effbd	python: improve automesh options Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-07 20:50:32 +01:00
Thorsten Liebig	a0e45f8869	python: fix numpy datatype Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-07 20:50:10 +01:00
Thorsten Liebig	55068629b0	cmake: drop support for vtk<=5 Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-06 20:05:13 +01:00
Thorsten Liebig	6673aefd70	engine: try to find optimal number of engine threads Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-06 20:01:07 +01:00
Thorsten Liebig	63c5fe561d	fix hdf5 search to not find opencv hdf5.h Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>	2023-01-06 16:43:21 +01:00
aWZHY0yQH81uOYvH	115aeb64e2	update setup.m to work on Mac	2023-01-06 03:31:01 -08:00
aWZHY0yQH81uOYvH	3162c487d9	detect `arm64` under macOS	2023-01-06 00:38:53 -08:00
Apostolos	638d875906	Expose Debugs to Python API	2023-01-03 21:30:20 +01:00
Apostolos	3772497901	Expose DebugMaterial, DebugOperator and DebugBox	2023-01-03 21:30:20 +01:00
Thorsten Liebig	9677c457e8	python: update Tutorials meta data, remove verbose call Signed-off-by: Thorsten Liebig <liebig@imst.de>	2023-01-02 12:32:28 +01:00
Thorsten Liebig	0777302f1f	python Tutorial: fix CRLH mesh hints Signed-off-by: Thorsten Liebig <liebig@imst.de>	2023-01-01 14:14:15 +01:00
Thorsten Liebig	164d3983e3	python: allow windows to find AppCSXCAD Signed-off-by: Thorsten Liebig <liebig@imst.de>	2023-01-01 14:13:53 +01:00
Thorsten Liebig	df7c58d961	info: update welcome screen Signed-off-by: Thorsten Liebig <liebig@imst.de>	2022-12-30 17:18:53 +01:00