diff --git a/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_prefix.md b/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_prefix.md
new file mode 100644
index 0000000000000000000000000000000000000000..a0c291ee737ad8f5ef18e4adc8ac42964dee7c4d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_prefix.md
@@ -0,0 +1,5 @@
+Sets the options prefix for the PETSc SNES solver.
+
+See the PETSc
+[SNESSetOptionsPrefix](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESSetOptionsPrefix.html)
+documentation for details.
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 87892a7155e93eb679f1d664479da903102d02f5..910fcbd3cfc4d00a1e25d9a78449c7256ae41ee8 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -439,9 +439,13 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
 #ifdef USE_PETSC
     if (boost::iequals(type, "PETScSNES"))
     {
+        auto prefix =
+            //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__prefix}
+            config.getConfigParameter<std::string>("prefix", "");
         auto const tag = NonlinearSolverTag::Newton;
         using ConcreteNLS = PETScNonlinearSolver;
-        return std::make_pair(std::make_unique<ConcreteNLS>(linear_solver),
+        return std::make_pair(std::make_unique<ConcreteNLS>(
+                                  linear_solver, max_iter, std::move(prefix)),
                               tag);
     }
 
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index 64d63be03924a937b9ecbb0655bb3e0336b70d38..df470e69a4108108926304d80bd2b9f55c89310d 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -71,6 +71,10 @@ public:
     virtual void applyKnownSolutionsNewton(
         GlobalMatrix& Jac, GlobalVector& res,
         GlobalVector& minus_delta_x) const = 0;
+
+    virtual void updateConstraints(GlobalVector& /*lower*/,
+                                   GlobalVector& /*upper*/,
+                                   int /*process_id*/) = 0;
 };
 
 /*! A System of nonlinear equations to be solved with the Picard fixpoint
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 71e75aa25b64cb1813921838dbe93d86c054afef..9620dce5475eaa99a243fad4182b3c30d57ba368 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -66,6 +66,10 @@ public:
     {
         return nullptr;  // by default there are no known solutions
     }
+
+    virtual void updateConstraints(GlobalVector& /*lower*/,
+                                   GlobalVector& /*upper*/,
+                                   int const /*process_id*/){};
 };
 
 /*! Interface for a first-order implicit quasi-linear ODE.
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.cpp b/NumLib/ODESolver/PETScNonlinearSolver.cpp
index b536f8481806d29d77d8f654b75232ee6ba5edbf..dbbdd3b78a8032028b396f18f180d7eac62f35b4 100644
--- a/NumLib/ODESolver/PETScNonlinearSolver.cpp
+++ b/NumLib/ODESolver/PETScNonlinearSolver.cpp
@@ -48,6 +48,7 @@ PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
     context->system->getResidual(*context->x[context->process_id],
                                  *context->x_prev[context->process_id],
                                  *context->r);
+    context->r->finalizeAssembly();
     context->J->finalizeAssembly();
 
     context->system->getJacobian(*context->J);
@@ -77,9 +78,25 @@ PetscErrorCode updateJacobian(SNES /*snes*/, Vec /*x*/, Mat J,
 namespace NumLib
 {
 PETScNonlinearSolver::PETScNonlinearSolver(
-    GlobalLinearSolver& /*linear_solver*/)
+    GlobalLinearSolver& /*linear_solver*/, int const maxiter,
+    std::string prefix)
 {
     SNESCreate(PETSC_COMM_WORLD, &_snes_solver);
+    if (!prefix.empty())
+    {
+        prefix = prefix + "_";
+        SNESSetOptionsPrefix(_snes_solver, prefix.c_str());
+    }
+    // force SNESSolve() to take at least one iteration regardless of the
+    // initial residual norm
+    SNESSetForceIteration(_snes_solver, PETSC_TRUE);
+
+    // Set the maximum iterations.
+    PetscReal atol, rtol, stol;
+    PetscInt maxf;
+    SNESGetTolerances(_snes_solver, &atol, &rtol, &stol, nullptr, &maxf);
+    SNESSetTolerances(_snes_solver, atol, rtol, stol, maxiter, maxf);
+
     SNESSetFromOptions(_snes_solver);
 #ifndef NDEBUG
     PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_WORLD);
@@ -139,6 +156,15 @@ NonlinearSolverStatus PETScNonlinearSolver::solve(
     // during the SNES' solve call.
     auto& J_snes = NumLib::GlobalMatrixProvider::provider.getMatrix(
         system->getMatrixSpecifications(process_id), _petsc_jacobian_id);
+    BaseLib::RunTime timer_dirichlet;
+    double time_dirichlet = 0.0;
+
+    timer_dirichlet.start();
+    system->computeKnownSolutions(*x[process_id], process_id);
+    system->applyKnownSolutions(*x[process_id]);
+    time_dirichlet += timer_dirichlet.elapsed();
+    INFO("[time] Applying Dirichlet BCs took {} s.", time_dirichlet);
+
     auto& r_snes = NumLib::GlobalVectorProvider::provider.getVector(
         system->getMatrixSpecifications(process_id), _petsc_residual_id);
     auto& x_snes = NumLib::GlobalVectorProvider::provider.getVector(
@@ -156,6 +182,29 @@ NonlinearSolverStatus PETScNonlinearSolver::solve(
     SNESSetJacobian(_snes_solver, J_snes.getRawMatrix(), J_snes.getRawMatrix(),
                     updateJacobian, &petsc_context);
 
+    std::unique_ptr<GlobalVector> xl = nullptr;
+    std::unique_ptr<GlobalVector> xu = nullptr;
+
+    SNESType snes_type;
+    SNESGetType(_snes_solver, &snes_type);
+    if ((std::strcmp(snes_type, SNESVINEWTONRSLS) == 0) ||
+        (std::strcmp(snes_type, SNESVINEWTONSSLS) == 0))
+    {
+        // Set optional constraints via callback.
+        DBUG("PETScNonlinearSolver: set constraints");
+        xl = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            system->getMatrixSpecifications(process_id));
+        xu = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            system->getMatrixSpecifications(process_id));
+
+        system->updateConstraints(*xl, *xu, process_id);
+        MathLib::finalizeVectorAssembly(*xl);
+        MathLib::finalizeVectorAssembly(*xu);
+
+        SNESVISetVariableBounds(_snes_solver, xl->getRawVector(),
+                                xu->getRawVector());
+    }
+
     DBUG("PETScNonlinearSolver: call SNESSolve");
     SNESSolve(_snes_solver, nullptr, x_snes.getRawVector());
 
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.h b/NumLib/ODESolver/PETScNonlinearSolver.h
index 85b71434d4e2a7c1c7c24e93ab69e509a009343f..7f423e62facd20e5c91f108bc18993a2b7ef6260 100644
--- a/NumLib/ODESolver/PETScNonlinearSolver.h
+++ b/NumLib/ODESolver/PETScNonlinearSolver.h
@@ -30,8 +30,13 @@ public:
     /*! Constructs a new instance.
      *
      * \param linear_solver the linear solver used by this nonlinear solver.
+     * \param maxiter the maximum number of iterations used to solve the
+     *                equation.
+     * \param prefix used to set the SNES options.
      */
-    explicit PETScNonlinearSolver(GlobalLinearSolver& linear_solver);
+    PETScNonlinearSolver(GlobalLinearSolver& linear_solver,
+                         int maxiter,
+                         std::string prefix);
 
     //! Set the nonlinear equation system that will be solved.
     void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit);
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 066911818a0d444a583fc2145ddf8244a021c641..0264082d0e1be4f441a711b0512e5004a6c0938a 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -100,6 +100,12 @@ public:
     void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
                                    GlobalVector& minus_delta_x) const override;
 
+    void updateConstraints(GlobalVector& lower, GlobalVector& upper,
+                           int const process_id) override
+    {
+        _ode.updateConstraints(lower, upper, process_id);
+    }
+
     bool isLinear() const override { return _ode.isLinear(); }
 
     void preIteration(const unsigned iter, GlobalVector const& x) override
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index ce27a7c5e4ea2b7621135952673bb6dfa5b4a6a4..089190945330bf53b851c5b1371733049c7dc269 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -52,6 +52,9 @@ endif()
 if (OGS_USE_MPI)
     # OgsTest(WRAPPER mpirun -np 4 PROJECTFILE Mechanics/Linear/disc_with_hole.prj)
     OgsTest(WRAPPER mpirun -np 2 PROJECTFILE Mechanics/InitialStates/soil_column_nonequilibrium_sigma_elementwise.prj)
+    OgsTest(PROJECTFILE Mechanics/Ehlers/cube_1e0_SNES.prj)
+    OgsTest(PROJECTFILE Mechanics/Linear/square_1e0_SNES.prj)
+    OgsTest(PROJECTFILE Mechanics/Linear/cube_1e0_SNES.prj)
 endif()
 
 if (OGS_USE_MFRONT)
diff --git a/Tests/Data/Mechanics/Ehlers/cube_1e0_SNES.prj b/Tests/Data/Mechanics/Ehlers/cube_1e0_SNES.prj
new file mode 100644
index 0000000000000000000000000000000000000000..6cb996ecada7e225c32534ce5e974e562c3d673f
--- /dev/null
+++ b/Tests/Data/Mechanics/Ehlers/cube_1e0_SNES.prj
@@ -0,0 +1,340 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>cube_1x1x1_hex_1e0.vtu</mesh>
+    <geometry>cube_1x1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>Ehlers</type>
+                <shear_modulus>G</shear_modulus>
+                <bulk_modulus>K</bulk_modulus>
+                <kappa>kappa</kappa>
+                <beta>beta</beta>
+                <gamma>gamma</gamma>
+                <hardening_modulus>hard</hardening_modulus>
+                <alpha>alpha</alpha>
+                <delta>delta</delta>
+                <eps>epsilon</eps>
+                <m>m</m>
+                <alphap>alphap</alphap>
+                <deltap>deltap</deltap>
+                <epsp>epsilonp</epsp>
+                <mp>mp</mp>
+                <betap>betap</betap>
+                <gammap>gammap</gammap>
+                <tangent_type>Plastic</tangent_type>
+                <nonlinear_solver>
+                    <maximum_iterations>100</maximum_iterations>
+                    <residuum_tolerance>1e-14</residuum_tolerance>
+                    <increment_tolerance>0</increment_tolerance>
+                </nonlinear_solver>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0 0</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="sigma" output_name="sigma"/>
+                <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>petsc_snes</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-14</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>5.1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>0.05</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>1000</repeat>
+                            <delta_t>0.025</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>cube_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>10000000</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>G</name>
+            <type>Constant</type>
+            <value>150.</value>
+        </parameter>
+        <parameter>
+            <name>K</name>
+            <type>Constant</type>
+            <value>200.</value>
+        </parameter>
+        <parameter>
+            <name>kappa</name>
+            <type>Constant</type>
+            <value>0.1</value>
+        </parameter>
+        <parameter>
+            <name>beta</name>
+            <type>Constant</type>
+            <value>0.095</value>
+        </parameter>
+        <parameter>
+            <name>gamma</name>
+            <type>Constant</type>
+            <value>1.</value>
+        </parameter>
+        <parameter>
+            <name>hard</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>0.01</value>
+        </parameter>
+        <parameter>
+            <name>delta</name>
+            <type>Constant</type>
+            <value>0.0078</value>
+        </parameter>
+        <parameter>
+            <name>epsilon</name>
+            <type>Constant</type>
+            <value>0.1</value>
+        </parameter>
+        <parameter>
+            <name>m</name>
+            <type>Constant</type>
+            <value>0.54</value>
+        </parameter>
+        <parameter>
+            <name>alphap</name>
+            <type>Constant</type>
+            <value>0.01</value>
+        </parameter>
+        <parameter>
+            <name>deltap</name>
+            <type>Constant</type>
+            <value>0.0078</value>
+        </parameter>
+        <parameter>
+            <name>epsilonp</name>
+            <type>Constant</type>
+            <value>0.1</value>
+        </parameter>
+        <parameter>
+            <name>mp</name>
+            <type>Constant</type>
+            <value>0.54</value>
+        </parameter>
+        <parameter>
+            <name>betap</name>
+            <type>Constant</type>
+            <value>0.0608</value>
+        </parameter>
+        <parameter>
+            <name>gammap</name>
+            <type>Constant</type>
+            <value>1.</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0 0</values>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_left</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_bottom</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_front</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_top_spatial</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>Neumann_spatial</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_top</name>
+            <type>CurveScaled</type>
+            <curve>Dirichlet_top_temporal</curve>
+            <parameter>Dirichlet_top_spatial</parameter>
+        </parameter>
+        <parameter>
+            <name>Neumann_force_right</name>
+            <type>CurveScaled</type>
+            <curve>Neumann_temporal_right</curve>
+            <parameter>Neumann_spatial</parameter>
+        </parameter>
+        <parameter>
+            <name>Neumann_force_top</name>
+            <type>CurveScaled</type>
+            <curve>Neumann_temporal_top</curve>
+            <parameter>Neumann_spatial</parameter>
+        </parameter>
+    </parameters>
+    <curves>
+        <curve>
+            <name>Dirichlet_top_temporal</name>
+            <coords>0.0 0.1 5.1</coords>
+            <values>0.0 0.0 -0.005</values>
+        </curve>
+        <curve>
+            <name>Neumann_temporal_right</name>
+            <coords>0.0 0.1 5.1</coords>
+            <values>0.0 -0.07 -0.02</values>
+        </curve>
+        <curve>
+            <name>Neumann_temporal_top</name>
+            <coords>0.0 0.1 5.1</coords>
+            <values>0 -0.07 -0.02</values>
+        </curve>
+    </curves>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>3</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <!-- fixed boundaries -->
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>Dirichlet_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>front</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>Dirichlet_front</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>2</component>
+                    <parameter>Dirichlet_bottom</parameter>
+                </boundary_condition>
+                <!-- force -->
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>right</geometry>
+                    <type>Neumann</type>
+                    <component>0</component>
+                    <parameter>Neumann_force_right</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>back</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>Dirichlet_top</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>2</component>
+                    <parameter>Neumann_force_top</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>petsc_snes</name>
+            <type>PETScSNES</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <petsc>
+                <parameters>-snes_type newtonls -snes_rtol 1e-16 -ksp_type cg -pc_type bjacobi -ksp_rtol 1e-16</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <file>cube_1e0_ts_101_t_2.550000.vtu</file>
+            <field>displacement</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>0</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_101_t_2.550000.vtu</file>
+            <field>sigma</field>
+            <absolute_tolerance>1e-15</absolute_tolerance>
+            <relative_tolerance>0</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_203_t_5.100000.vtu</file>
+            <field>displacement</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>0</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_203_t_5.100000.vtu</file>
+            <field>sigma</field>
+            <absolute_tolerance>2e-15</absolute_tolerance>
+            <relative_tolerance>0</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_2_550000_0.vtu b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_2_550000_0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2b438af075f2cd2d873b48dd62aecaee8a03f455
--- /dev/null
+++ b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_2_550000_0.vtu
@@ -0,0 +1,27 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.0026676085653"      offset="0"                   />
+        <DataArray type="Float64" Name="sigma" NumberOfComponents="6" format="appended" RangeMin="0.28704467504"        RangeMax="0.28704467504"        offset="104"                 />
+        <DataArray type="Float64" Name="epsilon" NumberOfComponents="6" format="appended" RangeMin="0.0026676085653"      RangeMax="0.0026676085653"      offset="340"                 />
+      </PointData>
+      <CellData>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.7320508076"         offset="544"                 />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="608"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="672"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="712"                 />
+        <DataArray type="Int64" Name="faces" format="appended" RangeMin=""                     RangeMax=""                     offset="748"                 />
+        <DataArray type="Int64" Name="faceoffsets" format="appended" RangeMin=""                     RangeMax=""                     offset="788"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAACAAADAAAAAOgAAAA==eJxjYMAO4vtPsM8s9rDHIc2QtJhxD6tQyn4YPxGqHl0cHcRA1cWimR+DxoeZAxOPQzMfJg4ATG4mZA==AQAAAACAAACAAQAAngAAAA==eJzT+/9++THv5fvrNOfqOzy/uF8Dyhed+9DQ91HEHgY0YAKVb4Kqt4TynS33B766VG6Drt4IKt8IVQ/Tz35Sp+7g7AAM9QlQ+Raoeisof7tYUOWDy9kY7rGByjdD1TtA+ZdCzy7NqLbDMF8OTb0MlJ/p7L/FsToFw3yY+2vQwueetxHL4Wh7DPUCaOEpC/NPQ7b6XDcrDPcAAJHIkuQ=AQAAAACAAACAAQAAhQAAAA==eJxTbTjBPrPYwz5mMeMeVqGU/epQvp+aTont7ae7GdCAHlQ+GapeB8rfLaMonpfCYoOu3hQqnwZVbwTlW83ccOFrznFrdPUGUPlUqHoY3weHe0xwqHfCYb4FVD4Dqt4cynfBYb4SVD4Kql4eyk8wcD/RKr0SZ/gkQdXD+Co43AMAF1loEQ==AQAAAACAAADAAAAAHAAAAA==eJxjYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==AQAAAACAAABAAAAAHgAAAA==eJxjYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAdAQAAAACAAAAIAAAACwAAAA==eJzjYIAAAABIAAk=AQAAAACAAAABAAAACQAAAA==eJzjAQAADQANAQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=AQAAAACAAAAIAAAACwAAAA==eJxjZIAAAAAQAAI=
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_5_100000_0.vtu b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_5_100000_0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a0d73f72d2024becb7f290ac53b150c9727fce31
--- /dev/null
+++ b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_101_t_5_100000_0.vtu
@@ -0,0 +1,60 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8" NumberOfCells="1">
+      <PointData>
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="0.0058160307858">
+          AQAAAACAAADAAAAAMwAAAA==eJxjYMAOMqzUmCaZJtrjkGYoEFnn/rCqZD+6enRxXOaim4/Oh5mDrh4mngnlAwAVfyGm
+        </DataArray>
+        <DataArray type="Float64" Name="epsilon_xx" format="binary" RangeMin="0.002100738832" RangeMax="0.002100738832">
+          AQAAAACAAABAAAAAIQAAAA==eJxbVjtt0kTTRPslUHoZlF4FpVdA6dVQeg2aOgDTfR2Y
+        </DataArray>
+        <DataArray type="Float64" Name="epsilon_xy" format="binary" RangeMin="0" RangeMax="0">
+          AQAAAACAAABAAAAADAAAAA==eJxjYKAMAAAAQAAB
+        </DataArray>
+        <DataArray type="Float64" Name="epsilon_yy" format="binary" RangeMin="-0.005" RangeMax="-0.005">
+          AQAAAACAAABAAAAAHgAAAA==eJzLFVnn/rCqZH8hlC6A0llo4vlo/CIoDQAGLCA0
+        </DataArray>
+        <DataArray type="Float64" Name="epsilon_zz" format="binary" RangeMin="0.002100738832" RangeMax="0.002100738832">
+          AQAAAACAAABAAAAAHgAAAA==eJxbXTtt0kTTRPtlUHoFlF6Lxl8NpVeh0QDWBR2m
+        </DataArray>
+        <DataArray type="Float64" Name="sigma_xx" format="binary" RangeMin="-2.8228291448e-07" RangeMax="-2.8228291381e-07">
+          AQAAAACAAABAAAAAMwAAAA==eJwrk/iQNPXjpH0F9fVg+ua9S2Ba5u0SMK3vchFMn6k4BKalhOPAtMZ7mWQQDQDitSXs
+        </DataArray>
+        <DataArray type="Float64" Name="sigma_xy" format="binary" RangeMin="-9.7424738173e-18" RangeMax="6.1554594212e-17">
+          AQAAAACAAABAAAAASwAAAA==eJwBQAC//4qDstVptSc8rxihBPF2Zrx/Pfl8s4orPBL2dzf2hSC8vfvBv+O8Y7xs7ghp7b2RPFMvJVdECkG8B6FAKMpTRrwx3B+d
+        </DataArray>
+        <DataArray type="Float64" Name="sigma_yy" format="binary" RangeMin="-0.22776609819" RangeMax="-0.22776609819">
+          AQAAAACAAABAAAAAIwAAAA==eJy7aFnQXKB+dv9RKH0VSl+H0neg9EUovQ9K34bSACqsIQE=
+        </DataArray>
+        <DataArray type="Float64" Name="sigma_zz" format="binary" RangeMin="-2.8228291437e-07" RangeMax="-2.8228291366e-07">
+          AQAAAACAAABAAAAAMgAAAA==eJxzET6SNPXjpH1MtkFgWrB4GZg+vu4imA53+QumuY9OA9NdhkZget/692AaAHrAJSQ=
+        </DataArray>
+      </PointData>
+      <CellData>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="1.7320508076">
+          AQAAAACAAADAAAAAHAAAAA==eJxjYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="7">
+          AQAAAACAAABAAAAAHgAAAA==eJxjYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAd
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="binary" RangeMin="8" RangeMax="8">
+          AQAAAACAAAAIAAAACwAAAA==eJzjYIAAAABIAAk=
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="binary" RangeMin="12" RangeMax="12">
+          AQAAAACAAAABAAAACQAAAA==eJzjAQAADQAN
+        </DataArray>
+        <DataArray type="Int64" Name="faces" format="binary" RangeMin="0" RangeMax="0">
+          AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
+        </DataArray>
+        <DataArray type="Int64" Name="faceoffsets" format="binary" RangeMin="1" RangeMax="1">
+          AQAAAACAAAAIAAAACwAAAA==eJxjZIAAAAAQAAI=
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_203_t_5_100000_0.vtu b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_203_t_5_100000_0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9901c70ccef11123f5b911b0004cfffdddc6f9a7
--- /dev/null
+++ b/Tests/Data/Mechanics/Ehlers/cube_1e0_ts_203_t_5_100000_0.vtu
@@ -0,0 +1,27 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.0057046735735"      offset="0"                   />
+        <DataArray type="Float64" Name="sigma" NumberOfComponents="6" format="appended" RangeMin="0.25344987942"        RangeMax="0.25344987942"        offset="100"                 />
+        <DataArray type="Float64" Name="epsilon" NumberOfComponents="6" format="appended" RangeMin="0.0057046735735"      RangeMax="0.0057046735735"      offset="344"                 />
+      </PointData>
+      <CellData>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.7320508076"         offset="580"                 />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="644"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="708"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="748"                 />
+        <DataArray type="Int64" Name="faces" format="appended" RangeMin=""                     RangeMax=""                     offset="784"                 />
+        <DataArray type="Int64" Name="faceoffsets" format="appended" RangeMin=""                     RangeMax=""                     offset="824"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAACAAADAAAAAOQAAAA==eJxjYMAOfv3TObL0Yrw9DmmGapF17g+rSvbD+H+h6tHFcZn7F83872h8mDk/0NTDxH9C+QANeTgsAQAAAACAAACAAQAApQAAAA==eJxTD97g/rBqyn6V6/ETpspd2P8YyvcI/fja+1WeDQMaaAqByJtA1TNB+SHf0yusb9fvQVd/D2qeNlT9TCi/pzxnBmNZAoZ6F6h51lD1u6H88tDapKrz5RjuCYHKm0LV20D5wjcv/bmd34xh/pNgVPdvgvIFtepYdbsnYaivQvNvA5Q/N7a38aHqRAz3fIaapwdVD7PvtNRRzvaVKzDMBwAvUoqKAQAAAACAAACAAQAAnQAAAA==eJw78UznyNKL8fblIuvcH1aV7D8G5Z/4VyR8VuGpNQMauAKVr4eqvwDlz7h5ufuQCpMNuvqbUPk6qPrrUH5Nmr8+1+Qfu9HVX0JTD+N7q+mU2N4WxDD/BlS+Aaoe5r61PhEvqrYx70FXfwvN/TD36QmL7Zr3VgFD/RGofBFU/X4oXw2sXgzDPddg/kNz/4fpszPYVhpgmA8AlHyJZA==AQAAAACAAADAAAAAHAAAAA==eJxjYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==AQAAAACAAABAAAAAHgAAAA==eJxjYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAdAQAAAACAAAAIAAAACwAAAA==eJzjYIAAAABIAAk=AQAAAACAAAABAAAACQAAAA==eJzjAQAADQANAQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=AQAAAACAAAAIAAAACwAAAA==eJxjZIAAAAAQAAI=
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Mechanics/Linear/cube_1e0_SNES.prj b/Tests/Data/Mechanics/Linear/cube_1e0_SNES.prj
new file mode 100644
index 0000000000000000000000000000000000000000..dea66c1eb8a7e07c9c27c1685d77cb0b6b9b2339
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/cube_1e0_SNES.prj
@@ -0,0 +1,191 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>cube_1x1x1_hex_1e0.vtu</mesh>
+    <geometry>cube_1x1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0 0</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="sigma" output_name="sigma"/>
+                <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>petsc_snes</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-15</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>4</repeat>
+                            <delta_t>0.25</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>cube_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10000000</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.3</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0 0</values>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_left</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_bottom</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Dirichlet_front</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>Neumann_force_top</name>
+            <type>Constant</type>
+            <value>0.01</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>3</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <!-- fixed boundaries -->
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>Dirichlet_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>front</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>Dirichlet_front</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>2</component>
+                    <parameter>Dirichlet_bottom</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>2</component>
+                    <parameter>Neumann_force_top</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>petsc_snes</name>
+            <type>PETScSNES</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <petsc>
+                <parameters>-snes_type newtonls -ksp_type cg -pc_type bjacobi -ksp_rtol 1e-16 -ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <file>cube_1e0_ts_4_t_1_000000.vtu</file>
+            <field>displacement</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_4_t_1_000000.vtu</file>
+            <field>NodalForces</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_4_t_1_000000.vtu</file>
+            <field>sigma</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>cube_1e0_ts_4_t_1_000000.vtu</file>
+            <field>epsilon</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Linear/cube_1e0_ts_4_t_1_000000_0.vtu b/Tests/Data/Mechanics/Linear/cube_1e0_ts_4_t_1_000000_0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..6b7576d86b6e86a487e99f8011921784689d2982
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/cube_1e0_ts_4_t_1_000000_0.vtu
@@ -0,0 +1,26 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="Float64" Name="NodalForces" NumberOfComponents="3" format="appended" RangeMin="0.0025"               RangeMax="0.0025"               offset="0"                   />
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.010862780491"       offset="252"                 />
+        <DataArray type="Float64" Name="sigma" NumberOfComponents="6" format="appended" RangeMin="0.01"                 RangeMax="0.01"                 offset="356"                 />
+        <DataArray type="Float64" Name="epsilon" NumberOfComponents="6" format="appended" RangeMin="0.010862780491"       RangeMax="0.010862780491"       offset="712"                 />
+      </PointData>
+      <CellData>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.7320508076"         offset="1116"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1180"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="1244"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="1284"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAACAAADAAAAAqgAAAA==eJz7nm+ykXnSR2sl1xezZG4x21SLrHN/WJWy/8HDz74/TXZa74xJvrfe7PXuKqj4K+H7D5tW/91tq5r9bkX5C2uY+udLglOmHfq4W2r1hs/7Zf/Cxb/s1mFXOXh3t1XChI7chw9310DE7bX/TegQyXq+e6f4WhElxXPWUPPty15aST45wmiT/C2uZHU96x6Y+InTleqLTjy33jRfxL3H7fNuqPn2AGVnY24=AQAAAACAAADAAAAAOgAAAA==eJxjYMAOftZl7SmZnLEfhzTDbzT5H1A+IX01IuvcH1a12KOrq4KKo9tfCRX/C+X/gtK1UHEAsJI05w==AQAAAACAAACAAQAA+QAAAA==eJy79mtf+IRYa5uJQRksbza52lSJrHN/WNVin/X5zdaP6hd2M6ABh9XR8w0X++6RKnusunWdzZ5yqPoEhqsOTwI2WqOrd3PlenbG1HnP9C37D8oLB+ypgapf1jB5z4+dKjbo6huemRrumhJko3bk+o8tDf42TVD1WjcXTfSUk8ZQHxPFxbDZx9/Gkivy0kSWAJsWqPqS7Db7M2IcGOqNZGUnTP3otGdP9PmMVg2XPXVQ9fMh7tmDrn6JmHh98OeAPcVK7duCrgTuyYGqf7r5qF2OnRaG+hk627mmS9jbqPE4G91ca27TAFXvuIJx7QkFHgz3AADIRHTbAQAAAACAAACAAQAAGwEAAA==eJz7Vpe1p2Ryxv7vULpSZJ37w6oW+weLfC13MVzdnfFpmh3nYsU939gnH6ydcGH3L6g6GF0OVf9B0iPkUNM2axmTK1wqR//sXutlE35S5ps1Qz1EHSOUroGqX9rbE334o5ZNRuiSxydlD+5+5mGS+fyi7J5/UHP/QOkmqPovN1fOq1mgaBOxq7TuyENJGxexWQaZzkJ7fkPV/UNT76k6NehaAL+N39UN99+oWuxRuWibe/GDx56/UHVMUPfUQdXPgbhnz8oFZ/V/1rvucWRliVJ/52TzEar+DZTOg6p/v8W4/yef4Z69by6rWR9ytTkjFK4U89gUbj7MPQ1Q9T/cSgqaPIRsJKTaa46quts07Pdwv/vinTUApF/FlQ==AQAAAACAAADAAAAAHAAAAA==eJxjYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==AQAAAACAAABAAAAAHgAAAA==eJxjYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAdAQAAAACAAAAIAAAACwAAAA==eJzjYIAAAABIAAk=AQAAAACAAAABAAAACQAAAA==eJzjAQAADQAN
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Mechanics/Linear/square_1e0_SNES.prj b/Tests/Data/Mechanics/Linear/square_1e0_SNES.prj
new file mode 100644
index 0000000000000000000000000000000000000000..73a2ba31902dbda020da2287e0c9054c050d3621
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/square_1e0_SNES.prj
@@ -0,0 +1,159 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>square_1x1_quad_1e0.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="sigma" output_name="sigma"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>petsc</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-15</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>4</repeat>
+                            <delta_t>0.25</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>square_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10000000</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>sigma</variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.3</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>dirichlet1</name>
+            <type>Constant</type>
+            <value>0.05</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet1</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>petsc</name>
+            <type>PETScSNES</type>
+            <max_iter>4</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <petsc>
+                <parameters>-snes_type newtonls -ksp_type cg -pc_type bjacobi -ksp_rtol 1e-16 -ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <file>square_1e0_ts_4_t_1_000000.vtu</file>
+            <field>displacement</field>
+            <absolute_tolerance>1e-16</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>square_1e0_ts_4_t_1_000000.vtu</file>
+            <field>sigma</field>
+            <absolute_tolerance>1e-16</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Linear/square_1e0_ts_4_t_1_000000_0.vtu b/Tests/Data/Mechanics/Linear/square_1e0_ts_4_t_1_000000_0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..56e701545f49b2d0a0e175fbe31d158d4f18f6cb
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/square_1e0_ts_4_t_1_000000_0.vtu
@@ -0,0 +1,39 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="4" NumberOfCells="1">
+      <PointData>
+        <DataArray type="Float64" Name="NodalForces" NumberOfComponents="2" format="binary" RangeMin="0.027472527473" RangeMax="0.027472527473">
+          AQAAAACAAABAAAAALAAAAA==eJyTOvi/KXIpg43MIUUgmrNfCsLfA+Mvnvoo9SS/PEzeHsqHydsDAKCfG/U=
+        </DataArray>
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="2" format="binary" RangeMin="0" RangeMax="0.054398379328">
+          AQAAAACAAABAAAAAIAAAAA==eJxjYEAFYvEfReM/Tt2PJswwayYIrLSHycP4AFfTET0=
+        </DataArray>
+        <DataArray type="Float64" Name="sigma" NumberOfComponents="4" format="binary" RangeMin="0.057364321478" RangeMax="0.057364321478">
+          AQAAAACAAACAAAAAQQAAAA==eJzjXHn06P0rMjZihxRlDimusefleyjA93CCPQMUCEDlRaDyPGjy0UuVWIufWe0RwiEfBZWXgcoLoskDAIPGJjc=
+        </DataArray>
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="binary" RangeMin="0" RangeMax="0">
+          AQAAAACAAAAEAAAADAAAAA==eJxjYGBgAAAABAAB
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="1.4142135624">
+          AQAAAACAAABgAAAAFQAAAA==eJxjYMAHPtjjlcaQh/ER4gCW5AS9
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="3">
+          AQAAAACAAAAgAAAAEwAAAA==eJxjYIAARijNDKWZoDQAAHgABw==
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="binary" RangeMin="4" RangeMax="4">
+          AQAAAACAAAAIAAAACwAAAA==eJxjYYAAAAAoAAU=
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="binary" RangeMin="9" RangeMax="9">
+          AQAAAACAAAABAAAACQAAAA==eJzjBAAACgAK
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>