diff --git a/FECore/FESPRProjection.cpp b/FECore/FESPRProjection.cpp
index 68225a62d..b131e239a 100644
--- a/FECore/FESPRProjection.cpp
+++ b/FECore/FESPRProjection.cpp
@@ -64,13 +64,8 @@ void FESPRProjection::Project(FESolidDomain& dom, const vector< vector<double> >
 	{
 	case ET_TET4   :
 	case ET_TET5   : { NDOF =  4; NCN = 4; } break;
-	case ET_TET10  : { NDOF =  4; NCN = 4; } break;
-	case ET_TET15  : 
-		{
-			NDOF = (m_p == 1 ? 4 : 10); 
-			NCN = 4; 
-		}
-		break;
+	case ET_TET10  : { NDOF = (m_p == 1 ? 4 : 10); NCN = 4;	} break;
+	case ET_TET15  : { NDOF = (m_p == 1 ? 4 : 10); NCN = 4; } break;
 	case ET_TET20 : { NDOF = 10; NCN = 4; } break;
 	case ET_HEX8  : { NDOF =  7; NCN = 8; } break;
 	case ET_HEX20 : { NDOF = (m_p == 1 ? 7 : 10); NCN = 8; } break;
diff --git a/FECore/writeplot.cpp b/FECore/writeplot.cpp
index ff39bad22..696b48258 100644
--- a/FECore/writeplot.cpp
+++ b/FECore/writeplot.cpp
@@ -73,6 +73,68 @@ void writeSPRElementValue(FESolidDomain& dom, FEDataStream& ar, std::function<do
 	}
 }
 
+//-------------------------------------------------------------------------------------------------
+void writeSPRElementValueVectorDouble(FESolidDomain& dom, FEDataStream& ar, std::function<std::vector<double>(const FEMaterialPoint&)> fnc, int interpolOrder, int n_fields)
+{
+
+	// get all nodes and elements
+	int NN = dom.Nodes();
+	int NE = dom.Elements();
+
+	// build the element data array
+	vector<vector< vector<double> > > ED(n_fields);
+	// for each component
+	for (int n = 0; n < n_fields; ++n)
+	{
+		// fill element data. for each element add the number of gauss points in the element.
+		ED[n].resize(NE);
+		for (int i = 0; i < NE; ++i)
+		{
+			FESolidElement& e = dom.Element(i);
+			int nint = e.GaussPoints();
+
+			ED[n][i].assign(nint, 0.0);
+		}
+	}
+
+	// this array will store the results
+	FESPRProjection map;
+	map.SetInterpolationOrder(interpolOrder);
+	vector<vector<double> >val(n_fields);
+
+	// fill the ED array
+	for (int i = 0; i < NE; ++i)
+	{
+		FESolidElement& el = dom.Element(i);
+		int nint = el.GaussPoints();
+		for (int j = 0; j < nint; ++j)
+		{
+			FEMaterialPoint& mp = *el.GetMaterialPoint(j);
+			vector<double> v = fnc(mp);
+
+			// loop over all solutes components
+			for (int n = 0; n < n_fields; ++n)
+			{
+				ED[n][i][j] = v[n];
+			}
+		}
+	}
+
+	// project to nodes
+	// loop over stress components
+	for (int n = 0; n < n_fields; ++n)
+	{
+		map.Project(dom, ED[n], val[n]);
+	}
+
+	// copy results to archive
+	for (int i_node = 0; i_node < NN; ++i_node)
+	{
+		for (int i_field = 0; i_field < n_fields; ++i_field)
+			ar.push_back((float)val[i_field][i_node]);
+	}
+}
+
 void writeSPRElementValueMat3dd(FESolidDomain& dom, FEDataStream& ar, std::function<mat3dd(const FEMaterialPoint&)> fnc, int interpolOrder)
 {
 	int NN = dom.Nodes();
diff --git a/FECore/writeplot.h b/FECore/writeplot.h
index b8c3e8f21..92867a0de 100644
--- a/FECore/writeplot.h
+++ b/FECore/writeplot.h
@@ -245,6 +245,7 @@ template <class T> void writeNodalProjectedElementValues(FESurface& dom, FEDataS
 // helper functions for writing SPR projected element values
 // TODO: I needed to give these functions a different name because of the implicit conversion between mat3ds and mat3dd
 FECORE_API void writeSPRElementValue(FESolidDomain& dom, FEDataStream& ar, std::function<double(const FEMaterialPoint&)> fnc, int interpolOrder = -1);
+FECORE_API void writeSPRElementValueVectorDouble(FESolidDomain& dom, FEDataStream& ar, std::function<std::vector<double>(const FEMaterialPoint&)> fnc, int interpolOrder = -1, int n_fields = -1);
 FECORE_API void writeSPRElementValueMat3dd(FESolidDomain& dom, FEDataStream& ar, std::function<mat3dd(const FEMaterialPoint&)> fnc, int interpolOrder = -1);
 FECORE_API void writeSPRElementValueMat3ds(FESolidDomain& dom, FEDataStream& ar, std::function<mat3ds(const FEMaterialPoint&)> fnc, int interpolOrder = -1);