diff --git a/examples/spectralLinearElastic2Materials.i b/examples/spectralLinearElastic2Materials.i
new file mode 100644
index 00000000..7641a7e8
--- /dev/null
+++ b/examples/spectralLinearElastic2Materials.i
@@ -0,0 +1,142 @@
+[Mesh]
+  type = MyTRIMMesh
+  dim = 3
+  xmax = 10
+  ymax = 10
+  zmax = 10
+  nx = 64
+  ny = 64
+  nz = 64
+[]
+
+[Problem]
+  type = FFTProblem
+[]
+
+[AuxVariables]
+  [./epsilon_aux_var]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+
+  [./stress_aux_var]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+
+  [./elastic_aux_var]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+
+  [./index_buffer_aux_var]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+
+  [./stiffness_ratio_aux]
+    order = CONSTANT
+    family = MONOMIAL
+    [./InitialCondition]
+      type = SmoothSuperellipsoidIC
+      x1 = 5
+      y1 = 5
+      z1 = 5
+      a = 2
+      b = 2
+      c = 2
+      n = 2
+      int_width = 0
+      invalue = 1
+      outvalue = 2.0
+    [../]
+  [../]
+[]
+
+[UserObjects]
+  # Buffers
+  [./epsilon]
+    type = RankTwoTensorFFTWBuffer
+  [../]
+  [./stress]
+    type = RankTwoTensorFFTWBuffer
+  [../]
+  # Reciprocal space: Elastic tensor
+  [./gamma]
+    type = RankFourTensorFFTWBuffer
+  [../]
+  [./elastic]
+    type = RankFourTensorFFTWBuffer
+  [../]
+  [./stiffness_ratio]
+     type = RealFFTWBuffer
+     moose_variable = stiffness_ratio_aux
+  []
+  [./index_buffer]
+     type = RealFFTWBuffer
+  []
+[]
+
+[AuxKernels]
+  [./epsilon_aux]
+    type = FFTBufferAux
+    variable = epsilon_aux_var
+    fft_buffer = epsilon
+    execute_on = final
+    component = '0 1'
+  [../]
+
+  [./stress_aux]
+    type = FFTBufferAux
+    variable = stress_aux_var
+    fft_buffer = stress
+    execute_on = final
+    component = '0 1'
+  [../]
+
+  [./stiffness_aux]
+    type = FFTBufferAux
+    variable = stiffness_ratio_aux
+    fft_buffer = stiffness_ratio
+    execute_on = final
+  [../]
+
+  [./index_aux]
+    type = FFTBufferAux
+    variable = index_buffer_aux_var
+    fft_buffer = index_buffer
+    execute_on = final
+  [../]
+[]
+
+[Executioner]
+  type = SpectralExecutionerLinearElastic
+
+  time_step = 1.0
+  number_steps = 150
+  initial_shear_strain = 0.01
+  young_modulus = 1e4
+  poisson_ratio = 0.3
+
+[]
+
+[VectorPostprocessors]
+  [./linevalue]
+    type = LineValueSampler
+    variable = 'stress_aux_var'
+    start_point = '0 0 0'
+    end_point = '9.9999999999 0 0'
+    num_points = 101
+    sort_by = x
+    execute_on = final
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+  execute_on = 'INITIAL FINAL'
+  perf_graph = true
+  [./comp]
+    type = CSV
+  [../]
+[]
diff --git a/include/executioners/SpectralExecutionerLinearElastic.h b/include/executioners/SpectralExecutionerLinearElastic.h
new file mode 100644
index 00000000..ef01074e
--- /dev/null
+++ b/include/executioners/SpectralExecutionerLinearElastic.h
@@ -0,0 +1,88 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#pragma once
+
+#include "SpectralExecutionerBase.h"
+
+// System includes
+#include <string>
+
+// Forward declarations
+class InputParameters;
+
+/**
+ * Executioner for diffusion spectral solver.
+ */
+class SpectralExecutionerLinearElastic : public SpectralExecutionerBase
+{
+public:
+  static InputParameters validParams();
+
+  SpectralExecutionerLinearElastic(const InputParameters & parameters);
+  /**
+   * Algorithm for incremental solution using forward/backward transforms of Green's function.
+   */
+  virtual void execute() final;
+
+protected:
+  /// Time step
+  const Real _dt;
+
+  /// Number of steps
+  const unsigned int _nsteps;
+
+  /// First parameter
+  const Real _young_modulus;
+
+  /// Second parameter
+  const Real _poisson_ratio;
+
+  /// Current time
+  Real _t_current;
+
+  /// Initial homogeneous shear deformation
+  const Real _initial_shear_strain;
+
+  /// Initial strain tensor
+  RankTwoTensor _initial_strain_tensor;
+
+private:
+  /**
+   * Helper function to get the diffusion equation Green's function corresponding to one time step.
+   */
+  void getGreensFunction(FFTBufferBase<RankFourTensor> & gamma_hat,
+                         FFTBufferBase<Real> & ratio_buffer);
+
+  /**
+   * Helper function to get the initial stress from strain and tensor of elastic coefficients.
+   */
+  FFTBufferBase<RankTwoTensor> & getInitialStress(FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+                                                  FFTBufferBase<RankFourTensor> & elastic_buffer);
+  /**
+   * Helper to populate initial strain buffer
+   */
+  void populateEpsilonBuffer(FFTBufferBase<RankTwoTensor> & epsilon_buffer);
+
+  /**
+   * Advance Fourier epsilon to next iteration
+   */
+  void advanceReciprocalEpsilon(FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+                                FFTBufferBase<RankTwoTensor> & stress_buffer,
+                                const FFTBufferBase<RankFourTensor> & gamma_hat);
+  /**
+   * Update sigma in real space for this iteration
+   */
+  void updateRealSigma(FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+                       FFTBufferBase<RankTwoTensor> & stress_buffer,
+                       FFTBufferBase<RankFourTensor> & elastic_tensor);
+
+  void filloutElasticTensor(const FFTBufferBase<Real> & ratio_buffer,
+                            FFTBufferBase<Real> & index_buffer,
+                            FFTBufferBase<RankFourTensor> & elastic_tensor_buffer);
+};
diff --git a/src/executioners/SpectralExecutionerLinearElastic.C b/src/executioners/SpectralExecutionerLinearElastic.C
new file mode 100644
index 00000000..3a7ae926
--- /dev/null
+++ b/src/executioners/SpectralExecutionerLinearElastic.C
@@ -0,0 +1,325 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "SpectralExecutionerLinearElastic.h"
+#include "InputParameters.h"
+#include "RankFourTensor.h"
+
+registerMooseObject("MagpieApp", SpectralExecutionerLinearElastic);
+
+Real
+deltaij(const int i, const int j)
+{
+  if (i == j)
+    return 1.0;
+  else
+    return 0.0;
+}
+
+InputParameters
+SpectralExecutionerLinearElastic::validParams()
+{
+  InputParameters params = SpectralExecutionerBase::validParams();
+  params.addClassDescription("Executioner for spectral solve of diffusion equation");
+  params.addParam<Real>("time_step", 1.0, "Time step for ODE integration");
+  params.addParam<unsigned int>("number_steps", 0.0, "Time step for ODE integration");
+  params.addParam<Real>("young_modulus", 1.0, "First parameter for isotropic materials");
+  params.addParam<Real>("poisson_ratio", 1.0, "Second parameter for isotropic materials");
+  params.addParam<Real>(
+      "initial_shear_strain", 0.001, "Homogeneous two-dimensional shear deformation");
+
+  return params;
+}
+
+SpectralExecutionerLinearElastic::SpectralExecutionerLinearElastic(
+    const InputParameters & parameters)
+  : SpectralExecutionerBase(parameters),
+    _dt(getParam<Real>("time_step")),
+    _nsteps(getParam<unsigned int>("number_steps")),
+    _young_modulus(getParam<Real>("young_modulus")),
+    _poisson_ratio(getParam<Real>("poisson_ratio")),
+    _initial_shear_strain(getParam<Real>("initial_shear_strain"))
+{
+  // Add check that's a 2D problem with LIBMESH_DIM == 2
+  _initial_strain_tensor(0, 0) = 0.0;
+  _initial_strain_tensor(0, 1) = _initial_strain_tensor(1, 0) = _initial_shear_strain;
+  _initial_strain_tensor(2, 1) = _initial_strain_tensor(1, 2) = 0.0;
+  _initial_strain_tensor(0, 2) = _initial_strain_tensor(2, 0) = 0.0;
+  _initial_strain_tensor(1, 1) = _initial_strain_tensor(2, 2) = 0.0;
+
+  _t_current = 0.0;
+}
+
+void
+SpectralExecutionerLinearElastic::populateEpsilonBuffer(
+    FFTBufferBase<RankTwoTensor> & epsilon_buffer)
+{
+  const auto & grid = epsilon_buffer.grid();
+  const int ni = grid[0];
+  const int nj = grid[1]; // Real space.
+  const int nk = grid[2]; // Real space.
+
+  std::size_t index = 0;
+
+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+  {
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+    {
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        epsilon_buffer.realSpace()[index] = _initial_strain_tensor;
+        index++;
+      }
+    }
+  }
+}
+
+void
+SpectralExecutionerLinearElastic::getGreensFunction(FFTBufferBase<RankFourTensor> & gamma_hat,
+                                                    FFTBufferBase<Real> & ratio_buffer)
+{
+  const auto & grid = gamma_hat.grid();
+  const int ndim = 3;
+
+  const int ni = grid[0];
+  const int nj = grid[1];            // Real space.
+  const int nk = (grid[2] >> 1) + 1; // Real space.
+
+  std::size_t index = 0;
+
+  // To plug the right frequencies
+  const auto & ivec = gamma_hat.kTable(0);
+  const auto & jvec = gamma_hat.kTable(1);
+  const auto & kvec = gamma_hat.kTable(2);
+
+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        const std::array<Real, 3> freq{ivec[freq_x], jvec[freq_y], kvec[freq_z]};
+
+        Real lambda = _young_modulus * ratio_buffer.realSpace()[index] * _poisson_ratio / ((1 + _poisson_ratio) * (1 - 2 * _poisson_ratio));
+        Real nu = _young_modulus * ratio_buffer.realSpace()[index] / (2 * (1 + _poisson_ratio));
+
+        for (int i = 0; i < ndim; i++)
+          for (int j = 0; j < ndim; j++)
+            for (int k = 0; k < ndim; k++)
+              for (int l = 0; l < ndim; l++)
+              {
+                Real q_square =
+                    std::pow(freq[0] * freq[0] + freq[1] * freq[1] + freq[2] * freq[2], 2.0);
+                gamma_hat.reciprocalSpace()[index](i, j, k, l) =
+                    -1.0 * ((lambda + nu) / (nu * (lambda + 2.0 * nu))) *
+                        (freq[i] * freq[j] * freq[k] * freq[l]) / (q_square * q_square) +
+                    (deltaij(j, k) * freq[i] * freq[l] + deltaij(j, l) * freq[i] * freq[k] +
+                     deltaij(i, k) * freq[j] * freq[l] + deltaij(i, l) * freq[j] * freq[k]) /
+                        (4.0 * nu * q_square);
+              }
+        // index refers to space
+        index++;
+      }
+}
+
+FFTBufferBase<RankTwoTensor> &
+SpectralExecutionerLinearElastic::getInitialStress(FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+                                                   FFTBufferBase<RankFourTensor> & elastic_tensor)
+{
+  auto & stress_buffer = getFFTBuffer<RankTwoTensor>("stress");
+  const auto & grid = epsilon_buffer.grid();
+  int ni = grid[0];
+  int nj = grid[1];
+  int nk = grid[2];
+
+  // Homogeneous initial state of strain
+  populateEpsilonBuffer(epsilon_buffer);
+
+  size_t index = 0;
+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        stress_buffer.realSpace()[index] =
+            elastic_tensor.realSpace()[index] * epsilon_buffer.realSpace()[index];
+        index++;
+      }
+
+  return stress_buffer;
+}
+
+void
+SpectralExecutionerLinearElastic::advanceReciprocalEpsilon(
+    FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+    FFTBufferBase<RankTwoTensor> & stress_buffer,
+    const FFTBufferBase<RankFourTensor> & gamma_hat)
+{
+  const auto & grid = epsilon_buffer.grid();
+  int ni = grid[0];
+  int nj = grid[1];
+  int nk = (grid[2] >> 1) + 1;
+
+  size_t index = 0;
+
+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        if (freq_x == 0 && freq_y == 0 && freq_z == 0)
+        {
+          epsilon_buffer.reciprocalSpace()[index] = _initial_strain_tensor;
+        }
+        else
+        {
+          epsilon_buffer.reciprocalSpace()[index] -=
+              gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index];
+        }
+        index++;
+      }
+}
+
+void
+SpectralExecutionerLinearElastic::updateRealSigma(FFTBufferBase<RankTwoTensor> & epsilon_buffer,
+                                                  FFTBufferBase<RankTwoTensor> & stress_buffer,
+                                                  FFTBufferBase<RankFourTensor> & elastic_tensor)
+{
+  const auto & grid = epsilon_buffer.grid();
+  int ni = grid[0];
+  int nj = grid[1];
+  int nk = grid[2];
+
+  size_t index = 0;
+
+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        stress_buffer.realSpace()[index] =
+            elastic_tensor.realSpace()[index] * epsilon_buffer.realSpace()[index];
+        index++;
+      }
+}
+
+
+void
+SpectralExecutionerLinearElastic::filloutElasticTensor(const FFTBufferBase<Real> & ratio_buffer,
+                                                  FFTBufferBase<Real> & index_buffer,
+                                                  FFTBufferBase<RankFourTensor> & elastic_tensor_buffer) {
+const auto & grid = elastic_tensor_buffer.grid();
+
+int ni = grid[0];
+int nj = grid[1];
+int nk = grid[2];
+
+size_t index = 0;
+
+for (int freq_x = 0; freq_x < ni; ++freq_x)
+  for (int freq_y = 0; freq_y < nj; ++freq_y)
+    for (int freq_z = 0; freq_z < nk; ++freq_z)
+    {
+      // Define elastic tensor
+      std::vector<Real> iso_const(2);
+      iso_const[0] = _young_modulus * ratio_buffer.realSpace()[index];
+      iso_const[1] = _poisson_ratio;
+
+      elastic_tensor_buffer.realSpace()[index].fillFromInputVector(
+          iso_const, RankFourTensor::symmetric_isotropic_E_nu);
+      index_buffer.realSpace()[index] = index;
+      index++;
+    }
+}
+
+void
+SpectralExecutionerLinearElastic::execute()
+{
+  unsigned int thisStep = 0;
+
+  _time_step = thisStep;
+  _time = _time_step;
+  _fe_problem.outputStep(EXEC_INITIAL);
+  _fe_problem.advanceState();
+
+  auto & epsilon_buffer = getFFTBuffer<RankTwoTensor>("epsilon");
+  populateEpsilonBuffer(epsilon_buffer);
+
+  /*  --------------------------------------------------- */
+  auto & ratio_buffer = getFFTBuffer<Real>("stiffness_ratio");
+  auto & elastic_tensor_buffer = getFFTBuffer<RankFourTensor>("elastic");
+  auto & index_buffer = getFFTBuffer<Real>("index_buffer");
+
+  filloutElasticTensor(ratio_buffer, index_buffer, elastic_tensor_buffer);
+  /*  --------------------------------------------------- */
+  //  Moose::out << "epsilon_buffer: " << epsilon_buffer.realSpace()[0] << "\n";
+  epsilon_buffer.realSpace()[0].print();
+  // Get corresponding initial stress
+  auto & stress_buffer = getInitialStress(epsilon_buffer, elastic_tensor_buffer);
+
+  //  Moose::out << "stress_buffer: " << stress_buffer.realSpace()[0] << "\n";
+  stress_buffer.realSpace()[1].print();
+  // Get specific Green's function
+  auto & gamma_hat = getFFTBuffer<RankFourTensor>("gamma");
+
+  thisStep++;
+  _t_current += _dt;
+  _time_step = thisStep;
+
+  _fe_problem.execute(EXEC_FINAL);
+  _time = _t_current;
+
+  Moose::out << "_t_current: " << _t_current << ". \n";
+
+  _fe_problem.outputStep(EXEC_FINAL);
+
+  if (gamma_hat.dim() != 3)
+    mooseError("Error: Dimension not implemented in SpectralExecutionerLinearElastic.");
+
+  getGreensFunction(gamma_hat, ratio_buffer);
+
+  // Our plans do not preserve the inputs (unfortunately)
+  FFTData<RankTwoTensor> stress_buffer_backup = stress_buffer.realSpace();
+  FFTData<ComplexType<RankTwoTensor>::type> epsilon_buffer_backup = stress_buffer.reciprocalSpace();
+  Real scale_factor = stress_buffer.backwardScale();
+
+  for (unsigned int step_no = 0; step_no < _nsteps; step_no++)
+  {
+    // (a)
+	stress_buffer_backup = stress_buffer.realSpace();
+    stress_buffer.forward();
+
+    // (b)
+    // Convergence test
+
+    // Compute new strain tensor in Fourier space
+    // (c)
+    // Preserve data
+    epsilon_buffer.reciprocalSpace() = epsilon_buffer_backup;
+    advanceReciprocalEpsilon(epsilon_buffer, stress_buffer, gamma_hat);
+
+    // (d)
+    epsilon_buffer_backup = epsilon_buffer.reciprocalSpace();
+    epsilon_buffer.backward();
+
+    // (e)
+    // Preserve data
+    stress_buffer.realSpace() = stress_buffer_backup;
+    updateRealSigma(epsilon_buffer, stress_buffer, elastic_tensor_buffer);
+    // End of fixed-point iterations
+
+    thisStep++;
+    _t_current += _dt;
+    _time_step = thisStep;
+
+    _fe_problem.execute(EXEC_FINAL);
+    _time = _t_current;
+
+    Moose::out << "_t_current: " << _t_current << ". \n";
+
+    _fe_problem.outputStep(EXEC_FINAL);
+
+    if (step_no != _nsteps - 1)
+      _fe_problem.advanceState();
+  }
+}
diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C
index 33a16e4d..c9e86096 100644
--- a/src/userobjects/FFTBufferBase.C
+++ b/src/userobjects/FFTBufferBase.C
@@ -162,6 +162,7 @@ void
 FFTBufferBase<T>::forward()
 {
   forwardRaw();
+  // _reciprocal_space_data *= std::sqrt(backwardScale());
 }
 
 template <typename T>
diff --git a/src/userobjects/FFTWBufferBase.C b/src/userobjects/FFTWBufferBase.C
index 7ed4cb60..796562e9 100644
--- a/src/userobjects/FFTWBufferBase.C
+++ b/src/userobjects/FFTWBufferBase.C
@@ -17,13 +17,15 @@
 template <typename T>
 FFTWBufferBase<T>::FFTWBufferBase(const InputParameters & parameters)
   : FFTBufferBase<T>(parameters),
-    _perf_plan(this->registerTimedSection("fftw_plan_r2r", 2)),
+    _perf_plan(this->registerTimedSection("fftw_plan_r2c", 2)),
     _perf_fft(this->registerTimedSection("fftw_execute", 2))
 {
   // create plans
   {
     TIME_SECTION(_perf_plan);
 
+    // We should think about preserving input here
+    // Option not available for generic "many call"
     _forward_plan =
         fftw_plan_many_dft_r2c(_dim,
                                _grid.data(),