diff --git a/src/ASM/NewmarkMats.C b/src/ASM/NewmarkMats.C
index a9d79bbe..dd211cf3 100644
--- a/src/ASM/NewmarkMats.C
+++ b/src/ASM/NewmarkMats.C
@@ -15,15 +15,26 @@
 #include "NewmarkMats.h"
 
 
-NewmarkMats::NewmarkMats (double a1, double a2, double b, double c)
+NewmarkMats::NewmarkMats (double a1, double a2, double b, double c,
+                          bool generalizedAlpha) : isPredictor(true), h(0.0)
 {
   alpha1 = a1;
   alpha2 = a2;
-  beta   = b;
-  gamma  = c;
 
-  isPredictor = true;
-  h = 0.0;
+  if (generalizedAlpha)
+  {
+    alpha_m = b;
+    alpha_f = c;
+    double alpha = alpha_f - alpha_m;
+    beta = 0.25*(1.0-alpha)*(1.0-alpha);
+    gamma = 0.5 - alpha;
+  }
+  else
+  {
+    alpha_m = alpha_f = 1.0;
+    beta  = b;
+    gamma = c;
+  }
 }
 
 
@@ -32,10 +43,8 @@ const Matrix& NewmarkMats::getNewtonMatrix () const
   Matrix& N = const_cast<Matrix&>(A.front());
 
   N = A[1];
-  if (alpha1 > 0.0)
-    N.multiply(1.0 + alpha1*gamma*h);    // [N]  = (1+alpha1*gamma*h)*[M]
-
-  N.add(A[2],(alpha2*gamma + beta*h)*h); // [N] += (alpha2*gamma+beta*h)*h*[K]
+  N.multiply(alpha_m + alpha_f*alpha1*gamma*h);
+  N.add(A[2],alpha_f*(alpha2*gamma + beta*h)*h);
 
 #if SP_DEBUG > 2
   std::cout <<"\nElement mass matrix"<< A[1];
@@ -56,13 +65,12 @@ const Vector& NewmarkMats::getRHSVector () const
     int ia = vec.size() - 1; // index to element acceleration vector (a)
     int iv = vec.size() - 2; // index to element velocity vector (v)
 #if SP_DEBUG > 2
-    std::cout <<"\nf_ext"<< dF;
+    std::cout <<"\nf_ext - f_s"<< dF;
     std::cout <<"f_i = M*a"<< A[1]*vec[ia];
     if (alpha1 > 0.0)
       std::cout <<"f_d1/alpha1 = M*v (alpha1="<< alpha1 <<")"<< A[1]*vec[iv];
     if (alpha2 > 0.0)
       std::cout <<"f_d2/alpha2 = K*v (alpha2="<< alpha2 <<")"<< A[2]*vec[iv];
-    std::cout <<"f_s = K*d"<< A[2]*vec.front();
 #endif
 
     dF.add(A[1]*vec[ia],-1.0);      // dF = Fext - M*a
@@ -72,8 +80,6 @@ const Vector& NewmarkMats::getRHSVector () const
 
     if (alpha2 > 0.0)
       dF.add(A[2]*vec[iv],-alpha2); // dF -= alpha2*K*v
-
-    dF.add(A[2]*vec.front(),-1.0);  // dF -= K*d
   }
 
 #if SP_DEBUG > 2
diff --git a/src/ASM/NewmarkMats.h b/src/ASM/NewmarkMats.h
index 5f20ad7a..620cb70a 100644
--- a/src/ASM/NewmarkMats.h
+++ b/src/ASM/NewmarkMats.h
@@ -26,11 +26,21 @@ class NewmarkMats : public ElmMats
 {
 public:
   //! \brief The constructor initializes the time integration parameters.
-  NewmarkMats(double a1 = 0.0, double a2 = 0.0, double b = 0.0, double c = 0.0);
+  //! \param[in] a1 Mass-proportional damping coefficient
+  //! \param[in] a2 Stiffness-proportional damping coefficient
+  //! \param[in] b Time integration parameter
+  //! \param[in] c Time integration parameter
+  //! \param[in] generalizedAlpha If \e true, iterpret \a b and \a c as the
+  //! generalized alpha parameters \a alpha_m and \a alpha_f, respectively,
+  //! otherwise as \a beta and \a gamma
+  NewmarkMats(double a1 = 0.0, double a2 = 0.0, double b = 0.0, double c = 0.0,
+              bool generalizedAlpha = false);
   //! \brief Empty destructor.
   virtual ~NewmarkMats() {}
 
   //! \brief Updates the time step size and the \a isPredictor flag.
+  //! \param[in] dt New time step size
+  //! \param[in] iter Newton-Raphson iteration counter
   void setStepSize(double dt, int iter) { h = dt; isPredictor = iter == 0; }
 
   //! \brief Returns the element-level Newton matrix.
@@ -39,13 +49,17 @@ public:
   virtual const Vector& getRHSVector() const;
 
 protected:
-  double alpha1; //!< Mass-proportional damping
-  double alpha2; //!< Stiffness-proportional damping
-  double beta;   //!< Time integration parameter
-  double gamma;  //!< Time integration parameter
-  double h;      //!< Time step size
+  bool  isPredictor; //!< If \e true, we are in the predictor step
+  double h;          //!< Time step size
 
-  bool isPredictor; //!< If \e true, we are in the predictor step
+  double alpha1; //!< Mass-proportional damping coefficient
+  double alpha2; //!< Stiffness-proportional damping coefficient
+  double beta;   //!< Newmark time integration parameter
+  double gamma;  //!< Newmark time integration parameter
+
+private:
+  double alpha_m; //!< Generalized-alpha parameter
+  double alpha_f; //!< Generalized-alpha parameter
 };
 
 #endif
diff --git a/src/SIM/GenAlphaSIM.C b/src/SIM/GenAlphaSIM.C
new file mode 100644
index 00000000..c2a32d24
--- /dev/null
+++ b/src/SIM/GenAlphaSIM.C
@@ -0,0 +1,189 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GenAlphaSIM.C
+//!
+//! \date Aug 21 2014
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Generalized-alpha driver for isogeometric dynamic FEM simulators.
+//!
+//==============================================================================
+
+#include "GenAlphaSIM.h"
+#include "SIMoutput.h"
+#include "TimeStep.h"
+#include "IFEM.h"
+#include "tinyxml.h"
+
+
+GenAlphaSIM::GenAlphaSIM (SIMbase& s) : NewmarkSIM(s), prevSol(3), tempSol(3)
+{
+  // Default Newmark parameters (alpha = -0.1)
+  alpha_m = 1.0;
+  alpha_f = 0.9;
+  beta = 0.3025;
+  gamma = 0.6;
+}
+
+
+bool GenAlphaSIM::parse (const TiXmlElement* elem)
+{
+  bool ok = this->NewmarkSIM::parse(elem);
+
+  if (!strcasecmp(elem->Value(),"newmarksolver"))
+  {
+    double alpha = -0.1;
+    const char* attr = elem->Attribute("alpha");
+    if (attr)
+    {
+      alpha = atof(attr);
+      alpha_f = alpha + 1.0;
+      alpha_m = 1.0;
+    }
+    else if ((attr = elem->Attribute("rho_inf")))
+    {
+      double rho = atof(attr);
+      alpha_f = 1.0/(1.0+rho);
+      alpha_m = (2.0-rho)/(1.0+rho);
+      alpha = alpha_f - alpha_m;
+    }
+    beta = 0.25*(1.0-alpha)*(1.0-alpha);
+    gamma = 0.5 - alpha;
+  }
+
+  return ok;
+}
+
+
+void GenAlphaSIM::initPrm ()
+{
+  model.setIntegrationPrm(0,alpha1);
+  model.setIntegrationPrm(1,fabs(alpha2));
+  model.setIntegrationPrm(2,alpha_m);
+  model.setIntegrationPrm(3,alpha_f);
+  model.setIntegrationPrm(4,2.0);
+}
+
+
+bool GenAlphaSIM::initSol (size_t nSol)
+{
+  this->MultiStepSIM::initSol(nSol);
+
+  size_t nDOFs = model.getNoDOFs();
+  for (size_t i = 0; i < 3; i++)
+  {
+    prevSol[i].resize(nDOFs,true);
+    tempSol[i].resize(nDOFs,true);
+  }
+
+  return true;
+}
+
+
+bool GenAlphaSIM::advanceStep (TimeStep& param, bool updateTime)
+{
+  // Update displacement solutions between time steps
+  for (int n = solution.size()-3; n > 0; n--)
+    std::copy(solution[n-1].begin(),solution[n-1].end(),solution[n].begin());
+
+  // Update the previous solution
+  std::copy(solution.front().begin(),solution.front().end(),prevSol[0].begin());
+  if (solution.size() > 2)
+  {
+    size_t iA = solution.size() - 1;
+    size_t iV = solution.size() - 2;
+    std::copy(solution[iV].begin(),solution[iV].end(),prevSol[1].begin());
+    std::copy(solution[iA].begin(),solution[iA].end(),prevSol[2].begin());
+  }
+
+  return this->NewmarkSIM::advanceStep(param,updateTime);
+}
+
+
+bool GenAlphaSIM::predictStep (TimeStep& param)
+{
+  if (!this->NewmarkSIM::predictStep(param))
+    return false;
+
+  size_t ipD = 0;
+  size_t ipA = solution.size() - 1;
+  size_t ipV = solution.size() - 2;
+  tempSol[0] = solution[ipD];
+  tempSol[1] = solution[ipV];
+  tempSol[2] = solution[ipA];
+
+  // Compute the intermediate solution
+  // {U}_(n+alpha) = (1-alpha)*{U}_n + alpha*{U}_(n+1)
+  solution[ipD].relax(alpha_f,prevSol[0]);
+  solution[ipV].relax(alpha_f,prevSol[1]);
+  solution[ipA].relax(alpha_m,prevSol[2]);
+  return true;
+}
+
+
+bool GenAlphaSIM::correctStep (TimeStep& param, bool converged)
+{
+  if (solution.size() < 3)
+  {
+    std::cerr <<" *** GenAlphaSIM::correctStep: Too few solution vectors "
+              << solution.size() << std::endl;
+    return false;
+  }
+
+#ifdef SP_DEBUG
+  std::cout <<"\nGenAlphaSIM::correctStep(converged="
+            << std::boolalpha << converged <<")";
+#endif
+
+  size_t ipD = 0;
+  size_t ipA = solution.size() - 1;
+  size_t ipV = solution.size() - 2;
+
+  // Update current displacement, velocity and acceleration solutions
+  tempSol[0].add(linsol,beta*param.time.dt*param.time.dt);
+  tempSol[1].add(linsol,gamma*param.time.dt);
+  tempSol[2].add(linsol,1.0);
+
+  if (converged)
+  {
+    // Converged solution
+    solution[ipD] = tempSol[0];
+    solution[ipV] = tempSol[1];
+    solution[ipA] = tempSol[2];
+  }
+  else
+  {
+    // Compute new intermediate solution
+    // {U}_(n+alpha) = (1-alpha)*{U}_n + alpha*{U}_(n+1)
+    solution[ipD].relax(alpha_f,prevSol[0],tempSol[0]);
+    solution[ipV].relax(alpha_f,prevSol[1],tempSol[1]);
+    solution[ipA].relax(alpha_m,prevSol[2],tempSol[2]);
+  }
+
+#if SP_DEBUG > 1
+  std::cout <<"\nCorrected displacement:"<< solution[ipD]
+            <<"Corrected velocity:"<< solution[ipV]
+            <<"Corrected acceleration:"<< solution[ipA];
+#elif defined(SP_DEBUG)
+  if (converged)
+    std::cout <<"\nConverged displacement:"<< solution[ipD].max()
+              <<"\nConverged velocity:"<< solution[ipV].max()
+              <<"\nConverged acceleration:"<< solution[ipA].max() << std::endl;
+#endif
+
+  model.updateRotations(linsol,alpha_f); //TODO look at this
+  return model.updateConfiguration(tempSol[0]);
+}
+
+
+void GenAlphaSIM::setSolution (const Vector& newSol, int idx)
+{
+  solution[idx] = newSol;
+
+  if (idx == 0)
+    tempSol[0] = newSol;
+  else if (solution.size() > 2)
+    tempSol[solution.size()-3+idx] = newSol;
+}
diff --git a/src/SIM/GenAlphaSIM.h b/src/SIM/GenAlphaSIM.h
new file mode 100644
index 00000000..15f5bba3
--- /dev/null
+++ b/src/SIM/GenAlphaSIM.h
@@ -0,0 +1,62 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GenAlphaSIM.h
+//!
+//! \date Aug 21 2014
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Generalized-alpha driver for isogeometric dynamic FEM simulators.
+//!
+//==============================================================================
+
+#ifndef _GEN_ALPHA_SIM_H
+#define _GEN_ALPHA_SIM_H
+
+#include "NewmarkSIM.h"
+
+
+/*!
+  \brief Generalized-alpha driver for dynamic isogeometric FEM simulators.
+*/
+
+class GenAlphaSIM : public NewmarkSIM
+{
+public:
+  //! \brief The constructor initializes default solution parameters.
+  GenAlphaSIM(SIMbase& s);
+  //! \brief Empty destructor.
+  virtual ~GenAlphaSIM() {}
+
+  //! \brief Parses a data section from an XML document.
+  virtual bool parse(const TiXmlElement* elem);
+
+  //! \brief Initializes time integration parameters for the integrand.
+  virtual void initPrm();
+  //! \brief Initializes primary solution vectors.
+  //! \param[in] nSol Number of consequtive solutions stored
+  virtual bool initSol(size_t nSol = 3);
+
+  //! \brief Advances the time step one step forward.
+  //! \param param Time stepping parameters
+  //! \param[in] updateTime If \e false, the time parameters are not incremented
+  virtual bool advanceStep(TimeStep& param, bool updateTime = true);
+
+  //! \brief Modifies the current solution vector (used by sub-iterations only).
+  virtual void setSolution(const Vector& newSol, int idx);
+
+protected:
+  //! \brief Calculates predicted velocities and accelerations.
+  virtual bool predictStep(TimeStep& param);
+  //! \brief Updates configuration variables (solution vector) in an iteration.
+  virtual bool correctStep(TimeStep& param, bool converged);
+
+private:
+  double alpha_m;  //!< Generalized-alpha parameter
+  double alpha_f;  //!< Generalized-alpha parameter
+  Vectors prevSol; //!< Solution of last converged time step
+  Vectors tempSol; //!< Intermediate solution
+};
+
+#endif