diff --git a/.gitignore b/.gitignore
index 7076a398..67e5c49e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -52,7 +52,6 @@ temp_print_trace.*
 *.poly
 *.mpx
 *.btr
-*.xml
 *.gmv
 *.plt
 *.slh
diff --git a/Makefile b/Makefile
index 32933054..8d85a071 100644
--- a/Makefile
+++ b/Makefile
@@ -45,6 +45,13 @@ POROUS_FLOW         := no
 include $(MOOSE_DIR)/modules/modules.mk
 ###############################################################################
 
+ifdef NEML
+	ADDITIONAL_INCLUDES := -I$(NEML)/src
+	ADDITIONAL_LIBS 	  := -L$(NEML)/lib -lneml
+endif
+
+##############################################################################
+
 # dep apps
 APPLICATION_DIR    := $(CURDIR)
 APPLICATION_NAME   := blackbear
diff --git a/include/materials/ComputeNEMLStress.h b/include/materials/ComputeNEMLStress.h
new file mode 100644
index 00000000..59ec8053
--- /dev/null
+++ b/include/materials/ComputeNEMLStress.h
@@ -0,0 +1,52 @@
+#ifndef COMPUTENEMLSTRESS_H
+#define COMPUTENEMLSTRESS_H
+
+#include "ComputeStressBase.h"
+
+#include "neml_interface.h"
+
+class ComputeNEMLStress;
+
+template <>
+InputParameters validParams<ComputeNEMLStress>();
+
+class ComputeNEMLStress: public ComputeStressBase
+{
+ public:
+  ComputeNEMLStress(const InputParameters & parameters);
+
+  virtual void computeQpStress() override;
+  virtual void initQpStatefulProperties() override;
+  bool isElasticityTensorGuaranteedIsotropic();
+
+ protected:
+  FileName _fname;
+  std::string _mname;
+  std::unique_ptr<neml::NEMLModel> _model;
+  MaterialProperty<std::vector<Real>> & _hist;
+  const MaterialProperty<std::vector<Real>> & _hist_old;
+  const MaterialProperty<RankTwoTensor> & _mechanical_strain_old;
+  const MaterialProperty<RankTwoTensor> & _stress_old;
+  MaterialProperty<Real> & _energy;
+  const MaterialProperty<Real> & _energy_old;
+  MaterialProperty<Real> & _dissipation;
+  const MaterialProperty<Real> & _dissipation_old;
+  const VariableValue & _temperature; // Will default to zero
+  const VariableValue & _temperature_old;
+  MaterialProperty<RankTwoTensor> & _inelastic_strain;
+  MaterialProperty<Real> & _shear_modulus;
+  MaterialProperty<Real> & _bulk_modulus;
+  MaterialProperty<RankFourTensor> & _elasticity_tensor;
+};
+
+/// Tensor -> my notation
+void tensor_neml(const RankTwoTensor & in, double * const out);
+
+/// Vector -> tensor
+void neml_tensor(const double * const in, RankTwoTensor & out);
+
+/// Tangent -> tensor
+void neml_tangent(const double * const in, RankFourTensor & out);
+
+
+#endif // COMPUTENEMLSTRESS_H
diff --git a/include/materials/ComputeThermalExpansionEigenstrainNEML.h b/include/materials/ComputeThermalExpansionEigenstrainNEML.h
new file mode 100644
index 00000000..3d76259a
--- /dev/null
+++ b/include/materials/ComputeThermalExpansionEigenstrainNEML.h
@@ -0,0 +1,48 @@
+#ifndef COMPUTETHERMALEXPANSIONEIGENSTRAINNEML_H
+#define COMPUTETHERMALEXPANSIONEIGENSTRAINNEML_H
+
+#include "ComputeThermalExpansionEigenstrainBase.h"
+
+#include "neml_interface.h"
+
+class ComputeThermalExpansionEigenstrainNEML;
+
+template <>
+InputParameters validParams<ComputeThermalExpansionEigenstrainNEML>();
+
+/**
+ *  ComputeThermalExpansionEigenstrainNEML computes the thermal expansion
+ *  strain from the instantaneous CTE provided by a NEML model
+ */
+class ComputeThermalExpansionEigenstrainNEML: public ComputeThermalExpansionEigenstrainBase
+{
+ public:
+  ComputeThermalExpansionEigenstrainNEML(const InputParameters & parameters);
+  virtual void initQpStatefulProperties() override;
+
+ protected:
+  virtual void computeThermalStrain(Real & thermal_strain, Real & instantaneous_cte) override;
+
+ protected:
+  FileName _fname;
+  std::string _mname;
+  std::unique_ptr<neml::NEMLModel> _model;
+
+  MaterialProperty<Real> & _tstrain;
+  const MaterialProperty<Real> & _tstrain_old;
+
+  const VariableValue & _temperature_old;
+};
+
+
+
+
+
+
+
+
+
+
+
+
+#endif // COMPUTETHERMALEXPANSIONEIGENSTRAINNEML_H
diff --git a/src/materials/ComputeNEMLStress.C b/src/materials/ComputeNEMLStress.C
new file mode 100644
index 00000000..ca1e5f3f
--- /dev/null
+++ b/src/materials/ComputeNEMLStress.C
@@ -0,0 +1,187 @@
+#include "ComputeNEMLStress.h"
+
+registerMooseObject("BlackBearApp", ComputeNEMLStress);
+
+template <>
+InputParameters
+validParams<ComputeNEMLStress>()
+{
+  InputParameters params = validParams<ComputeStressBase>();
+  params.addRequiredParam<FileName>("database", "Path to NEML XML database.");
+  params.addRequiredParam<std::string>("model", "Model name in NEML database.");
+  params.addCoupledVar("temperature", 0.0, "Coupled temperature");
+  return params;
+}
+
+ComputeNEMLStress::ComputeNEMLStress(const InputParameters & parameters) :
+    ComputeStressBase(parameters),
+    _fname(getParam<FileName>("database")),
+    _mname(getParam<std::string>("model")),
+    _hist(declareProperty<std::vector<Real>>(_base_name + "hist")),
+    _hist_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "hist")),
+    _mechanical_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "mechanical_strain")),
+    _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
+    _energy(declareProperty<Real>(_base_name + "energy")),
+    _energy_old(getMaterialPropertyOld<Real>(_base_name + "energy")),
+    _dissipation(declareProperty<Real>(_base_name + "dissipation")),
+    _dissipation_old(getMaterialPropertyOld<Real>(_base_name + "dissipation")),
+    _temperature(coupledValue("temperature")),
+    _temperature_old(coupledValueOld("temperature")),
+    _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "inelastic_strain")),
+    _shear_modulus(declareProperty<Real>(_base_name + "shear_modulus")),
+    _bulk_modulus(declareProperty<Real>(_base_name + "bulk_modulus")),
+    _elasticity_tensor(declareProperty<RankFourTensor>(_base_name + "elasticity_tensor"))
+{
+  // I strongly hesitate to put this here, may change later
+  _model = neml::parse_xml_unique(_fname, _mname);
+}
+
+void ComputeNEMLStress::computeQpStress()
+{
+  // We must update:
+  // 1) _stress
+  // 2) _Jacobian_mult
+  // 3) _elastic_strain
+  // 4) _history
+
+  // First do some declaration and translation
+  double s_np1[6];
+  double s_n[6];
+  tensor_neml(_stress_old[_qp], s_n);
+
+  double e_np1[6];
+  tensor_neml(_mechanical_strain[_qp], e_np1);
+  double e_n[6];
+  tensor_neml(_mechanical_strain_old[_qp], e_n);
+
+  double t_np1 = _t;
+  double t_n = _t - _dt;
+
+  double T_np1 = _temperature[_qp];
+  double T_n = _temperature_old[_qp];
+  
+  double * h_np1 = &(_hist[_qp][0]);
+  const double * const h_n = &(_hist_old[_qp][0]);
+
+  double A_np1[36];
+  
+  double u_np1;
+  double u_n = _energy_old[_qp];
+
+  double p_np1;
+  double p_n = _dissipation_old[_qp];
+
+  double estrain[6];
+
+  int ier;
+
+  // Actually call the update
+  ier = _model->update_sd(e_np1, e_n, T_np1, T_n, t_np1, t_n,
+                    s_np1, s_n, h_np1, h_n, A_np1, u_np1, u_n,
+                    p_np1, p_n);
+
+  if (ier != neml::SUCCESS)
+    throw MooseException("NEML stress update failed!");
+
+  // Do more translation, now back to tensors
+  neml_tensor(s_np1, _stress[_qp]);
+  neml_tangent(A_np1, _Jacobian_mult[_qp]);
+
+  // Get the elastic strain
+  ier = _model->elastic_strains(s_np1, T_np1, h_np1, estrain);
+
+  if (ier != neml::SUCCESS)
+    mooseError("Error in NEML call for elastic strain!");
+
+  // Translate
+  neml_tensor(estrain, _elastic_strain[_qp]);
+
+  // For EPP purposes calculate the inelastic strain
+  double pstrain[6];
+  for (int i=0; i<6; i++) {
+    pstrain[i] = e_np1[i] - estrain[i];
+  }
+  neml_tensor(pstrain, _inelastic_strain[_qp]);
+
+  // Store dissipation
+  _energy[_qp] = u_np1;
+  _dissipation[_qp] = p_np1;
+
+  // Store elastic properties at current time
+  double mu = _model->shear(T_np1);
+  double K = _model->bulk(T_np1);
+  double l = K - 2.0 * mu / 3.0;
+  
+  _shear_modulus[_qp] = mu;
+  _bulk_modulus[_qp] = K;
+
+  std::vector<Real> props({l, mu});
+  _elasticity_tensor[_qp].fillFromInputVector(props, RankFourTensor::FillMethod::symmetric_isotropic);
+
+}
+
+void ComputeNEMLStress::initQpStatefulProperties()
+{
+  ComputeStressBase::initQpStatefulProperties();
+
+  // Figure out initial history
+  int ier;
+  _hist[_qp].resize(_model->nhist());
+  ier = _model->init_hist(&(_hist[_qp][0]));
+  
+
+  // Various other junk
+  _energy[_qp] = 0.0;
+  _dissipation[_qp] = 0.0;
+
+  if (ier != neml::SUCCESS)
+    mooseError("Error initializing NEML history!");
+}
+
+bool ComputeNEMLStress::isElasticityTensorGuaranteedIsotropic()
+{
+  // Technically we could query the model but it's just safer to say
+  // no as it's not in general for a NEML material model.
+  return false;
+}
+
+void tensor_neml(const RankTwoTensor & in, double * const out)
+{
+  double inds[6][2] = {{0,0}, {1,1}, {2,2}, {1,2}, {0,2}, {0,1}};
+  double mults[6] = {1.0, 1.0, 1.0, sqrt(2.0), sqrt(2.0), sqrt(2.0)};
+
+  for (int i = 0; i < 6; i++) {
+    out[i] = in(inds[i][0], inds[i][1]) * mults[i];
+  }
+}
+
+
+void neml_tensor(const double * const in, RankTwoTensor & out)
+{
+  double inds[6][2] = {{0,0}, {1,1}, {2,2}, {1,2}, {0,2}, {0,1}};
+  double mults[6] = {1.0, 1.0, 1.0, sqrt(2.0), sqrt(2.0), sqrt(2.0)};
+
+  for (int i = 0; i < 6; i++) {
+    out(inds[i][0], inds[i][1]) = in[i] / mults[i];
+    out(inds[i][1], inds[i][0]) = in[i] / mults[i];
+  }
+}
+
+void neml_tangent(const double * const in, RankFourTensor & out)
+{
+  double inds[6][2] = {{0,0}, {1,1}, {2,2}, {1,2}, {0,2}, {0,1}};
+  double mults[6] = {1.0, 1.0, 1.0, sqrt(2.0), sqrt(2.0), sqrt(2.0)};
+
+  for (int i = 0; i < 6; i++) {
+    for (int j = 0; j < 6; j++) {
+      out(inds[i][0], inds[i][1], inds[j][0], inds[j][1]) = in[i*6+j] / (
+          mults[i] * mults[j]);
+      out(inds[i][1], inds[i][0], inds[j][0], inds[j][1]) = in[i*6+j] / (
+          mults[i] * mults[j]);
+      out(inds[i][0], inds[i][1], inds[j][1], inds[j][0]) = in[i*6+j] / (
+          mults[i] * mults[j]);
+      out(inds[i][1], inds[i][0], inds[j][1], inds[j][0]) = in[i*6+j] / (
+          mults[i] * mults[j]);
+    }
+  }
+}
diff --git a/src/materials/ComputeThermalExpansionEigenstrainNEML.C b/src/materials/ComputeThermalExpansionEigenstrainNEML.C
new file mode 100644
index 00000000..d5649072
--- /dev/null
+++ b/src/materials/ComputeThermalExpansionEigenstrainNEML.C
@@ -0,0 +1,48 @@
+#include "ComputeThermalExpansionEigenstrainNEML.h"
+#include <string>
+
+registerMooseObject("BlackBearApp", ComputeThermalExpansionEigenstrainNEML);
+
+template <>
+InputParameters
+validParams<ComputeThermalExpansionEigenstrainNEML>()
+{
+  InputParameters params = validParams<ComputeThermalExpansionEigenstrainBase>();
+  params.addRequiredParam<FileName>("database", "Path to NEML XML database.");
+  params.addRequiredParam<std::string>("model", "Model name in NEML database.");
+  return params;
+}
+
+ComputeThermalExpansionEigenstrainNEML::ComputeThermalExpansionEigenstrainNEML(
+    const InputParameters & parameters) :
+    ComputeThermalExpansionEigenstrainBase(parameters),
+    _fname(getParam<FileName>("database")),
+    _mname(getParam<std::string>("model")),
+    _tstrain(declareProperty<Real>(_base_name + "tstrain")),
+    _tstrain_old(getMaterialPropertyOld<Real>(_base_name + "tstrain")),
+    _temperature_old(coupledValueOld("temperature"))
+{
+  // I strongly hesitate to put this here, may change later
+  _model = neml::parse_xml_unique(_fname, _mname);
+}
+
+void
+ComputeThermalExpansionEigenstrainNEML::computeThermalStrain(
+    Real & thermal_strain,
+    Real & instantaneous_cte)
+{
+  double nemlCTE = _model->alpha(_temperature[_qp]);
+  double nemlCTE_old = _model->alpha(_temperature_old[_qp]);
+
+  thermal_strain = _tstrain_old[_qp] + (nemlCTE + nemlCTE_old) / 2 * (
+      _temperature[_qp] - _temperature_old[_qp]);
+  
+  instantaneous_cte = nemlCTE;
+  _tstrain[_qp] = thermal_strain;
+}
+
+void ComputeThermalExpansionEigenstrainNEML::initQpStatefulProperties()
+{
+  ComputeThermalExpansionEigenstrainBase::initQpStatefulProperties();
+  _tstrain[_qp] = 0.0;
+}
diff --git a/test/tests/neml_complex/complex.xml b/test/tests/neml_complex/complex.xml
new file mode 100644
index 00000000..b64d2f1b
--- /dev/null
+++ b/test/tests/neml_complex/complex.xml
@@ -0,0 +1,31 @@
+<materials>
+ <example type="GeneralIntegrator">
+   <elastic type="IsotropicLinearElasticModel">
+     <m1>150000.0</m1>
+     <m1_type>youngs</m1_type>
+     <m2>0.3</m2>
+     <m2_type>poissons</m2_type>
+   </elastic>
+
+   <rule type="TVPFlowRule">
+     <elastic type="IsotropicLinearElasticModel">
+       <m1>150000.0</m1>
+       <m1_type>youngs</m1_type>
+       <m2>0.3</m2>
+       <m2_type>poissons</m2_type>
+     </elastic>
+
+     <flow type="PerzynaFlowRule">
+       <surface type="IsoJ2"/>
+       <hardening type="LinearIsotropicHardeningRule">
+         <s0>100.0</s0>
+         <K>2500.0</K>
+       </hardening>
+       <g type="GPowerLaw">
+         <n>5.0</n>
+         <eta>100.0</eta>
+       </g>
+     </flow>
+   </rule>
+ </example>
+</materials>
diff --git a/test/tests/neml_complex/gold/inelastic_out.csv b/test/tests/neml_complex/gold/inelastic_out.csv
new file mode 100644
index 00000000..aa8ca02f
--- /dev/null
+++ b/test/tests/neml_complex/gold/inelastic_out.csv
@@ -0,0 +1,102 @@
+time,stress_xx,stress_xy,stress_xz,stress_yy,stress_yz,stress_zz
+0,0,0,0,0,0,0
+100,110.42985068536,1.4174894412611e-15,1.4174894412611e-15,-1.4265818021405e-14,-1.0371873960447e-15,-7.5354630203067e-15
+200,115.73384612199,-1.400303161872e-15,-1.400303161872e-15,7.5788209648075e-09,-2.8653950688245e-15,7.5788190502008e-09
+300,118.24064980173,-6.5419115555319e-16,-6.5419115555319e-16,6.9735235224327e-09,5.3962378782659e-15,6.9735225129551e-09
+400,120.70046277077,-5.4208257476714e-15,-5.4208257476714e-15,6.9638758867792e-09,4.7526025552273e-15,6.9638838200393e-09
+500,123.15949244081,1.4167765725607e-15,1.4167765725607e-15,6.9544704917125e-09,-4.2969484568331e-15,6.9544801608946e-09
+600,125.61850905552,1.1854920882584e-15,1.1854920882584e-15,6.9637127169829e-09,7.58244002336e-15,6.963711161811e-09
+700,128.07752545265,8.3181677365844e-15,8.3181677365844e-15,6.9637594441538e-09,-6.488491450658e-15,6.9637648003911e-09
+800,130.53654184615,-5.3366662079378e-16,-5.3366662079378e-16,6.9637163894836e-09,6.710227182111e-15,6.9637171326369e-09
+900,132.9955582396,-8.5509751214222e-15,-8.5509751214222e-15,6.963738510209e-09,6.0844207906319e-16,6.9637360350571e-09
+1000,135.45457463304,1.827936817047e-14,1.827936817047e-14,6.9636914539423e-09,-1.4986872712827e-14,6.9636692937934e-09
+1100,137.91359102648,3.4341307864674e-15,3.4341307864674e-15,6.9637463712158e-09,-8.3513099200082e-15,6.9637357960684e-09
+1200,140.37260741992,8.0444348447578e-15,8.0444348447578e-15,6.9638322751182e-09,3.0123127124439e-14,6.9638330661885e-09
+1300,142.83162381337,-1.5368070540169e-14,-1.5368070540169e-14,6.9637451066126e-09,2.0616133546e-14,6.9637477800791e-09
+1400,145.29064020681,7.1133321836589e-15,7.1133321836589e-15,6.9636497280197e-09,1.8873192410983e-14,6.9636874754175e-09
+1500,147.74965660025,4.4204767582075e-15,4.4204767582075e-15,6.9636874713613e-09,7.4226789150016e-15,6.9637165077997e-09
+1600,150.20867299369,-6.2623415816541e-15,-6.2623415816541e-15,6.9637282831328e-09,-2.7001284955608e-14,6.9637489900135e-09
+1700,152.66768938714,-1.1785934019684e-14,-1.1785934019684e-14,6.9636617140929e-09,1.3863020085759e-14,6.9636787741439e-09
+1800,155.12670578058,-2.3453893151023e-15,-2.3453893151023e-15,6.9636778976569e-09,7.574052512022e-16,6.9636605742635e-09
+1900,157.58572217402,-1.8832362445588e-15,-1.8832362445588e-15,6.9638383918472e-09,-1.6034910352277e-14,6.9638272064937e-09
+2000,160.04473856747,-1.6883607459456e-16,-1.6883607459456e-16,6.9637053639138e-09,-1.8361673736551e-14,6.9637770107948e-09
+2100,162.50375496091,1.3177343490256e-14,1.3177343490256e-14,6.963894693666e-09,4.5964977810784e-15,6.9639302370994e-09
+2200,164.96277135435,-4.1542736991989e-15,-4.1542736991989e-15,6.9637942335018e-09,-1.5380914537945e-15,6.9638348915987e-09
+2300,167.42178774779,-4.7693767433386e-15,-4.7693767433386e-15,6.9634648319261e-09,1.6161914793874e-14,6.9635092836065e-09
+2400,169.88080414124,-1.5581085250685e-14,-1.5581085250685e-14,6.963746738717e-09,2.6741643186401e-14,6.9636963069794e-09
+2500,172.33982053468,2.0222652517556e-14,2.0222652517556e-14,6.9639201106559e-09,-2.2377918365695e-14,6.9639682596077e-09
+2600,174.79883692812,-3.0590134540558e-14,-3.0590134540558e-14,6.9637526116336e-09,2.6343933390704e-14,6.963701031394e-09
+2700,177.25785332156,-2.2419580671697e-15,-2.2419580671697e-15,6.9635944737489e-09,-1.895166254099e-14,6.9636435430007e-09
+2800,179.71686971501,7.6154899959538e-15,7.6154899959538e-15,6.9639816152359e-09,-1.537968402988e-14,6.963934886983e-09
+2900,182.17588610845,-1.8523869206836e-14,-1.8523869206836e-14,6.9640710763608e-09,-1.9243172541494e-14,6.9640698324491e-09
+3000,184.63490250189,-7.2713432677579e-15,-7.2713432677579e-15,6.963740263659e-09,-8.4562251552775e-15,6.9636892205663e-09
+3100,187.09391889533,-1.6847095456241e-14,-1.6847095456241e-14,6.9641093370878e-09,-1.6657545691793e-14,6.9641651179112e-09
+3200,189.55293528878,2.3823860021361e-14,2.3823860021361e-14,6.9632705538063e-09,-9.6446743835143e-15,6.9631964892197e-09
+3300,192.01195168222,2.3407597675145e-14,2.3407597675145e-14,6.963971379123e-09,3.2756565203238e-14,6.9639265199149e-09
+3400,194.47096807566,5.6478195843732e-14,5.6478195843732e-14,6.9637178737399e-09,1.4244329152432e-14,6.9636417127509e-09
+3500,196.9299844691,3.3783950933873e-14,3.3783950933873e-14,6.963371339804e-09,1.2608250688095e-13,6.9632283648268e-09
+3600,199.38900086255,-1.8587884877537e-14,-1.8587884877537e-14,6.9637332818012e-09,2.3902037531594e-14,6.9637100225812e-09
+3700,201.84801725599,2.2765463438567e-14,2.2765463438567e-14,6.9633330623874e-09,7.2920341902026e-14,6.9633277357738e-09
+3800,204.30703364943,-1.1275718360616e-14,-1.1275718360616e-14,6.9641435154832e-09,-5.1652112793475e-14,6.9639213172666e-09
+3900,206.76605004288,-4.506356964851e-15,-4.506356964851e-15,6.9638445865614e-09,2.9545874228886e-14,6.9639371769378e-09
+4000,209.22506643632,1.2084017214168e-14,1.2084017214168e-14,6.963612128577e-09,-4.1098265496684e-14,6.9636082395075e-09
+4100,211.68408282976,-1.4828657626715e-14,-1.4828657626715e-14,6.9635779711239e-09,4.0926162683013e-14,6.9635211576631e-09
+4200,214.1430992232,-3.0830448882258e-14,-3.0830448882258e-14,6.9638242165409e-09,3.6077647985645e-16,6.9639005037071e-09
+4300,216.60211561665,4.4468082982282e-15,4.4468082982282e-15,6.9634679448631e-09,-9.3801655019895e-15,6.9635109481199e-09
+4400,219.06113201009,1.0392344820131e-14,1.0392344820131e-14,6.9637534036181e-09,-5.6958476343321e-14,6.9636031681321e-09
+4500,221.52014840353,-5.5402657756632e-15,-5.5402657756632e-15,6.9638001383473e-09,1.7825664952173e-14,6.963878671229e-09
+4600,223.97916479697,-8.9813279046149e-14,-8.9813279046149e-14,6.9636781806244e-09,9.1513972059357e-14,6.963441905654e-09
+4700,226.43818119042,7.4707420930803e-14,7.4707420930803e-14,6.963577558665e-09,5.2428414563788e-14,6.9635003845163e-09
+4800,228.89719758386,-4.4485548962352e-14,-4.4485548962352e-14,6.9634227660137e-09,9.0872428047483e-14,6.963548523835e-09
+4900,231.3562139773,4.7751563746633e-15,4.7751563746633e-15,6.9633346627348e-09,2.5626492794367e-14,6.9635790723455e-09
+5000,233.81523037074,-4.7676049808003e-14,-4.7676049808003e-14,6.9637722886357e-09,-8.7820036976226e-14,6.9636783901121e-09
+5100,236.27424676419,2.2874358967347e-14,2.2874358967347e-14,6.9635503728131e-09,-2.2251329992502e-14,6.9636149736131e-09
+5200,238.73326315763,1.3394746495651e-14,1.3394746495651e-14,6.9636274845009e-09,1.0512521917358e-13,6.9636632892644e-09
+5300,241.19227955107,1.5179047623407e-14,1.5179047623407e-14,6.9635718280071e-09,-6.2327558230912e-14,6.9636502758307e-09
+5400,243.65129594451,-5.70409833069e-14,-5.70409833069e-14,6.9641516646837e-09,-5.632737051666e-15,6.9642545313127e-09
+5500,246.11031233796,2.5099749810567e-14,2.5099749810567e-14,6.9638142723791e-09,-8.2507257553391e-14,6.9638748851968e-09
+5600,248.5693287314,2.2334211508786e-14,2.2334211508786e-14,6.9641157859043e-09,-1.4709007591772e-13,6.964152226212e-09
+5700,251.02834512484,8.0169471943644e-14,8.0169471943644e-14,6.9633788765908e-09,-6.1904530115086e-14,6.9634629189392e-09
+5800,253.48736151829,1.2680150782889e-14,1.2680150782889e-14,6.963741602821e-09,9.3321758436442e-14,6.9638554996438e-09
+5900,255.94637791173,5.7641253539028e-15,5.7641253539028e-15,6.9639141530996e-09,8.7998458030306e-14,6.9640390159361e-09
+6000,258.40539430517,1.099667892919e-14,1.099667892919e-14,6.9631314891911e-09,-7.8290604449698e-14,6.9632100961353e-09
+6100,260.86441069861,6.6708256441928e-15,6.6708256441928e-15,6.963887108517e-09,-5.1263577848148e-14,6.9636331110095e-09
+6200,263.32342709206,-3.5085136974658e-15,-3.5085136974658e-15,6.9635449149246e-09,-1.5043696445693e-13,6.963815583394e-09
+6300,265.7824434855,6.8137125557923e-14,6.8137125557923e-14,6.9637212568786e-09,-8.0522231650014e-14,6.9637629411426e-09
+6400,268.24145987894,-6.6885369475721e-14,-6.6885369475721e-14,6.9635938808674e-09,2.4975782470794e-14,6.9631281846131e-09
+6500,270.70047627238,-4.1155282035822e-14,-4.1155282035822e-14,6.9635864048504e-09,-1.3330985810229e-14,6.9636650932038e-09
+6600,273.15949266583,2.0670427609053e-14,2.0670427609053e-14,6.9638457762971e-09,-2.3591095262453e-14,6.9640153228481e-09
+6700,275.61850905927,1.4732347258755e-14,1.4732347258755e-14,6.9634948435982e-09,1.74204013843e-13,6.9635986073261e-09
+6800,278.07752545271,-3.917632402377e-14,-3.917632402377e-14,6.9638716010759e-09,-5.7985805224599e-14,6.9638128563747e-09
+6900,280.53654184615,5.3164456864709e-14,5.3164456864709e-14,6.9645605982426e-09,-1.2345693116575e-13,6.9650182386501e-09
+7000,282.9955582396,-4.4592612846928e-14,-4.4592612846928e-14,6.9637919802837e-09,7.2629928558483e-14,6.9637047459068e-09
+7100,285.45457463304,4.2311498648118e-14,4.2311498648118e-14,6.964029431288e-09,-3.8383256550835e-14,6.9641039607568e-09
+7200,287.91359102648,-4.780586893781e-14,-4.780586893781e-14,6.9630838014753e-09,1.5607187506123e-13,6.9630080370769e-09
+7300,290.37260741992,-1.0684434941169e-13,-1.0684434941169e-13,6.9643931706747e-09,7.1058805595711e-14,6.9639697505934e-09
+7400,292.83162381337,3.2878212273166e-15,3.2878212273166e-15,6.9633973573211e-09,1.1205262641925e-13,6.9633821969787e-09
+7500,295.29064020681,1.1977945582102e-14,1.1977945582102e-14,6.9640507450511e-09,-1.9578992740771e-13,6.9641662338152e-09
+7600,297.74965660025,-1.3718388627065e-13,-1.3718388627065e-13,6.9637751523229e-09,-1.3867447250206e-13,6.9632134120176e-09
+7700,300.2086729937,4.7614805492503e-14,4.7614805492503e-14,6.9642206247798e-09,4.8651788032949e-14,6.964246095742e-09
+7800,302.66768938714,-1.1117363997486e-13,-1.1117363997486e-13,6.9637562914516e-09,1.5224403985293e-14,6.9644129678085e-09
+7900,305.12670578058,6.2632652195461e-14,6.2632652195461e-14,6.9637305173433e-09,7.7844487715741e-14,6.9636426444392e-09
+8000,307.58572217402,1.8301006671931e-14,1.8301006671931e-14,6.9637988228377e-09,6.0228627886447e-14,6.9637427695284e-09
+8100,310.04473856747,1.2457420085348e-13,1.2457420085348e-13,6.9637300838572e-09,-4.1806319086532e-13,6.9639516972278e-09
+8200,312.50375496091,-6.4595566459128e-14,-6.4595566459128e-14,6.9638193877879e-09,-1.4647717899192e-14,6.9639624267966e-09
+8300,314.96277135435,3.4939260173118e-14,3.4939260173118e-14,6.9640961899445e-09,3.0968961968318e-14,6.9643159629373e-09
+8400,317.42178774779,4.7138230443507e-14,4.7138230443507e-14,6.9638682680056e-09,2.6278469334863e-14,6.9633602959928e-09
+8500,319.88080414124,-1.1939930540245e-14,-1.1939930540245e-14,6.9639875061236e-09,1.9211596816672e-14,6.963780067291e-09
+8600,322.33982053468,-1.262136299607e-13,-1.262136299607e-13,6.9633536546583e-09,1.2589891961715e-14,6.9633415280813e-09
+8700,324.79883692812,5.01039924636e-14,5.01039924636e-14,6.9636984618808e-09,2.6205201802436e-13,6.9639510020091e-09
+8800,327.25785332156,-4.9075220229714e-14,-4.9075220229714e-14,6.9633767038732e-09,6.5913669594263e-14,6.9631454865076e-09
+8900,329.71686971501,1.1983186199014e-13,1.1983186199014e-13,6.9632924981844e-09,-7.7852613596752e-14,6.9638051435457e-09
+9000,332.17588610845,9.2052339637937e-15,9.2052339637937e-15,6.9646393675214e-09,-1.587940795611e-13,6.9644488629779e-09
+9100,334.63490250189,-6.2974480934262e-15,-6.2974480934262e-15,6.9637021337148e-09,-7.230184905894e-14,6.9646354624068e-09
+9200,337.09391889533,6.2144363193548e-14,6.2144363193548e-14,6.9641542309573e-09,2.5869319502838e-13,6.9638478817622e-09
+9300,339.55293528878,3.7287374061872e-14,3.7287374061872e-14,6.9626133392013e-09,1.6836387349155e-13,6.962944531119e-09
+9400,342.01195168222,7.9086701452755e-14,7.9086701452755e-14,6.9640325339146e-09,7.0885577905181e-14,6.9642539920217e-09
+9500,344.47096807566,-1.5164378589391e-13,-1.5164378589391e-13,6.9640047008042e-09,4.3320883439851e-14,6.9640210211773e-09
+9600,346.92998446911,-1.8122348012162e-14,-1.8122348012162e-14,6.9649133955625e-09,-4.6256348392741e-14,6.9645328249185e-09
+9700,349.38900086255,2.1368925065329e-14,2.1368925065329e-14,6.9633506289687e-09,1.5348685473794e-13,6.9632357841787e-09
+9800,351.84801725599,-2.4056791241032e-15,-2.4056791241032e-15,6.9632751230894e-09,-5.2057345980243e-14,6.9634625816101e-09
+9900,354.30703364943,-1.0701199978223e-13,-1.0701199978223e-13,6.9642347754984e-09,2.0086161190743e-13,6.9639536501859e-09
+10000,356.76605004287,6.2676404172943e-14,6.2676404172943e-14,6.9624001407291e-09,-2.229974578608e-13,6.962886817562e-09
diff --git a/test/tests/neml_complex/inelastic.i b/test/tests/neml_complex/inelastic.i
new file mode 100644
index 00000000..69fee12a
--- /dev/null
+++ b/test/tests/neml_complex/inelastic.i
@@ -0,0 +1,283 @@
+#
+#
+#
+[GlobalParams]
+      displacements = 'disp_x disp_y disp_z'
+[]
+
+[Variables]
+      [./disp_x]
+      [../]
+      [./disp_y]
+      [../]
+      [./disp_z]
+      [../]
+[]
+
+[Mesh]
+      type = GeneratedMesh
+      dim = 3
+[]
+
+[Kernels]
+      [./TensorMechanics]
+      [../]
+[]
+
+[AuxVariables]
+  [./strain_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+[]
+
+[AuxKernels]
+  [./strain_xx]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_zz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_yz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_xx]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./stress_yy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_zz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_yz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+[]
+
+[BCs]
+      [./fx]
+            type = PresetBC
+            variable = disp_x
+            boundary = 4
+            value = 0.0
+      [../]
+
+      [./fy]
+            type = PresetBC
+            variable = disp_y
+            boundary = 1
+            value = 0.0
+      [../]
+
+      [./fz]
+            type = PresetBC
+            variable = disp_z
+            boundary = 0
+            value = 0.0
+      [../]
+
+      [./px]
+            type = FunctionPresetBC
+            variable = disp_x
+            boundary = 2
+            function = pull
+      [../]
+[]
+
+[Functions]
+      [./pull]
+            type = PiecewiseLinear
+            x = '0 10000.0'
+            y = '0 0.1'
+      [../]
+[]
+
+[Materials]
+  [./strain]
+      type = ComputeSmallStrain
+  [../]
+
+  [./stress]
+    type = ComputeNEMLStress
+    database = 'complex.xml'
+    model = 'example'
+  [../]
+[]
+
+[Preconditioning]
+  [./pc]
+    type = SMP
+    full = True
+  [../]
+[]
+
+[Postprocessors]
+  [./stress_xx]
+      type = ElementAverageValue
+      variable = stress_xx
+  [../]
+  [./stress_yy]
+      type = ElementAverageValue
+      variable = stress_yy
+  [../]
+  [./stress_zz]
+      type = ElementAverageValue
+      variable = stress_zz
+  [../]
+  [./stress_xy]
+      type = ElementAverageValue
+      variable = stress_xy
+  [../]
+  [./stress_xz]
+      type = ElementAverageValue
+      variable = stress_xz
+  [../]
+  [./stress_yz]
+      type = ElementAverageValue
+      variable = stress_yz
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+
+  l_max_its = 5
+  l_tol = 1e-14
+  nl_max_its = 10
+  nl_rel_tol = 1e-10
+  nl_abs_tol = 1e-12
+
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  end_time = 10000.0
+  dt = 100.0
+[]
+
+[Outputs]
+  exodus = false
+  console = true
+  csv = true
+[]
diff --git a/test/tests/neml_complex/tests b/test/tests/neml_complex/tests
new file mode 100644
index 00000000..54f2f4c1
--- /dev/null
+++ b/test/tests/neml_complex/tests
@@ -0,0 +1,7 @@
+[Tests]
+  [./neml_inelastic]
+    type = 'CSVDiff'
+    input = 'inelastic.i'
+    csvdiff = 'inelastic_out.csv'
+  [../]
+[]
diff --git a/test/tests/neml_simple/examples.xml b/test/tests/neml_simple/examples.xml
new file mode 100644
index 00000000..2331d285
--- /dev/null
+++ b/test/tests/neml_simple/examples.xml
@@ -0,0 +1,14 @@
+<!-- Units are MPa, mm, s, and Kelvin -->
+<materials>
+  <!-- Nominally valid over 273 to 1073 K -->
+  <!-- Source: MCM fit to ASME and collected data -->
+  <le type="SmallStrainElasticity">
+    <elastic type="IsotropicLinearElasticModel">
+      <m1_type>youngs</m1_type>
+      <m1>100000</m1>
+      <m2_type>poissons</m2_type>
+      <m2>0.30</m2>
+    </elastic>
+    <alpha>1.0e-4</alpha>
+  </le>
+</materials>
diff --git a/test/tests/neml_simple/gold/le_stress_out.csv b/test/tests/neml_simple/gold/le_stress_out.csv
new file mode 100644
index 00000000..bcf3b2dc
--- /dev/null
+++ b/test/tests/neml_simple/gold/le_stress_out.csv
@@ -0,0 +1,3 @@
+time,stress_xx,stress_xy,stress_xz,stress_yy,stress_yz,stress_zz
+0,0,0,0,0,0,0
+1,1.477928890381e-12,-2.0850041778567e-15,-2.0850041778567e-15,1.2505552149378e-12,3.7530075201421e-14,10000
diff --git a/test/tests/neml_simple/gold/le_stress_thermal_out.csv b/test/tests/neml_simple/gold/le_stress_thermal_out.csv
new file mode 100644
index 00000000..5ff4224d
--- /dev/null
+++ b/test/tests/neml_simple/gold/le_stress_thermal_out.csv
@@ -0,0 +1,3 @@
+time,stress_xx,stress_xy,stress_xz,stress_yy,stress_yz,stress_zz
+0,0,0,0,0,0,0
+1,-3.4106051316485e-13,-7.6102652491771e-14,-7.6102652491771e-14,-7.105427357601e-14,-1.1988774022676e-13,-1000
diff --git a/test/tests/neml_simple/le_stress.i b/test/tests/neml_simple/le_stress.i
new file mode 100644
index 00000000..79197e98
--- /dev/null
+++ b/test/tests/neml_simple/le_stress.i
@@ -0,0 +1,283 @@
+#
+#
+#
+[GlobalParams]
+      displacements = 'disp_x disp_y disp_z'
+[]
+
+[Variables]
+      [./disp_x]
+      [../]
+      [./disp_y]
+      [../]
+      [./disp_z]
+      [../]
+[]
+
+[Mesh]
+      type = GeneratedMesh
+      dim = 3
+[]
+
+[Kernels]
+      [./TensorMechanics]
+      [../]
+[]
+
+[AuxVariables]
+  [./strain_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+[]
+
+[AuxKernels]
+  [./strain_xx]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_zz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_yz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_xx]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./stress_yy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_zz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_yz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+[]
+
+[BCs]
+      [./fx]
+            type = PresetBC
+            variable = disp_x
+            boundary = 4
+            value = 0.0
+      [../]
+
+      [./fy]
+            type = PresetBC
+            variable = disp_y
+            boundary = 1
+            value = 0.0
+      [../]
+
+      [./fz]
+            type = PresetBC
+            variable = disp_z
+            boundary = 0
+            value = 0.0
+      [../]
+
+      [./pz]
+            type = FunctionPresetBC
+            variable = disp_z
+            boundary = 5
+            function = pull
+      [../]
+[]
+
+[Functions]
+      [./pull]
+            type = PiecewiseLinear
+            x = '0 1'
+            y = '0 0.1'
+      [../]
+[]
+
+[Materials]
+  [./strain]
+      type = ComputeSmallStrain
+  [../]
+
+  [./stress]
+    type = ComputeNEMLStress
+    database = 'examples.xml'
+    model = 'le'
+  [../]
+[]
+
+[Preconditioning]
+  [./pc]
+    type = SMP
+    full = True
+  [../]
+[]
+
+[Postprocessors]
+  [./stress_xx]
+      type = ElementAverageValue
+      variable = stress_xx
+  [../]
+  [./stress_yy]
+      type = ElementAverageValue
+      variable = stress_yy
+  [../]
+  [./stress_zz]
+      type = ElementAverageValue
+      variable = stress_zz
+  [../]
+  [./stress_xy]
+      type = ElementAverageValue
+      variable = stress_xy
+  [../]
+  [./stress_xz]
+      type = ElementAverageValue
+      variable = stress_xz
+  [../]
+  [./stress_yz]
+      type = ElementAverageValue
+      variable = stress_yz
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+
+  l_max_its = 2
+  l_tol = 1e-3
+  nl_max_its = 10
+  nl_rel_tol = 1e-6
+  nl_abs_tol = 1e-8
+
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  end_time = 1.0
+  dt = 1.0
+[]
+
+[Outputs]
+  exodus = false
+  console = true
+  csv = true
+[]
diff --git a/test/tests/neml_simple/le_stress_thermal.i b/test/tests/neml_simple/le_stress_thermal.i
new file mode 100644
index 00000000..6ba893e7
--- /dev/null
+++ b/test/tests/neml_simple/le_stress_thermal.i
@@ -0,0 +1,302 @@
+#
+#
+#
+[GlobalParams]
+      displacements = 'disp_x disp_y disp_z'
+[]
+
+[Variables]
+      [./disp_x]
+      [../]
+      [./disp_y]
+      [../]
+      [./disp_z]
+      [../]
+[]
+
+[Mesh]
+      type = GeneratedMesh
+      dim = 3
+[]
+
+[Kernels]
+      [./TensorMechanics]
+      [../]
+[]
+
+[AuxVariables]
+  [./strain_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./strain_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./temp]
+      initial_condition = 0.0
+  [../]
+[]
+
+[AuxKernels]
+  [./strain_xx]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_zz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xy]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./strain_yz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./strain_xz]
+    type = RankTwoAux
+    rank_two_tensor = total_strain
+    variable = strain_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_xx]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xx
+    index_i = 0
+    index_j = 0
+    execute_on = timestep_end
+  [../]
+  [./stress_yy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_zz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_zz
+    index_i = 2
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xy
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./stress_yz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_yz
+    index_i = 1
+    index_j = 2
+    execute_on = timestep_end
+  [../]
+  [./stress_xz]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    variable = stress_xz
+    index_i = 0
+    index_j = 1
+    execute_on = timestep_end
+  [../]
+  [./temp]
+      type = FunctionAux
+      variable = temp
+      function = temperature
+  [../]
+[]
+
+[BCs]
+      [./fx]
+            type = PresetBC
+            variable = disp_x
+            boundary = 4
+            value = 0.0
+      [../]
+
+      [./fy]
+            type = PresetBC
+            variable = disp_y
+            boundary = 1
+            value = 0.0
+      [../]
+
+      [./fz]
+            type = PresetBC
+            variable = disp_z
+            boundary = 0
+            value = 0.0
+      [../]
+
+      [./pz]
+            type = PresetBC
+            variable = disp_z
+            boundary = 5
+            value = 0.0
+      [../]
+[]
+
+[Functions]
+      [./temperature]
+            type = PiecewiseLinear
+            x = '0 1'
+            y = '0 100.0'
+      [../]
+[]
+
+[Materials]
+  [./strain]
+      type = ComputeSmallStrain
+      eigenstrain_names = thermal
+  [../]
+
+  [./stress]
+    type = ComputeNEMLStress
+    database = 'examples.xml'
+    model = 'le'
+    temperature = temp
+  [../]
+
+  [./thermal_strain]
+      type = ComputeThermalExpansionEigenstrainNEML
+      database = 'examples.xml'
+      model = 'le'
+      temperature = temp
+      eigenstrain_name = thermal
+      stress_free_temperature = 0.0
+  [../]
+[]
+
+[Preconditioning]
+  [./pc]
+    type = SMP
+    full = True
+  [../]
+[]
+
+[Postprocessors]
+  [./stress_xx]
+      type = ElementAverageValue
+      variable = stress_xx
+  [../]
+  [./stress_yy]
+      type = ElementAverageValue
+      variable = stress_yy
+  [../]
+  [./stress_zz]
+      type = ElementAverageValue
+      variable = stress_zz
+  [../]
+  [./stress_xy]
+      type = ElementAverageValue
+      variable = stress_xy
+  [../]
+  [./stress_xz]
+      type = ElementAverageValue
+      variable = stress_xz
+  [../]
+  [./stress_yz]
+      type = ElementAverageValue
+      variable = stress_yz
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+
+  l_max_its = 2
+  l_tol = 1e-3
+  nl_max_its = 10
+  nl_rel_tol = 1e-6
+  nl_abs_tol = 1e-8
+
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  end_time = 1.0
+  dt = 1.0
+[]
+
+[Outputs]
+  exodus = false
+  console = true
+  csv = true
+[]
diff --git a/test/tests/neml_simple/tests b/test/tests/neml_simple/tests
new file mode 100644
index 00000000..398b519a
--- /dev/null
+++ b/test/tests/neml_simple/tests
@@ -0,0 +1,12 @@
+[Tests]
+  [./neml_linear_elastic]
+    type = 'CSVDiff'
+    input = 'le_stress.i'
+    csvdiff = 'le_stress_out.csv'
+  [../]
+  [./neml_linear_elastic_thermal]
+    type = 'CSVDiff'
+    input = 'le_stress_thermal.i'
+    csvdiff = 'le_stress_thermal_out.csv'
+  [../]
+[]