diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index 800ccbce..2b226af9 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -970,10 +970,16 @@ void ASMbase::extractElmRes (const Matrix& globRes, Matrix& elmRes) const
 }
 
 
-void ASMbase::extractNodeVec (const Vector& globRes, Vector& nodeVec,
-			      const int* madof) const
+bool ASMbase::extractNodalVec (const Vector& globRes, Vector& nodeVec,
+                               const int* madof, int ngnod) const
 {
   nodeVec.clear();
+  if (nnod == 0)
+  {
+    std::cerr <<" *** ASMbase::extractNodalVec: Empty MADOF array."<< std::endl;
+    return false;
+  }
+
   nodeVec.reserve(nf*MLGN.size());
   for (size_t i = 0; i < MLGN.size(); i++)
   {
@@ -981,17 +987,67 @@ void ASMbase::extractNodeVec (const Vector& globRes, Vector& nodeVec,
     int idof = madof[inod-1] - 1;
     int jdof = madof[inod] - 1;
 #ifdef INDEX_CHECK
-    if (inod < 1 || jdof > (int)globRes.size())
-      std::cerr <<" *** ASMbase::extractNodeVec: Global DOF "<< jdof
+    bool ok = false;
+    if (inod < 1 || (inod > ngnod && ngnod > 0))
+      std::cerr <<" *** ASMbase::extractNodalVec: Global node "<< inod
+                <<" is out of range [1,"<< std::abs(ngnod) <<"]."<< std::endl;
+    else if (jdof > (int)globRes.size())
+      std::cerr <<" *** ASMbase::extractNodalVec: Global DOF "<< jdof
                 <<" is out of range [1,"<< globRes.size() <<"]."<< std::endl;
+    else
+      ok = true;
+    if (!ok) continue;
 #endif
     nodeVec.insert(nodeVec.end(),globRes.ptr()+idof,globRes.ptr()+jdof);
   }
+
+  return true;
+}
+
+
+bool ASMbase::injectNodalVec (const Vector& nodeVec, Vector& globVec,
+                              const std::vector<int>& madof, int basis) const
+{
+  if (madof.empty()) 
+  { 
+    std::cerr <<" *** ASMbase::injectNodalVec: Empty MADOF array."<< std::endl;
+    return false;
+  }
+
+  size_t ldof = 0;
+  char bType = basis == 1 ? 'D' : 'P'+basis-2;
+  for (size_t i = 0; i < MLGN.size(); i++)
+    if (basis == 0 || this->getNodeType(i+1) == bType)
+    {
+      int inod = MLGN[i];
+      int idof = madof[inod-1] - 1;
+      int ndof = madof[inod] - 1 - idof;
+#ifdef INDEX_CHECK
+      bool ok = false;
+      if (inod < 1 || inod > (int)madof.size())
+        std::cerr <<" *** ASMbase::injectNodalVec: Node "<< inod
+                  <<" is out of range [1,"<< madof.size() <<"]."<< std::endl;
+      else if (ldof+ndof > nodeVec.size())
+        std::cerr <<" *** ASMbase::injectNodalVec: Local DOF "<< ldof+ndof
+                  <<" is out of range [1,"<< nodeVec.size() <<"]"<< std::endl;
+      else if (idof+ndof > (int)globVec.size())
+        std::cerr <<" *** ASMbase::injectNodalVec: Global DOF "<< idof+ndof
+                  <<" is out of range [1,"<< globVec.size() <<"]"<< std::endl;
+      else
+        ok = true;
+      if (!ok) continue;
+#endif
+      std::copy(nodeVec.begin()+ldof, nodeVec.begin()+ldof+ndof,
+                globVec.begin()+idof);
+      ldof += ndof;
+    }
+
+  return true;
 }
 
 
 void ASMbase::extractNodeVec (const Vector& globRes, Vector& nodeVec,
-			      unsigned char nndof, int basis) const
+                              unsigned char nndof, int basis) const
 {
   if (nndof == 0) nndof = nf;
 
@@ -1019,7 +1075,7 @@ void ASMbase::extractNodeVec (const Vector& globRes, Vector& nodeVec,
 
 
 bool ASMbase::injectNodeVec (const Vector& nodeVec, Vector& globRes,
-			     unsigned char nndof, int) const
+                             unsigned char nndof, int) const
 {
   if (nndof == 0) nndof = nf;
 
diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h
index 7be24f4f..95f98c99 100644
--- a/src/ASM/ASMbase.h
+++ b/src/ASM/ASMbase.h
@@ -541,8 +541,9 @@ public:
   //! \param[in] globVec Global solution vector in DOF-order
   //! \param[out] nodeVec Nodal result vector for this patch
   //! \param[in] madof Global Matrix of Accumulated DOFs
-  void extractNodeVec(const Vector& globVec, Vector& nodeVec,
-                      const int* madof) const;
+  //! \param[in] ngnod Dimension of madof (the default -1 means unknown)
+  bool extractNodalVec(const Vector& globVec, Vector& nodeVec,
+                       const int* madof, int ngnod = -1) const;
 
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
@@ -560,6 +561,14 @@ public:
   virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
 			     unsigned char nndof = 0, int basis = 0) const;
 
+  //! \brief Injects nodal results for this patch into the global vector.
+  //! \param[in] nodeVec Nodal result vector for this patch
+  //! \param globVec Global solution vector in DOF-order
+  //! \param[in] madof Global Matrix of Accumulated DOFs
+  //! \param[in] basis Which basis to inject nodal values for (mixed methods)
+  bool injectNodalVec(const Vector& nodeVec, Vector& globVec,
+                      const std::vector<int>& madof, int basis = 0) const;
+
   //! \brief Creates and adds a two-point constraint to this patch.
   //! \param[in] slave Global node number of the node to constrain
   //! \param[in] dir Which local DOF to constrain (1, 2, 3)
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index b0f9a54f..3697bdc4 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -2502,7 +2502,7 @@ bool SIMbase::extractPatchSolution (IntegrandBase* problem,
   problem->initNodeMap(pch->getGlobalNodeNums());
   for (size_t i = 0; i < sol.size() && i < problem->getNoSolutions(); i++)
     if (!sol[i].empty())
-      pch->extractNodeVec(sol[i],problem->getSolution(i),mySam->getMADOF());
+      pch->extractNodalVec(sol[i],problem->getSolution(i),mySam->getMADOF());
 
   return this->extractPatchDependencies(problem,myModel,pindx);
 }
diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C
index e700594b..89d4e852 100644
--- a/src/SIM/SIMoutput.C
+++ b/src/SIM/SIMoutput.C
@@ -1097,7 +1097,7 @@ void SIMoutput::dumpPrimSol (const Vector& psol, utl::LogStream& os,
     if (myModel[i]->empty()) continue; // skip empty patches
 
     Vector patchSol;
-    myModel[i]->extractNodeVec(psol,patchSol,mySam->getMADOF());
+    myModel[i]->extractNodalVec(psol,patchSol,mySam->getMADOF());
 
     if (withID)
     {
@@ -1141,7 +1141,7 @@ bool SIMoutput::dumpSolution (const Vector& psol, utl::LogStream& os) const
     // Extract and write primary solution
     size_t nf = myModel[i]->getNoFields(1);
     Vector& patchSol = myProblem->getSolution();
-    myModel[i]->extractNodeVec(psol,patchSol,mySam->getMADOF());
+    myModel[i]->extractNodalVec(psol,patchSol,mySam->getMADOF());
     for (k = 0; k < nf; k++)
     {
       os << myProblem->getField1Name(k,"# FE");