diff --git a/Src/MultiGridOctreeData.IsoSurface.inl b/Src/MultiGridOctreeData.IsoSurface.inl
index 124ddd1..41c7e80 100644
--- a/Src/MultiGridOctreeData.IsoSurface.inl
+++ b/Src/MultiGridOctreeData.IsoSurface.inl
@@ -247,7 +247,7 @@ Real Octree< Real >::GetIsoValue( ConstPointer( Real ) solution , const std::vec
 	for( int d=maxDepth ; d>=_minDepth ; d-- )
 	{
 		std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-		for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( d );
+		for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( d );
 #pragma omp parallel for num_threads( threads ) reduction( + : isoValue , weightSum )
 		for( int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
 		{
@@ -279,8 +279,7 @@ Real Octree< Real >::GetIsoValue( ConstPointer( Real ) solution , const std::vec
 			if( w!=0 ) isoValue += value * w , weightSum += w;
 		}
 	}
-	if( _boundaryType==-1 ) return isoValue/weightSum - Real(0.5);
-	else                    return isoValue/weightSum;
+	return isoValue / weightSum;
 }
 
 template< class Real >
@@ -296,7 +295,7 @@ void Octree< Real >::SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPo
 {
 	typename Octree< Real >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice );
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ )
 	{
@@ -364,7 +363,7 @@ void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeig
 {
 	typename Octree< Real >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice );
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ )
 	{
@@ -456,7 +455,7 @@ void Octree< Real >::SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWei
 	typename Octree< Real >::template XSliceValues< Vertex >& xValues = slabValues[depth].xSliceValues( slab   );
 
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab+1] ; i++ )
 	{
@@ -661,7 +660,7 @@ void Octree< Real >::SetSliceIsoEdges( int depth , int slice , int z , std::vect
 {
 	typename Octree< Real >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice );
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ )
 	{
@@ -721,7 +720,7 @@ void Octree< Real >::SetXSliceIsoEdges( int depth , int slab , std::vector< Slab
 	typename Octree< Real >::template XSliceValues< Vertex >& xValues = slabValues[depth].xSliceValues( slab   );
 
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab+1] ; i++ )
 	{
@@ -797,7 +796,7 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 	std::vector< std::pair< int , Vertex > > polygon;
 	std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
 	std::vector< std::vector< IsoEdge > > edgess( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][offset] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][offset+1] ; i++ )
 	{
@@ -834,7 +833,7 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 							if( iter!=sValues.faceEdgeMap.end() )
 							{
 								const std::vector< IsoEdge >& _edges = iter->second;
-								for( int j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) );
+								for( size_t j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) );
 							}
 							else fprintf( stderr , "[ERROR] Invalid faces: %d  %d %d\n" , i , d , o ) , exit( 0 );
 						}
@@ -854,7 +853,7 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 							if( iter!=xValues.faceEdgeMap.end() )
 							{
 								const std::vector< IsoEdge >& _edges = iter->second;
-								for( int j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) );
+								for( size_t j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) );
 							}
 							else fprintf( stderr , "[ERROR] Invalid faces: %d  %d %d\n" , i , d , o ) , exit( 0 );
 						}
@@ -871,7 +870,7 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 					while( current!=start )
 					{
 						int idx;
-						for( idx=0 ; idx<edges.size() ; idx++ ) if( edges[idx][0]==current ) break;
+						for( idx=0 ; idx<(int)edges.size() ; idx++ ) if( edges[idx][0]==current ) break;
 						if( idx==edges.size() )
 						{
 							typename hash_map< long long , long long >::const_iterator iter;
@@ -890,10 +889,10 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 					loops.back().push_back( start );
 				}
 				// Add the loops to the mesh
-				for( int j=0 ; j<loops.size() ; j++ )
+				for( size_t j=0 ; j<loops.size() ; j++ )
 				{
 					std::vector< std::pair< int , Vertex > > polygon( loops[j].size() );
-					for( int k=0 ; k<loops[j].size() ; k++ )
+					for( size_t k=0 ; k<loops[j].size() ; k++ )
 					{
 						long long key = loops[j][k];
 						typename hash_map< long long , std::pair< int , Vertex > >::const_iterator iter;
diff --git a/Src/MultiGridOctreeData.SortedTreeNodes.inl b/Src/MultiGridOctreeData.SortedTreeNodes.inl
index 0fe89e1..4d5d6ec 100644
--- a/Src/MultiGridOctreeData.SortedTreeNodes.inl
+++ b/Src/MultiGridOctreeData.SortedTreeNodes.inl
@@ -151,7 +151,7 @@ void SortedTreeNodes::setSliceTableData( SliceTableData& sData , int depth , int
 	sData._cMap.resize( sData.nodeCount * Square::CORNERS , 0 ) , sData._eMap.resize( sData.nodeCount * Square::EDGES , 0 ) , sData._fMap.resize( sData.nodeCount * Square::FACES , 0 );
 	sData.cTable.resize( sData.nodeCount ) , sData.eTable.resize( sData.nodeCount ) , sData.fTable.resize( sData.nodeCount );
 	std::vector< TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=span.first ; i<span.second ; i++ )
 	{
@@ -246,9 +246,9 @@ void SortedTreeNodes::setSliceTableData( SliceTableData& sData , int depth , int
 	}
 	int cCount = 0 , eCount = 0 , fCount = 0;
 
-	for( int i=0 ; i<sData._cMap.size() ; i++ ) if( sData._cMap[i] ) sData._cMap[i] = cCount++;
-	for( int i=0 ; i<sData._eMap.size() ; i++ ) if( sData._eMap[i] ) sData._eMap[i] = eCount++;
-	for( int i=0 ; i<sData._fMap.size() ; i++ ) if( sData._fMap[i] ) sData._fMap[i] = fCount++;
+	for( size_t i=0 ; i<sData._cMap.size() ; i++ ) if( sData._cMap[i] ) sData._cMap[i] = cCount++;
+	for( size_t i=0 ; i<sData._eMap.size() ; i++ ) if( sData._eMap[i] ) sData._eMap[i] = eCount++;
+	for( size_t i=0 ; i<sData._fMap.size() ; i++ ) if( sData._fMap[i] ) sData._fMap[i] = fCount++;
 #pragma omp parallel for num_threads( threads )
 	for( int i=0 ; i<sData.nodeCount ; i++ )
 	{
@@ -272,7 +272,7 @@ void SortedTreeNodes::setXSliceTableData( XSliceTableData& sData , int depth , i
 	sData._eMap.resize( sData.nodeCount * Square::CORNERS , 0 ) , sData._fMap.resize( sData.nodeCount * Square::EDGES , 0 );
 	sData.eTable.resize( sData.nodeCount ) , sData.fTable.resize( sData.nodeCount );
 	std::vector< TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=span.first ; i<span.second ; i++ )
 	{
@@ -335,8 +335,8 @@ void SortedTreeNodes::setXSliceTableData( XSliceTableData& sData , int depth , i
 	}
 	int eCount = 0 , fCount = 0;
 
-	for( int i=0 ; i<sData._eMap.size() ; i++ ) if( sData._eMap[i] ) sData._eMap[i] = eCount++;
-	for( int i=0 ; i<sData._fMap.size() ; i++ ) if( sData._fMap[i] ) sData._fMap[i] = fCount++;
+	for( size_t i=0 ; i<sData._eMap.size() ; i++ ) if( sData._eMap[i] ) sData._eMap[i] = eCount++;
+	for( size_t i=0 ; i<sData._fMap.size() ; i++ ) if( sData._fMap[i] ) sData._fMap[i] = fCount++;
 #pragma omp parallel for num_threads( threads )
 	for( int i=0 ; i<sData.nodeCount ; i++ )
 	{
diff --git a/Src/MultiGridOctreeData.h b/Src/MultiGridOctreeData.h
index fe02875..784073c 100644
--- a/Src/MultiGridOctreeData.h
+++ b/Src/MultiGridOctreeData.h
@@ -30,7 +30,6 @@ DAMAGE.
 #define MULTI_GRID_OCTREE_DATA_INCLUDED
 
 #define NEW_CODE 1
-#define NEW_NEW_CODE 1
 
 //#define MAX_MEMORY_GB 15
 #define MAX_MEMORY_GB 0
@@ -71,8 +70,6 @@ DAMAGE.
 #define FORCE_NEUMANN_FIELD 1		// This flag forces the normal component across the boundary of the integration domain to be zero.
 									// This should be enabled if GRADIENT_DOMAIN_SOLUTION is not, so that CG doesn't run into trouble.
 
-#define ROBERTO_TOLDO_FIX 1
-
 #if !FORCE_NEUMANN_FIELD
 #pragma message( "[WARNING] Not zeroing out normal component on boundary" )
 #endif // !FORCE_NEUMANN_FIELD
@@ -216,16 +213,16 @@ public:
 	struct PointData
 	{
 		Point3D< Real > position;
-		Real weightedCoarserValue;
+		Real weightedCoarserDValue;
 		Real weight;
-		PointData( Point3D< Real > p=Point3D< Real >() , Real w=0 ) { position = p , weight = w , weightedCoarserValue = Real(0); }
+		PointData( Point3D< Real > p=Point3D< Real >() , Real w=0 ) { position = p , weight = w , weightedCoarserDValue = Real(0); }
 	};
 	template< class Data >
 	struct SparseNodeData
 	{
 		std::vector< int > indices;
 		std::vector< Data > data;
-		int index( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=indices.size() ? -1 : indices[ node->nodeData.nodeIndex ]; }
+		int index( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=(int)indices.size() ? -1 : indices[ node->nodeData.nodeIndex ]; }
 	};
 protected:
 	SortedTreeNodes _sNodes;
@@ -293,8 +290,8 @@ protected:
 	void SetPointValuesFromCoarser( SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients );
 	// Evalutes the solution @(depth) at the points @(depth-1) and updates the met constraints @(depth-1)
 	void SetPointConstraintsFromFiner ( const SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) metConstraints ) const;
-	Real _WeightedCoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) coarseCoefficients ) const;
-	Real _WeightedFinerFunctionValue  ( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real )  finerCoefficients ) const;
+	Real _CoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) coarseCoefficients ) const;
+	Real _FinerFunctionValue  ( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real )  finerCoefficients ) const;
 
 	// Down samples constraints @(depth) to constraints @(depth-1)
 	template< class C > void DownSample( int depth , const SortedTreeNodes& sNodes , ConstPointer( C ) fineConstraints    , Pointer( C ) coarseConstraints ) const;
diff --git a/Src/MultiGridOctreeData.inl b/Src/MultiGridOctreeData.inl
index ad02751..36d4e78 100644
--- a/Src/MultiGridOctreeData.inl
+++ b/Src/MultiGridOctreeData.inl
@@ -201,7 +201,7 @@ int Octree< Real >::SplatPointData( TreeOctNode* node , const Point3D< Real >& p
 			{
 				dxdydz = dxdy * dx[2][k];
 				TreeOctNode* _node = neighbors.neighbors[i][j][k];
-				if( dataInfo.indices.size()<TreeNodeData::NodeCount ) dataInfo.indices.resize( TreeNodeData::NodeCount , -1 );
+				if( (int)dataInfo.indices.size()<TreeNodeData::NodeCount ) dataInfo.indices.resize( TreeNodeData::NodeCount , -1 );
 				int idx = dataInfo.index( _node );
 				if( idx<0 )
 				{
@@ -499,7 +499,7 @@ template< class Real >
 int Octree< Real >::UpdateWeightContribution( std::vector< Real >& kernelDensityWeights , TreeOctNode* node , const Point3D<Real>& position , typename TreeOctNode::NeighborKey3& neighborKey , Real weight )
 {
 	typename TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
-	if( kernelDensityWeights.size()<TreeNodeData::NodeCount ) kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 );
+	if( (int)kernelDensityWeights.size()<TreeNodeData::NodeCount ) kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 );
 	double x , dxdy , dx[DIMENSION][3] , width;
 	Point3D< Real > center;
 	Real w;
@@ -724,7 +724,7 @@ int Octree< Real >::SetTree( OrientedPointStream< PointReal >* pointStream , int
 			myWidth = Real(1.0);
 			while( 1 )
 			{
-				if( pointInfo.indices.size()<TreeNodeData::NodeCount ) pointInfo.indices.resize( TreeNodeData::NodeCount , -1 );
+				if( (int)pointInfo.indices.size()<TreeNodeData::NodeCount ) pointInfo.indices.resize( TreeNodeData::NodeCount , -1 );
 				int idx = pointInfo.index( temp );
 
 				if( idx==-1 )
@@ -803,30 +803,30 @@ int Octree< Real >::SetTree( OrientedPointStream< PointReal >* pointStream , int
 		{
 			std::vector< int > temp = pointInfo.indices;
 			pointInfo.indices.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) pointInfo.indices[i] = temp[ indexMap[i] ];
-				else                          pointInfo.indices[i] = -1;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) pointInfo.indices[i] = temp[ indexMap[i] ];
+				else                               pointInfo.indices[i] = -1;
 		}
 		{
 			std::vector< int > temp = normalInfo.indices;
 			normalInfo.indices.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) normalInfo.indices[i] = temp[ indexMap[i] ];
-				else                          normalInfo.indices[i] = -1;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) normalInfo.indices[i] = temp[ indexMap[i] ];
+				else                               normalInfo.indices[i] = -1;
 		}
 		{
 			std::vector< Real > temp = centerWeights;
 			centerWeights.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) centerWeights[i] = temp[ indexMap[i] ];
-				else                          centerWeights[i] = (Real)0;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) centerWeights[i] = temp[ indexMap[i] ];
+				else                               centerWeights[i] = (Real)0;
 		}
 		{
 			std::vector< Real > temp = kernelDensityWeights;
 			kernelDensityWeights.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) kernelDensityWeights[i] = temp[ indexMap[i] ];
-				else                          kernelDensityWeights[i] = (Real)0;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) kernelDensityWeights[i] = temp[ indexMap[i] ];
+				else                               kernelDensityWeights[i] = (Real)0;
 		}
 	}
 	return cnt;
@@ -1028,7 +1028,7 @@ int Octree< Real >::SetTree( OrientedPointStreamWithData< PointReal , Data >* po
 			myWidth = Real(1.0);
 			while( 1 )
 			{
-				if( pointInfo.indices.size()<TreeNodeData::NodeCount ) pointInfo.indices.resize( TreeNodeData::NodeCount , -1 );
+				if( (int)pointInfo.indices.size()<TreeNodeData::NodeCount ) pointInfo.indices.resize( TreeNodeData::NodeCount , -1 );
 				int idx = pointInfo.index( temp );
 
 #if POINT_DATA_RES
@@ -1129,37 +1129,37 @@ int Octree< Real >::SetTree( OrientedPointStreamWithData< PointReal , Data >* po
 		{
 			std::vector< int > temp = pointInfo.indices;
 			pointInfo.indices.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) pointInfo.indices[i] = temp[ indexMap[i] ];
-				else                          pointInfo.indices[i] = -1;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) pointInfo.indices[i] = temp[ indexMap[i] ];
+				else                               pointInfo.indices[i] = -1;
 		}
 		{
 			std::vector< int > temp = normalInfo.indices;
 			normalInfo.indices.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) normalInfo.indices[i] = temp[ indexMap[i] ];
-				else                          normalInfo.indices[i] = -1;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) normalInfo.indices[i] = temp[ indexMap[i] ];
+				else                               normalInfo.indices[i] = -1;
 		}
 		{
 			std::vector< Real > temp = centerWeights;
 			centerWeights.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) centerWeights[i] = temp[ indexMap[i] ];
-				else                          centerWeights[i] = (Real)0;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) centerWeights[i] = temp[ indexMap[i] ];
+				else                               centerWeights[i] = (Real)0;
 		}
 		{
 			std::vector< Real > temp = kernelDensityWeights;
 			kernelDensityWeights.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) kernelDensityWeights[i] = temp[ indexMap[i] ];
-				else                          kernelDensityWeights[i] = (Real)0;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) kernelDensityWeights[i] = temp[ indexMap[i] ];
+				else                               kernelDensityWeights[i] = (Real)0;
 		}
 		{
 			std::vector< int > temp = dataValues.indices;
 			dataValues.indices.resize( indexMap.size() );
-			for( int i=0 ; i<indexMap.size() ; i++ )
-				if( indexMap[i]<temp.size() ) dataValues.indices[i] = temp[ indexMap[i] ];
-				else                          dataValues.indices[i] = -1;
+			for( size_t i=0 ; i<indexMap.size() ; i++ )
+				if( indexMap[i]<(int)temp.size() ) dataValues.indices[i] = temp[ indexMap[i] ];
+				else                               dataValues.indices[i] = -1;
 		}
 	}
 	return cnt;
@@ -1346,15 +1346,9 @@ int Octree< Real >::SetMatrixRow( const SparseNodeData< PointData >& pointInfo ,
 						Point3D< Real > p = pData.position;
 						for( int s=0 ; s<3 ; s++ )
 						{
-#if ROBERTO_TOLDO_FIX
 							if( idx[0]+j-s>=0 && idx[0]+j-s<((2<<d)-1) ) _splineValues[3*0+s] = Real( _fData.baseBSplines[ idx[0]+j-s][s]( p[0] ) );
 							if( idx[1]+k-s>=0 && idx[1]+k-s<((2<<d)-1) ) _splineValues[3*1+s] = Real( _fData.baseBSplines[ idx[1]+k-s][s]( p[1] ) );
 							if( idx[2]+l-s>=0 && idx[2]+l-s<((2<<d)-1) ) _splineValues[3*2+s] = Real( _fData.baseBSplines[ idx[2]+l-s][s]( p[2] ) );
-#else // !ROBERTO_TOLDO_FIX
-							_splineValues[3*0+s] = Real( _fData.baseBSplines[ idx[0]+j-s][s]( p[0] ) );
-							_splineValues[3*1+s] = Real( _fData.baseBSplines[ idx[1]+k-s][s]( p[1] ) );
-							_splineValues[3*2+s] = Real( _fData.baseBSplines[ idx[2]+l-s][s]( p[2] ) );
-#endif // ROBERTO_TOLDO_FIX
 						}
 						Real value = _splineValues[3*0+j] * _splineValues[3*1+k] * _splineValues[3*2+l];
 						Real weightedValue = value * weight;
@@ -1666,7 +1660,7 @@ void Octree< Real >::UpdateConstraintsFromFiner( const typename BSplineData< 2 >
 
 	// Iterate over the nodes @( depth )
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
 #pragma omp parallel for num_threads( threads )
 	for( int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
 	{
@@ -1761,13 +1755,12 @@ void Octree< Real >::UpdateConstraintsFromCoarser( const SparseNodeData< PointDa
 			if( neighbors5.neighbors[x][y][z] && pointInfo.index( neighbors5.neighbors[x][y][z] )!=-1 )
 			{
 				const PointData& pData = points[ pointInfo.index( neighbors5.neighbors[x][y][z] ) ];
-				Real weightedPointValue = pData.weightedCoarserValue;
 				Point3D< Real > p = pData.position;
 				constraint += 
 					_fData.baseBSplines[idx[0]][x-1]( p[0] ) *
 					_fData.baseBSplines[idx[1]][y-1]( p[1] ) *
 					_fData.baseBSplines[idx[2]][z-1]( p[2] ) * 
-					weightedPointValue;
+					pData.weightedCoarserDValue;
 			}
 		constraints[ node->nodeData.nodeIndex ] -= Real( constraint );
 	}
@@ -1789,7 +1782,7 @@ void Octree< Real >::DownSample( int depth , const SortedTreeNodes& sNodes , Con
 	else if( _boundaryType== 1 ) cornerValue = 1.00;
 	else                         cornerValue = 0.75;
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
 	{
@@ -1837,7 +1830,7 @@ void Octree< Real >::UpSample( int depth , const SortedTreeNodes& sNodes , Const
 	if( depth<=_minDepth ) return;
 
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
 #pragma omp parallel for num_threads( threads )
 	for( int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
 	{
@@ -1885,7 +1878,7 @@ void Octree< Real >::SetPointValuesFromCoarser( SparseNodeData< PointData >& poi
 	std::vector< PointData >& points = pointInfo.data;
 	// For every node at the current depth
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
 	for( int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
 	{
@@ -1894,20 +1887,18 @@ void Octree< Real >::SetPointValuesFromCoarser( SparseNodeData< PointData >& poi
 		if( pIdx!=-1 )
 		{
 			neighborKey.getNeighbors( sNodes.treeNodes[i] );
-			points[ pIdx ].weightedCoarserValue = _WeightedCoarserFunctionValue( points[pIdx] , neighborKey , sNodes.treeNodes[i] , coarseCoefficients-_sNodes.nodeCount[depth-1] );
+			points[ pIdx ].weightedCoarserDValue = (Real)( _CoarserFunctionValue( points[pIdx] , neighborKey , sNodes.treeNodes[i] , coarseCoefficients-_sNodes.nodeCount[depth-1] ) - 0.5 ) * points[pIdx].weight;
 		}
 	}
 }
 template< class Real >
-Real Octree< Real >::_WeightedCoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) coarseCoefficients ) const
+Real Octree< Real >::_CoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) coarseCoefficients ) const
 {
 	double pointValue = 0;
 	int depth = pointNode->depth();
-	if( _boundaryType==-1 && depth==0 ) return Real(-0.5) * pointData.weight;
 
 	if( depth<=_minDepth ) return Real(0.);
 
-	Real weight       = pointData.weight;
 	Point3D< Real > p = pointData.position;
 
 	// Iterate over all basis functions that overlap the point at the coarser resolutions
@@ -1921,40 +1912,26 @@ Real Octree< Real >::_WeightedCoarserFunctionValue( const PointData& pointData ,
 
 		for( int j=0 ; j<3 ; j++ )
 		{
-#if ROBERTO_TOLDO_FIX
 			double xValue = 0;
 			if( _idx[0]+j>=0 && _idx[0]+j<((1<<depth)-1) ) xValue = _fData.baseBSplines[ _idx[0]+j ][2-j]( p[0] );
 			else continue;
-#else // !ROBERTO_TOLDO_FIX
-			double xValue = _fData.baseBSplines[ _idx[0]+j ][2-j]( p[0] );
-#endif // ROBERTO_TOLDO_FIX
 			for( int k=0 ; k<3 ; k++ )
 			{
-#if ROBERTO_TOLDO_FIX
 				double xyValue = 0;
 				if( _idx[1]+k>=0 && _idx[1]+k<((1<<depth)-1) ) xyValue = xValue * _fData.baseBSplines[ _idx[1]+k ][2-k]( p[1] );
 				else continue;
-#else // !ROBERTO_TOLDO_FIX
-				double xyValue = xValue * _fData.baseBSplines[ _idx[1]+k ][2-k]( p[1] );
-#endif // ROBERTO_TOLDO_FIX
 				double _pointValue = 0;
 				for( int l=0 ; l<3 ; l++ )
 				{
 					const TreeOctNode* basisNode = neighbors.neighbors[j][k][l];
-#if ROBERTO_TOLDO_FIX
 					if( basisNode && basisNode->nodeData.nodeIndex>=0 && _idx[2]+l>=0 && _idx[2]+l<((1<<depth)-1) )
 						_pointValue += _fData.baseBSplines[ _idx[2]+l ][2-l]( p[2] ) * double( coarseCoefficients[basisNode->nodeData.nodeIndex] );
-#else // !ROBERTO_TOLDO_FIX
-					if( basisNode && basisNode->nodeData.nodeIndex>=0 )
-						_pointValue += _fData.baseBSplines[ _idx[2]+l ][2-l]( p[2] ) * double( coarseCoefficients[basisNode->nodeData.nodeIndex] );
-#endif // ROBERTO_TOLDO_FIX
 				}
 				pointValue += _pointValue * xyValue;
 			}
 		}
 	}
-	if( _boundaryType==-1 ) pointValue -= 0.5;
-	return Real( pointValue * weight );
+	return Real( pointValue );
 }
 template< class Real >
 void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) coarserConstraints ) const
@@ -1967,7 +1944,7 @@ void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointDa
 	size_t start = sNodes.nodeCount[depth-1] , end = sNodes.nodeCount[depth] , range = end-start;
 	memset( coarserConstraints , 0 , sizeof( Real ) * ( sNodes.nodeCount[depth]-sNodes.nodeCount[depth-1] ) );
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth-1 );
 #pragma omp parallel for num_threads( threads )
 	for( int i=sNodes.nodeCount[depth-1] ; i<sNodes.nodeCount[depth] ; i++ )
 	{
@@ -1978,7 +1955,7 @@ void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointDa
 			typename TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( sNodes.treeNodes[i] );
 			// Evaluate the solution @( depth ) at the current point @( depth-1 )
 			{
-				Real finerPointValue = _WeightedFinerFunctionValue( points[pIdx] , neighborKey , sNodes.treeNodes[i] , finerCoefficients-sNodes.nodeCount[depth] );
+				Real finerPointDValue = (Real)( _FinerFunctionValue( points[pIdx] , neighborKey , sNodes.treeNodes[i] , finerCoefficients-sNodes.nodeCount[depth] ) - 0.5 ) * points[pIdx].weight;
 				Point3D< Real > p = points[ pIdx ].position;
 				// Update constraints for all nodes @( depth-1 ) that overlap the point
 				int d , idx[3];
@@ -1997,7 +1974,7 @@ void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointDa
 							_fData.baseBSplines[idx[0]+x][2-x]( p[0] ) *
 							_fData.baseBSplines[idx[1]+y][2-y]( p[1] ) *
 							_fData.baseBSplines[idx[2]+z][2-z]( p[2] ) * 
-							finerPointValue
+							finerPointDValue
 							);
 					}
 			}
@@ -2005,7 +1982,7 @@ void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointDa
 	}
 }
 template< class Real >
-Real Octree< Real >::_WeightedFinerFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) finerCoefficients ) const
+Real Octree< Real >::_FinerFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) finerCoefficients ) const
 {
 	typename TreeOctNode::Neighbors3 childNeighbors;
 	double pointValue = 0;
@@ -2042,8 +2019,7 @@ Real Octree< Real >::_WeightedFinerFunctionValue( const PointData& pointData , c
 			_fData.baseBSplines[ idx[2]+l ][2-l]( p[2] ) *
 			double( finerCoefficients[ basisNode->nodeData.nodeIndex ] );
 	}
-	if( _boundaryType==-1 ) pointValue -= Real(0.5);
-	return Real( pointValue * weight );
+	return Real( pointValue );
 }
 template< class Real >
 int Octree< Real >::GetSliceMatrixAndUpdateConstraints( const SparseNodeData< PointData >& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd )
@@ -2054,9 +2030,9 @@ int Octree< Real >::GetSliceMatrixAndUpdateConstraints( const SparseNodeData< Po
 	SetLaplacianStencils( depth , integrator , stencils );
 	matrix.Resize( (int)range );
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
-	for( int i=0 ; i<range ; i++ )
+	for( int i=0 ; i<(int)range ; i++ )
 	{
 		typename TreeOctNode::NeighborKey3& neighborKey = neighborKeys[ omp_get_thread_num() ];
 		TreeOctNode* node = sNodes.treeNodes[i+nStart];
@@ -2109,9 +2085,9 @@ int Octree< Real >::GetMatrixAndUpdateConstraints( const SparseNodeData< PointDa
 	SetLaplacianStencils( depth , integrator , stencils );
 	matrix.Resize( (int)range );
 	std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-	for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
+	for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth );
 #pragma omp parallel for num_threads( threads )
-	for( int i=0 ; i<range ; i++ )
+	for( int i=0 ; i<(int)range ; i++ )
 	{
 		typename TreeOctNode::NeighborKey3& neighborKey = neighborKeys[ omp_get_thread_num() ];
 		TreeOctNode* node = sNodes.treeNodes[i+start];
@@ -2502,7 +2478,7 @@ Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const SparseNodeData< P
 		SetDivergenceStencils( d , integrator , stencils , true );
 
 		std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-		for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( _fData.depth );
+		for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( _fData.depth );
 #pragma omp parallel for num_threads( threads )
 		for( int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
 		{
@@ -2630,7 +2606,7 @@ Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const SparseNodeData< P
 		Stencil< Point3D< double > , 5 > stencils[2][2][2];
 		SetDivergenceStencils( d , integrator , stencils , false );
 		std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) );
-		for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( maxDepth );
+		for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( maxDepth );
 #pragma omp parallel for num_threads( threads )
 		for( int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
 		{
@@ -2750,9 +2726,6 @@ template< class V >
 V Octree< Real >::getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey , const TreeOctNode* node , int corner , ConstPointer( V ) solution , ConstPointer( V ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 > stencils[8] , bool isInterior ) const
 {
 	V value(0);
-	// [NOTE] For screening, we force the samples to have value 0.0, so Dirichlet boundary constraints are implemented
-	// by having the points on the boundary have value -0.5. Hence the offsetting.
-	if( _boundaryType==-1 ) value = -0.5;
 	int d , off[3];
 	node->depthAndOffset( d , off );
 
@@ -2819,7 +2792,6 @@ std::pair< Real , Point3D< Real > > Octree< Real >::getCornerValueAndNormal( con
 {
 	double value = 0;
 	Point3D< double > normal;
-	if( _boundaryType==-1 ) value = -0.5;
 	int d , off[3];
 	node->depthAndOffset( d , off );
 
@@ -2989,26 +2961,30 @@ V Octree< Real >::_Evaluate( const SparseNodeData< V >& coefficients , Point3D<
 {
 	V value = V(0);
 
-	for( int d=0 ; d<=neighborKey3.depth() ; d++ ) for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) for( int k=0 ; k<3 ; k++ )
+	double functionValues[3][3];
+	for( int d=0 ; d<=neighborKey3.depth() ; d++ )
 	{
-		const TreeOctNode* n = neighborKey3.neighbors[d].neighbors[i][j][k];
-		if( n )
 		{
+			const TreeOctNode* n = neighborKey3.neighbors[d].neighbors[1][1][1];
+			if( !n ) fprintf( stderr , "[ERROR] _Evaluate: Center node needs to exist\n" ) , exit( 0 );
 			int dd , off[3];
 			n->depthAndOffset( dd , off );
-	
-			int idx = coefficients.index( n );
-			if( idx>=0 )
-				value += 
-				(
-				coefficients.data[idx] *
-				(Real)
-				(
-					_fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[0] ) ]( p[0] ) *
-					_fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[1] ) ]( p[1] ) *
-					_fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[2] ) ]( p[2] )
-					)
-				);
+			for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) if( off[i]-1+j>=0 && off[i]-1+j<(1<<dd) )
+				functionValues[i][j] = _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[i]-1+j ) ]( p[i] );
+		}
+
+		for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ )
+		{
+			double functionValue = functionValues[0][i] * functionValues[1][j];
+			for( int k=0 ; k<3 ; k++ )
+			{
+				const TreeOctNode* n = neighborKey3.neighbors[d].neighbors[i][j][k];
+				if( n )
+				{
+					int idx = coefficients.index( n );
+					if( idx>=0 ) value += coefficients.data[idx] * (Real)( functionValue * functionValues[2][k] );
+				}
+			}
 		}
 	}
 	return value;
@@ -3143,9 +3119,6 @@ Pointer( V ) Octree< Real >::Evaluate( ConstPointer( V ) coefficients , int& res
 						vTables.valueTable[ idx[2] + z*vTables.functionCount ];
 				}
 	}
-	// [NOTE] For screening, we force the samples to have value 0.0, so Dirichlet boundary constraints are implemented
-	// by having the points on the boundary have value -0.5. Hence the offsetting.
-	if( _boundaryType==-1 ) for( int i=0 ; i<res*res*res ; i++ ) values[i] -= Real(0.5);
 	for( int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
 
 	return values;
diff --git a/Src/SparseMatrix.inl b/Src/SparseMatrix.inl
index 9f9874b..3a8b1eb 100644
--- a/Src/SparseMatrix.inl
+++ b/Src/SparseMatrix.inl
@@ -265,8 +265,7 @@ int SparseMatrix<T>::SolveSymmetric( const SparseMatrix<T>& M , const Vector<T2>
 	M.Multiply( solution , r );
 	r = b - r;
 	Vector< T2 > d = r;
-	double delta_new = 0.0;
-	double delta_0 = 0.0;
+	double delta_new = 0. , delta_0 = 0.;
 	for( int i=0 ; i<r.Dimensions() ; i++ ) delta_new += r.m_pV[i] * r.m_pV[i];
 	delta_0 = delta_new;
 	if( delta_new<eps ) return 0;