diff --git a/Common/operator_base.h b/Common/operator_base.h
index 0473ff6..063c863 100644
--- a/Common/operator_base.h
+++ b/Common/operator_base.h
@@ -41,8 +41,8 @@ public:
 	//! Get the number of timesteps satisfying the nyquist condition (may depend on the excitation)
 	virtual unsigned int GetNumberOfNyquistTimesteps() const =0;
 
-	//! Returns the number of lines as needed for post-processing etc. (for the engine, use GetOriginalNumLines())
-	virtual unsigned int GetNumberOfLines(int ny) const =0;
+	//! Returns the number of lines as needed for post-processing etc.
+	virtual unsigned int GetNumberOfLines(int ny, bool full=false) const =0;
 
 	//! Get the name for the given direction: 0 -> x, 1 -> y, 2 -> z
 	virtual std::string GetDirName(int ny) const;
@@ -76,7 +76,7 @@ public:
 	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0;
 
 	//! Snap the given coodinates to mesh indices, return box dimension
-	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const =0;
+	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const =0;
 
 	//! Snap a given box to the operator mesh, uiStart will be always <= uiStop
 	/*!
@@ -88,7 +88,7 @@ public:
 	  \param[in] SnapMethod Snapping method, 0=snap to closest line, 1/(2)=snap such that given box is inside (outside) the snapped lines
 	  \return returns the box dimension or -1 if box is not inside the simulation domain
 	  */
-	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const =0;
+	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const =0;
 
 	//! Set the boundary conditions
 	virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];}
diff --git a/Common/processing.cpp b/Common/processing.cpp
index f25720b..7cab6e6 100644
--- a/Common/processing.cpp
+++ b/Common/processing.cpp
@@ -170,7 +170,7 @@ void Processing::AddFrequency(vector<double> *freqs)
 
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	m_Dimension = Op->SnapBox2Mesh(dstart,dstop,start,stop,m_dualMesh,m_SnapMethod, m_start_inside, m_stop_inside);
+	m_Dimension = Op->SnapBox2Mesh(dstart,dstop,start,stop,m_dualMesh,false,m_SnapMethod, m_start_inside, m_stop_inside);
 	if (m_Dimension<0)
 	{
 		cerr << "Processing::DefineStartStopCoord: Warning in " << m_Name << " (" << GetProcessingName() << ") : Box is outside the field domain!! Disabling" << endl;
diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp
index 8a3aae6..26725e6 100644
--- a/FDTD/engine.cpp
+++ b/FDTD/engine.cpp
@@ -36,7 +36,7 @@ Engine::Engine(const Operator* op)
 	numTS = 0;
 	Op = op;
 	for (int n=0; n<3; ++n)
-		numLines[n] = Op->GetOriginalNumLines(n);
+		numLines[n] = Op->GetNumberOfLines(n, true);
 	volt=NULL;
 	curr=NULL;
 }
diff --git a/FDTD/engine_interface_fdtd.cpp b/FDTD/engine_interface_fdtd.cpp
index 2da5513..f59fd69 100644
--- a/FDTD/engine_interface_fdtd.cpp
+++ b/FDTD/engine_interface_fdtd.cpp
@@ -81,7 +81,7 @@ double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos,
 		{
 			nP = (n+1)%3;
 			nPP = (n+2)%3;
-			if ((pos[0]==m_Op->GetOriginalNumLines(0)-1) || (pos[1]==m_Op->GetOriginalNumLines(1)-1) || (pos[2]==m_Op->GetOriginalNumLines(2)-1))
+			if ((pos[0]==m_Op->GetNumberOfLines(0,true)-1) || (pos[1]==m_Op->GetNumberOfLines(1,true)-1) || (pos[2]==m_Op->GetNumberOfLines(2,true)-1))
 			{
 				out[n] = 0; //electric field outside the field domain is always zero
 				continue;
@@ -119,7 +119,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 		{
 			nP = (n+1)%3;
 			nPP = (n+2)%3;
-			if ((pos[0]==m_Op->GetOriginalNumLines(0)-1) || (pos[1]==m_Op->GetOriginalNumLines(1)-1) || (pos[2]==m_Op->GetOriginalNumLines(2)-1) || (pos[nP]==0) || (pos[nPP]==0))
+			if ((pos[0]==m_Op->GetNumberOfLines(0,true)-1) || (pos[1]==m_Op->GetNumberOfLines(1,true)-1) || (pos[2]==m_Op->GetNumberOfLines(2,true)-1) || (pos[nP]==0) || (pos[nPP]==0))
 			{
 				out[n] = 0;
 				continue;
@@ -140,7 +140,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 		{
 			delta = m_Op->GetEdgeLength(n,iPos,true);
 			out[n] = m_Eng->GetCurr(n,iPos);
-			if ((pos[n]>=m_Op->GetOriginalNumLines(n)-1))
+			if ((pos[n]>=m_Op->GetNumberOfLines(n,true)-1))
 			{
 				out[n] = 0; //magnetic field on the outer boundaries is always zero
 				continue;
diff --git a/FDTD/extensions/engine_ext_cylinder.cpp b/FDTD/extensions/engine_ext_cylinder.cpp
index af17c52..e240f1b 100644
--- a/FDTD/extensions/engine_ext_cylinder.cpp
+++ b/FDTD/extensions/engine_ext_cylinder.cpp
@@ -29,7 +29,7 @@ Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine
 	CC_R0_included = op_ext->CC_R0_included;
 
 	for (int n=0; n<3; ++n)
-		numLines[n] = op_ext->m_Op->GetOriginalNumLines(n);
+		numLines[n] = op_ext->m_Op->GetNumberOfLines(n,true);
 
 	//this cylindrical extension should be executed first?
 	m_Priority = ENG_EXT_PRIO_CYLINDER;
diff --git a/FDTD/extensions/operator_ext_conductingsheet.cpp b/FDTD/extensions/operator_ext_conductingsheet.cpp
index 6e1a7ae..74e9a76 100644
--- a/FDTD/extensions/operator_ext_conductingsheet.cpp
+++ b/FDTD/extensions/operator_ext_conductingsheet.cpp
@@ -44,7 +44,7 @@ bool Operator_Ext_ConductingSheet::BuildExtension()
 	double dT = m_Op->GetTimestep();
 	unsigned int pos[] = {0,0,0};
 	double coord[3];
-	unsigned int numLines[3] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)};
+	unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
 
 	m_Order = 0;
 	vector<unsigned int> v_pos[3];
diff --git a/FDTD/extensions/operator_ext_cylinder.cpp b/FDTD/extensions/operator_ext_cylinder.cpp
index dead36b..277cad2 100644
--- a/FDTD/extensions/operator_ext_cylinder.cpp
+++ b/FDTD/extensions/operator_ext_cylinder.cpp
@@ -50,19 +50,19 @@ bool Operator_Ext_Cylinder::BuildExtension()
 	if (CC_R0_included==false)
 		return true;
 
-	vv_R0 = new FDTD_FLOAT[m_Op->GetOriginalNumLines(2)];
-	vi_R0 = new FDTD_FLOAT[m_Op->GetOriginalNumLines(2)];
+	vv_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)];
+	vi_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)];
 
 	unsigned int pos[3];
 	double coord[3];
 	double inEC[4];
 	double dT = m_Op->GetTimestep();
 	pos[0]=0;
-	for (pos[2]=0; pos[2]<m_Op->GetOriginalNumLines(2); ++pos[2])
+	for (pos[2]=0; pos[2]<m_Op->GetNumberOfLines(2,true); ++pos[2])
 	{
 		double C=0;
 		double G=0;
-		for (pos[1]=0; pos[1]<m_Op->GetOriginalNumLines(1)-2; ++pos[1])
+		for (pos[1]=0; pos[1]<m_Op->GetNumberOfLines(1,true)-2; ++pos[1])
 		{
 			m_Op_Cyl->Calc_ECPos(2,pos,inEC);
 			C+=inEC[0];
@@ -72,7 +72,7 @@ bool Operator_Ext_Cylinder::BuildExtension()
 		vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);
 		vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C);
 
-		for (unsigned int i=0; i<m_Op->GetOriginalNumLines(1); ++i)
+		for (unsigned int i=0; i<m_Op->GetNumberOfLines(1,true); ++i)
 		{
 			m_Op->EC_C[2][m_Op->MainOp->SetPos(0,i,pos[2])] = C;
 			m_Op->EC_G[2][m_Op->MainOp->SetPos(0,i,pos[2])] = G;
diff --git a/FDTD/extensions/operator_ext_excitation.cpp b/FDTD/extensions/operator_ext_excitation.cpp
index 0fd7109..44a3182 100644
--- a/FDTD/extensions/operator_ext_excitation.cpp
+++ b/FDTD/extensions/operator_ext_excitation.cpp
@@ -137,7 +137,7 @@ bool Operator_Ext_Excitation::BuildExtension()
 	CSProperties* prop=NULL;
 	int priority=0;
 
-	unsigned int numLines[] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)};
+	unsigned int numLines[] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
 
 	for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
 	{
@@ -243,8 +243,8 @@ bool Operator_Ext_Excitation::BuildExtension()
 			{
 				for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
 				{
-					curv->GetPoint(i-1,p1);
-					curv->GetPoint(i,p2);
+					curv->GetPoint(i-1,p1,m_Op->m_MeshType);
+					curv->GetPoint(i,p2,m_Op->m_MeshType);
 					path = m_Op->FindPath(p1,p2);
 					if (path.dir.size()>0)
 						prim->SetPrimitiveUsed(true);
diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
index b030c7e..723bc36 100644
--- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp
+++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
@@ -88,7 +88,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 	double dT = m_Op->GetTimestep();
 	unsigned int pos[] = {0,0,0};
 	double coord[3];
-	unsigned int numLines[3] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)};
+	unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
 	CSPropLorentzMaterial* mat = NULL;
 
 	double w_plasma,t_relax;
diff --git a/FDTD/extensions/operator_ext_mur_abc.cpp b/FDTD/extensions/operator_ext_mur_abc.cpp
index d1bbf70..a9697d2 100644
--- a/FDTD/extensions/operator_ext_mur_abc.cpp
+++ b/FDTD/extensions/operator_ext_mur_abc.cpp
@@ -104,12 +104,12 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
 	}
 	else
 	{
-		m_LineNr = m_Op->GetOriginalNumLines(m_ny)-1;
-		m_LineNr_Shift = m_Op->GetOriginalNumLines(m_ny) - 2;
+		m_LineNr = m_Op->GetNumberOfLines(m_ny,true)-1;
+		m_LineNr_Shift = m_Op->GetNumberOfLines(m_ny,true) - 2;
 	}
 
-	m_numLines[0] = m_Op->GetOriginalNumLines(m_nyP);
-	m_numLines[1] = m_Op->GetOriginalNumLines(m_nyPP);
+	m_numLines[0] = m_Op->GetNumberOfLines(m_nyP,true);
+	m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP,true);
 
 	m_Mur_Coeff_nyP = Create2DArray<FDTD_FLOAT>(m_numLines);
 	m_Mur_Coeff_nyPP = Create2DArray<FDTD_FLOAT>(m_numLines);
diff --git a/FDTD/extensions/operator_ext_tfsf.cpp b/FDTD/extensions/operator_ext_tfsf.cpp
index 1b743e0..7736b9f 100644
--- a/FDTD/extensions/operator_ext_tfsf.cpp
+++ b/FDTD/extensions/operator_ext_tfsf.cpp
@@ -176,7 +176,7 @@ bool Operator_Ext_TFST::BuildExtension()
 				m_ActiveDir[n][0]=false;
 			else
 				m_ActiveDir[n][0]=true;
-			if (m_Stop[n]==m_Op->GetOriginalNumLines(n)-1)
+			if (m_Stop[n]==m_Op->GetNumberOfLines(n,true)-1)
 				m_ActiveDir[n][1]=false;
 			else
 				m_ActiveDir[n][1]=true;
diff --git a/FDTD/extensions/operator_ext_upml.cpp b/FDTD/extensions/operator_ext_upml.cpp
index 29e28bf..487e04e 100644
--- a/FDTD/extensions/operator_ext_upml.cpp
+++ b/FDTD/extensions/operator_ext_upml.cpp
@@ -78,7 +78,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 
 	//check if mesh is large enough to support the pml
 	for (int n=0; n<3; ++n)
-		if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetOriginalNumLines(n) )
+		if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetNumberOfLines(n,true) )
 		{
 			cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl;
 			BC[2*n]=0;
@@ -131,7 +131,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 
 	Operator_Ext_UPML* op_ext_upml=NULL;
 	unsigned int start[3]={0 ,0 ,0};
-	unsigned int stop[3] ={op->GetOriginalNumLines(0)-1,op->GetOriginalNumLines(1)-1,op->GetOriginalNumLines(2)-1};
+	unsigned int stop[3] ={op->GetNumberOfLines(0,true)-1,op->GetNumberOfLines(1,true)-1,op->GetNumberOfLines(2,true)-1};
 
 	//create a pml in x-direction over the full width of yz-space
 	if (BC[0]==3)
@@ -148,8 +148,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 	{
 		op_ext_upml = new Operator_Ext_UPML(op);
 		op_ext_upml->SetGradingFunction(gradFunc);
-		start[0]=op->GetOriginalNumLines(0)-1-size[1];
-		stop[0] =op->GetOriginalNumLines(0)-1;
+		start[0]=op->GetNumberOfLines(0,true)-1-size[1];
+		stop[0] =op->GetNumberOfLines(0,true)-1;
 		op_ext_upml->SetBoundaryCondition(BC, size);
 		op_ext_upml->SetRange(start,stop);
 		op->AddExtension(op_ext_upml);
@@ -157,7 +157,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 
 	//create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions)
 	start[0]=(size[0]+1)*(BC[0]==3);
-	stop[0] =op->GetOriginalNumLines(0)-1-(size[0]+1)*(BC[1]==3);
+	stop[0] =op->GetNumberOfLines(0,true)-1-(size[0]+1)*(BC[1]==3);
 
 	if (BC[2]==3)
 	{
@@ -173,8 +173,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 	{
 		op_ext_upml = new Operator_Ext_UPML(op);
 		op_ext_upml->SetGradingFunction(gradFunc);
-		start[1]=op->GetOriginalNumLines(1)-1-size[3];
-		stop[1] =op->GetOriginalNumLines(1)-1;
+		start[1]=op->GetNumberOfLines(1,true)-1-size[3];
+		stop[1] =op->GetNumberOfLines(1,true)-1;
 		op_ext_upml->SetBoundaryCondition(BC, size);
 		op_ext_upml->SetRange(start,stop);
 		op->AddExtension(op_ext_upml);
@@ -182,7 +182,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 
 	//create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions)
 	start[1]=(size[2]+1)*(BC[2]==3);
-	stop[1] =op->GetOriginalNumLines(1)-1-(size[3]+1)*(BC[3]==3);
+	stop[1] =op->GetNumberOfLines(1,true)-1-(size[3]+1)*(BC[3]==3);
 
 	//exclude x-lines that does not belong to the base multi-grid operator;
 	Operator_CylinderMultiGrid* op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op);
@@ -203,8 +203,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 	{
 		op_ext_upml = new Operator_Ext_UPML(op);
 		op_ext_upml->SetGradingFunction(gradFunc);
-		start[2]=op->GetOriginalNumLines(2)-1-size[5];
-		stop[2] =op->GetOriginalNumLines(2)-1;
+		start[2]=op->GetNumberOfLines(2,true)-1-size[5];
+		stop[2] =op->GetNumberOfLines(2,true)-1;
 		op_ext_upml->SetBoundaryCondition(BC, size);
 		op_ext_upml->SetRange(start,stop);
 		op->AddExtension(op_ext_upml);
@@ -220,7 +220,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 		for (int n=0; n<2; ++n)
 		{
 			start[n]=0;
-			stop[n]=op_child->GetOriginalNumLines(n)-1;
+			stop[n]=op_child->GetNumberOfLines(n,true)-1;
 		}
 
 		if (op_cyl_MG)
@@ -240,8 +240,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
 		{
 			op_ext_upml = new Operator_Ext_UPML(op_child);
 			op_ext_upml->SetGradingFunction(gradFunc);
-			start[2]=op->GetOriginalNumLines(2)-1-size[5];
-			stop[2] =op->GetOriginalNumLines(2)-1;
+			start[2]=op->GetNumberOfLines(2,true)-1-size[5];
+			stop[2] =op->GetNumberOfLines(2,true)-1;
 			op_ext_upml->SetBoundaryCondition(BC, size);
 			op_ext_upml->SetRange(start,stop);
 			op_child->AddExtension(op_ext_upml);
@@ -320,10 +320,10 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm,
 			else
 				kappa_i[n] = 0;
 		}
-		else if ((pos[n] >= m_Op->GetOriginalNumLines(n) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3))  //upper pml in n-dir
+		else if ((pos[n] >= m_Op->GetNumberOfLines(n,true) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3))  //upper pml in n-dir
 		{
-			width = (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-m_Size[2*n+1]-1))*m_Op->GetGridDelta();
-			depth = width - (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta();
+			width = (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-m_Size[2*n+1]-1))*m_Op->GetGridDelta();
+			depth = width - (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta();
 
 			if ((m_Op_Cyl) && (n==1))
 			{
diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 6a05dcc..27cee11 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -227,14 +227,14 @@ double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
 	return GetNodeArea(ny, uiPos, dualMesh);
 }
 
-unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh) const
+unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
 {
 	inside = false;
 	if ((ny<0) || (ny>2))
 		return 0;
 	if (coord<GetDiscLine(ny,0))
 		return 0;
-	unsigned int numLines = GetNumberOfLines(ny);
+	unsigned int numLines = GetNumberOfLines(ny, fullMesh);
 	if (coord>GetDiscLine(ny,numLines-1))
 		return numLines-1;
 	inside=true;
@@ -258,24 +258,21 @@ unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool d
 	return 0;
 }
 
-bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool* inside) const
+bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool fullMesh, bool* inside) const
 {
 	bool meshInside=false;
 	bool ok=true;
 	for (int n=0; n<3; ++n)
 	{
-		uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh);
+		uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh,fullMesh);
 		ok &= meshInside;
 		if (inside)
 			inside[n]=meshInside;
 	}
-//	cerr << "Operator::SnapToMesh Wish: " << dcoord[0] << " " << dcoord[1] << " " << dcoord[2] << endl;
-//	cerr << "Operator::SnapToMesh Found: " << discLines[0][uicoord[0]] << " " << discLines[1][uicoord[1]] << " " << discLines[2][uicoord[2]] << endl;
-//	cerr << "Operator::SnapToMesh Index: " << uicoord[0] << " " << uicoord[1] << " " << uicoord[2] << endl;
 	return ok;
 }
 
-int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
+int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
 {
 	double l_start[3], l_stop[3];
 	for (int n=0;n<3;++n)
@@ -283,15 +280,15 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
 		l_start[n] = fmin(start[n],stop[n]);
 		l_stop[n] = fmax(start[n], stop[n]);
 		double min = GetDiscLine(n,0);
-		double max = GetDiscLine(n,GetNumberOfLines(n)-1);
+		double max = GetDiscLine(n,GetNumberOfLines(n, fullMesh)-1);
 		if ( ((l_start[n]<min) && (l_stop[n]<min)) || ((l_start[n]>max) && (l_stop[n]>max)) )
 		{
 			return -2;
 		}
 	}
 
-	SnapToMesh(l_start, uiStart, dualMesh, bStartIn);
-	SnapToMesh(l_stop, uiStop, dualMesh, bStopIn);
+	SnapToMesh(l_start, uiStart, dualMesh, fullMesh, bStartIn);
+	SnapToMesh(l_stop, uiStop, dualMesh, fullMesh, bStopIn);
 	int iDim = 0;
 
 	if (SnapMethod==0)
@@ -309,7 +306,7 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
 			{
 				if ((GetDiscLine( n, uiStart[n], dualMesh ) > l_start[n]) && (uiStart[n]>0))
 					--uiStart[n];
-				if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n)-1))
+				if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n, fullMesh)-1))
 					++uiStop[n];
 			}
 			if (uiStop[n]>uiStart[n])
@@ -323,7 +320,7 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
 		{
 			if (uiStop[n]>uiStart[n])
 			{
-				if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n)-1))
+				if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n, fullMesh)-1))
 					++uiStart[n];
 				if ((GetDiscLine( n, uiStop[n], dualMesh ) > l_stop[n]) && (uiStop[n]>0))
 					--uiStop[n];
@@ -341,24 +338,23 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
 struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
 {
 	struct Grid_Path path;
-//	double dV[] = {stop[0]-start[0],stop[1]-start[1],stop[2]-start[2]};
-
 	unsigned int uiStart[3],uiStop[3],currPos[3];
-	SnapToMesh(start,uiStart);
-	SnapToMesh(stop,uiStop);
+
+	SnapToMesh(start,uiStart,false,true);
+	SnapToMesh(stop,uiStop,false,true);
+
 	currPos[0]=uiStart[0];
 	currPos[1]=uiStart[1];
 	currPos[2]=uiStart[2];
-	double meshStart[] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]};
-	double meshStop[] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
-
+	double meshStart[3] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]};
+	double meshStop[3] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
 	bool UpDir = false;
 	double foot=0,dist=0,minFoot=0,minDist=0;
 	int minDir=0;
 	unsigned int minPos[3];
 	double startFoot,stopFoot,currFoot;
-	Point_Line_Distance(meshStart,start,stop,startFoot,dist);
-	Point_Line_Distance(meshStop,start,stop,stopFoot,dist);
+	Point_Line_Distance(meshStart,start,stop,startFoot,dist, m_MeshType);
+	Point_Line_Distance(meshStop,start,stop,stopFoot,dist, m_MeshType);
 	currFoot=startFoot;
 	minFoot=startFoot;
 	double P[3];
@@ -374,7 +370,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
 			if (((int)currPos[n]-1)>=0)
 			{
 				P[n] = discLines[n][currPos[n]-1];
-				Point_Line_Distance(P,start,stop,foot,dist);
+				Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
 				if ((foot>currFoot) && (dist<minDist))
 				{
 					minFoot=foot;
@@ -386,7 +382,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
 			if ((currPos[n]+1)<numLines[n])
 			{
 				P[n] = discLines[n][currPos[n]+1];
-				Point_Line_Distance(P,start,stop,foot,dist);
+				Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
 				if ((foot>currFoot) && (dist<minDist))
 				{
 					minFoot=foot;
@@ -408,6 +404,14 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
 			currPos[minDir]+=-1;
 			minPos[minDir]-=1;
 		}
+		//check validity of current postion
+		for (int n=0;n<3;++n)
+			if (currPos[n]>=numLines[n])
+			{
+				cerr << __func__ << ": Error, path went out of simulation domain, skipping path!" << endl;
+				Grid_Path empty;
+				return empty;
+			}
 		path.posPath[0].push_back(minPos[0]);
 		path.posPath[1].push_back(minPos[1]);
 		path.posPath[2].push_back(minPos[2]);
@@ -1324,7 +1328,7 @@ bool Operator::Calc_LumpedElements()
 				unsigned int uiStart[3];
 				unsigned int uiStop[3];
 				// snap to the native coordinate system
-				int Snap_Dimension = Operator::SnapBox2Mesh(box->GetStartCoord()->GetCoords(m_MeshType), box->GetStopCoord()->GetCoords(m_MeshType), uiStart, uiStop);
+				int Snap_Dimension = Operator::SnapBox2Mesh(box->GetStartCoord()->GetCoords(m_MeshType), box->GetStopCoord()->GetCoords(m_MeshType), uiStart, uiStop, false, true);
 				if (Snap_Dimension<=0)
 				{
 					if (Snap_Dimension>=-1)
@@ -1717,8 +1721,8 @@ void Operator::CalcPEC_Curves()
 			{
 				for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
 				{
-					curv->GetPoint(i-1,p1);
-					curv->GetPoint(i,p2);
+					curv->GetPoint(i-1,p1,m_MeshType);
+					curv->GetPoint(i,p2,m_MeshType);
 					path = FindPath(p1,p2);
 					if (path.dir.size()>0)
 						prim->SetPrimitiveUsed(true);
diff --git a/FDTD/operator.h b/FDTD/operator.h
index b164027..1901998 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -85,10 +85,7 @@ public:
 
 	virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();}
 
-	virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];}
-
-	//! Returns the number of lines as needed for the engine etc. (for post-processing etc, use GetNumLines())
-	virtual unsigned int GetOriginalNumLines(int ny) const {return numLines[ny];}
+	virtual unsigned int GetNumberOfLines(int ny, bool full=false) const {UNUSED(full);return numLines[ny];}
 
 	virtual void ShowStat() const;
 	virtual void ShowExtStat() const;
@@ -127,13 +124,13 @@ public:
 	*/
 	virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,pos,dualMesh);}
 
-	virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false) const;
+	virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const;
 
 	//! Snap the given coodinates to mesh indices
-	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const;
+	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const;
 
 	//! Snap a given box to the FDTD mesh
-	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
+	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
 
 	virtual void AddExtension(Operator_Extension* op_ext);
 	virtual void DeleteExtension(Operator_Extension* op_ext);
@@ -168,7 +165,7 @@ protected:
 		vector<unsigned int> posPath[3];
 		vector<unsigned short> dir;
 	};
-	struct Grid_Path FindPath(double start[], double stop[]);
+	virtual struct Grid_Path FindPath(double start[], double stop[]);
 
 	// debug
 	virtual void DumpOperator2File(string filename);
diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp
index d455196..58c437a 100644
--- a/FDTD/operator_cylinder.cpp
+++ b/FDTD/operator_cylinder.cpp
@@ -88,13 +88,16 @@ int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags )
 	return rc;
 }
 
-inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
+inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny, bool full) const
 {
+	if (full)
+		return Operator_Multithread::GetNumberOfLines(ny, full);
+
 	//this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1
 	if (CC_closedAlpha && ny==1)
-		return Operator_Multithread::GetNumberOfLines(ny)-2;
+		return Operator_Multithread::GetNumberOfLines(ny, true)-2;
 
-	return Operator_Multithread::GetNumberOfLines(ny);
+	return Operator_Multithread::GetNumberOfLines(ny, full);
 }
 
 string Operator_Cylinder::GetDirName(int ny) const
@@ -201,10 +204,10 @@ double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool du
 	return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh);
 }
 
-double Operator_Cylinder::FitToAlphaRange(double a_coord) const
+double Operator_Cylinder::FitToAlphaRange(double a_coord, bool fullMesh) const
 {
 	double min = GetDiscLine(1,0);
-	double max = GetDiscLine(1,GetOriginalNumLines(1)-1);
+	double max = GetDiscLine(1,GetNumberOfLines(1, fullMesh)-1);
 	if ((a_coord>=min) && (a_coord<=max))
 		return a_coord;
 	while (a_coord<min)
@@ -227,17 +230,17 @@ double Operator_Cylinder::FitToAlphaRange(double a_coord) const
 	return a_coord;
 }
 
-unsigned int Operator_Cylinder::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh) const
+unsigned int Operator_Cylinder::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
 {
 	if (ny==1)
 		coord=FitToAlphaRange(coord);
-	return Operator_Multithread::SnapToMeshLine(ny, coord, inside, dualMesh);
+	return Operator_Multithread::SnapToMeshLine(ny, coord, inside, dualMesh, fullMesh);
 }
 
-int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
+int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
 {
 	double a_min = GetDiscLine(1,0);
-	double a_max = GetDiscLine(1,GetNumberOfLines(1)-1);
+	double a_max = GetDiscLine(1,GetNumberOfLines(1,fullMesh)-1);
 
 	double a_size = stop[1] - start[1];
 	double a_center = FitToAlphaRange(0.5*(stop[1]+start[1]));
@@ -254,7 +257,7 @@ int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, uns
 
 	double l_start[3] = {start[0], a_start, start[2]};
 	double l_stop[3]  = {stop[0] , a_stop , stop[2] };
-	return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, SnapMethod, bStartIn, bStopIn);
+	return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, fullMesh, SnapMethod, bStartIn, bStopIn);
 }
 
 bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid)
diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h
index b6bea3f..e7a9a87 100644
--- a/FDTD/operator_cylinder.h
+++ b/FDTD/operator_cylinder.h
@@ -40,7 +40,7 @@ public:
 
 	virtual void ApplyMagneticBC(bool* dirs);
 
-	virtual unsigned int GetNumberOfLines(int ny) const;
+	virtual unsigned int GetNumberOfLines(int ny, bool full=false) const;
 
 	//! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z
 	virtual string GetDirName(int ny) const;
@@ -71,12 +71,12 @@ public:
 	*/
 	virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const;
 
-	virtual double FitToAlphaRange(double a_coord) const;
+	virtual double FitToAlphaRange(double a_coord, bool fullMesh=false) const;
 
-	virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false) const;
+	virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const;
 
 	//! Snap a given box to the FDTD mesh
-	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
+	virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
 
 	bool GetClosedAlpha() const {return CC_closedAlpha;}
 	bool GetR0Included() const {return CC_R0_included;}
diff --git a/FDTD/operator_mpi.cpp b/FDTD/operator_mpi.cpp
index 8aeeb0b..040c6d0 100644
--- a/FDTD/operator_mpi.cpp
+++ b/FDTD/operator_mpi.cpp
@@ -156,10 +156,13 @@ void Operator_MPI::SetOriginalMesh(CSRectGrid* orig_Mesh)
 	}
 }
 
-unsigned int Operator_MPI::GetNumberOfLines(int ny) const
+unsigned int Operator_MPI::GetNumberOfLines(int ny, bool fullMesh) const
 {
+	if (fullMesh)
+		return Operator_SSE_Compressed::GetNumberOfLines(ny,fullMesh);
+
 	if ((!m_MPI_Enabled) || (m_NeighborUp[ny]<0))
-		return Operator_SSE_Compressed::GetNumberOfLines(ny);
+		return Operator_SSE_Compressed::GetNumberOfLines(ny,fullMesh);
 
 	return Operator_SSE_Compressed::GetNumberOfLines(ny)-1;
 }
diff --git a/FDTD/operator_mpi.h b/FDTD/operator_mpi.h
index 341c6a6..a97f79b 100644
--- a/FDTD/operator_mpi.h
+++ b/FDTD/operator_mpi.h
@@ -50,7 +50,7 @@ public:
 	virtual void SetSplitPos(int ny, unsigned int pos) {m_SplitPos[ny]=pos;}
 	virtual void SetOriginalMesh(CSRectGrid* orig_Mesh);
 
-	virtual unsigned int GetNumberOfLines(int ny) const;
+	virtual unsigned int GetNumberOfLines(int ny, bool fullMesh=false) const;
 
 	virtual void AddExtension(Operator_Extension* op_ext);