diff --git a/examples/structural/CMakeLists.txt b/examples/structural/CMakeLists.txt
index 13cda55..bca1e3a 100644
--- a/examples/structural/CMakeLists.txt
+++ b/examples/structural/CMakeLists.txt
@@ -1,4 +1,5 @@
 add_subdirectory(example_1) # membrane extension
 add_subdirectory(example_2) # mindlin plate bending
+add_subdirectory(example_4) # elastoplastic continuum analysis
 add_subdirectory(example_5) # SIMP min-compliance topology
 add_subdirectory(example_6) # SIMP min-compliance topology with MPI
diff --git a/examples/structural/example_1/example_1.cpp b/examples/structural/example_1/example_1.cpp
index 35a280e..49a38b8 100644
--- a/examples/structural/example_1/example_1.cpp
+++ b/examples/structural/example_1/example_1.cpp
@@ -19,17 +19,16 @@
 
 // MAST includes
 #include <mast/base/exceptions.hpp>
-#include <mast/util/perf_log.hpp>
+#include <mast/base/scalar_constant.hpp>
+#include <mast/base/assembly/libmesh/residual_and_jacobian.hpp>
+#include <mast/base/assembly/libmesh/residual_sensitivity.hpp>
 #include <mast/fe/eval/fe_basis_derivatives.hpp>
 #include <mast/fe/libmesh/fe_data.hpp>
 #include <mast/fe/libmesh/fe_side_data.hpp>
 #include <mast/fe/fe_var_data.hpp>
 #include <mast/physics/elasticity/isotropic_stiffness.hpp>
-#include <mast/base/scalar_constant.hpp>
 #include <mast/physics/elasticity/linear_strain_energy.hpp>
 #include <mast/physics/elasticity/pressure_load.hpp>
-#include <mast/base/assembly/libmesh/residual_and_jacobian.hpp>
-#include <mast/base/assembly/libmesh/residual_sensitivity.hpp>
 #include <mast/numerics/libmesh/sparse_matrix_initialization.hpp>
 
 // libMesh includes
@@ -37,6 +36,7 @@
 #include <libmesh/elem.h>
 #include <libmesh/mesh_generation.h>
 #include <libmesh/equation_systems.h>
+#include <libmesh/nonlinear_implicit_system.h>
 #include <libmesh/boundary_info.h>
 #include <libmesh/dirichlet_boundaries.h>
 #include <libmesh/zero_function.h>
@@ -94,9 +94,9 @@ class Context {
         delete mesh;
     }
     
-    uint_t elem_dim() const {return elem->dim();}
-    uint_t  n_nodes() const {return elem->n_nodes();}
-    real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
+    inline uint_t elem_dim() const {return elem->dim();}
+    inline uint_t  n_nodes() const {return elem->n_nodes();}
+    inline real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
     inline bool elem_is_quad() const {return (elem->type() == libMesh::QUAD4 ||
                                               elem->type() == libMesh::QUAD8 ||
                                               elem->type() == libMesh::QUAD9);}
@@ -133,9 +133,9 @@ struct Traits {
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
     using press_t           = typename MAST::Base::ScalarConstant<SolScalarType>;
     using area_t            = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t, Context>;
-    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, Dim, Context>;
-    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, Dim, Context>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, Dim>;
+    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, Dim>;
     using element_vector_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
     using element_matrix_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
     using assembled_vector_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
diff --git a/examples/structural/example_2/example_2.cpp b/examples/structural/example_2/example_2.cpp
index f2c4f6d..ac82428 100644
--- a/examples/structural/example_2/example_2.cpp
+++ b/examples/structural/example_2/example_2.cpp
@@ -19,17 +19,16 @@
 
 // MAST includes
 #include <mast/base/exceptions.hpp>
-#include <mast/util/perf_log.hpp>
+#include <mast/base/scalar_constant.hpp>
+#include <mast/base/assembly/libmesh/residual_and_jacobian.hpp>
+#include <mast/base/assembly/libmesh/residual_sensitivity.hpp>
 #include <mast/fe/eval/fe_basis_derivatives.hpp>
 #include <mast/fe/libmesh/fe_data.hpp>
 #include <mast/fe/fe_var_data.hpp>
 #include <mast/physics/elasticity/isotropic_stiffness.hpp>
-#include <mast/base/scalar_constant.hpp>
 #include <mast/physics/elasticity/mindlin_plate_strain_energy.hpp>
 #include <mast/physics/elasticity/plate_bending_section_property.hpp>
 #include <mast/physics/elasticity/shell_face_pressure_load.hpp>
-#include <mast/base/assembly/libmesh/residual_and_jacobian.hpp>
-#include <mast/base/assembly/libmesh/residual_sensitivity.hpp>
 #include <mast/numerics/libmesh/sparse_matrix_initialization.hpp>
 
 // libMesh includes
@@ -95,9 +94,9 @@ class Context {
         delete mesh;
     }
     
-    uint_t elem_dim() const {return elem->dim();}
-    uint_t  n_nodes() const {return elem->n_nodes();}
-    real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
+    inline uint_t elem_dim() const {return elem->dim();}
+    inline uint_t  n_nodes() const {return elem->n_nodes();}
+    inline real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
     inline bool elem_is_quad() const {return (elem->type() == libMesh::QUAD4 ||
                                               elem->type() == libMesh::QUAD8 ||
                                               elem->type() == libMesh::QUAD9);}
@@ -130,10 +129,10 @@ struct Traits {
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
     using press_t           = typename MAST::Base::ScalarConstant<SolScalarType>;
     using thickness_t       = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, 2, modulus_t, nu_t, Context>;
-    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<SolScalarType, material_t, thickness_t, Context>;
-    using energy_t          = typename MAST::Physics::Elasticity::MindlinPlate::StrainEnergy<fe_var_t, section_t, Context>;
-    using press_load_t      = typename MAST::Physics::Elasticity::ShellFacePressureLoad<fe_var_t, press_t, Context>;
+    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, 2, modulus_t, nu_t>;
+    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<SolScalarType, material_t, thickness_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::MindlinPlate::StrainEnergy<fe_var_t, section_t>;
+    using press_load_t      = typename MAST::Physics::Elasticity::ShellFacePressureLoad<fe_var_t, press_t>;
     using element_vector_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
     using element_matrix_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
     using assembled_vector_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
diff --git a/examples/structural/example_4/CMakeLists.txt b/examples/structural/example_4/CMakeLists.txt
new file mode 100644
index 0000000..74f2417
--- /dev/null
+++ b/examples/structural/example_4/CMakeLists.txt
@@ -0,0 +1,25 @@
+add_executable(structural_example_4
+                example_4.cpp)
+
+target_include_directories(structural_example_4 PRIVATE
+                            ${PROJECT_SOURCE_DIR}/include)
+
+target_link_libraries(structural_example_4 mast)
+
+install(TARGETS structural_example_4
+    RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/examples)
+
+# Test on single processor, PETSc built-in LU direct linear solver (sequential).
+add_test(NAME structural_example_4
+    COMMAND $<TARGET_FILE:structural_example_4> -ksp_type preonly -pc_type lu -options_view)
+set_tests_properties(structural_example_4
+                     PROPERTIES
+                     LABELS "SHORT;SEQ")
+
+# Test multiple processors, parallel direct linear solver using external MUMPS package.
+#add_test(NAME structural_example_4_mpi
+#    COMMAND ${MPIEXEC_EXECUTABLE} -np 2 $<TARGET_FILE:structural_example_4>
+#            -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -options_view)
+#set_tests_properties(structural_example_4_mpi
+#    PROPERTIES
+#        LABELS "SHORT;MPI")
diff --git a/examples/structural/example_4/example_4.cpp b/examples/structural/example_4/example_4.cpp
index 48c1b58..e6a8208 100644
--- a/examples/structural/example_4/example_4.cpp
+++ b/examples/structural/example_4/example_4.cpp
@@ -1,4 +1,570 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
 
-// BEGIN_TRANSLATE Structural Example 4
+// MAST includes
+#include <mast/base/exceptions.hpp>
+#include <mast/base/scalar_constant.hpp>
+#include <mast/base/material_point_system.hpp>
+#include <mast/base/material_point_data_storage.hpp>
+#include <mast/base/assembly/libmesh/residual_and_jacobian.hpp>
+#include <mast/base/assembly/libmesh/residual_sensitivity.hpp>
+#include <mast/fe/eval/fe_basis_derivatives.hpp>
+#include <mast/fe/libmesh/fe_data.hpp>
+#include <mast/fe/libmesh/fe_side_data.hpp>
+#include <mast/fe/fe_var_data.hpp>
+#include <mast/physics/elasticity/isotropic_stiffness.hpp>
+#include <mast/physics/elasticity/material_nonlinear_continuum_strain_energy.hpp>
+#include <mast/physics/elasticity/von_mises_yield_criterion.hpp>
+#include <mast/physics/elasticity/pressure_load.hpp>
+#include <mast/numerics/libmesh/sparse_matrix_initialization.hpp>
 
-// END_TRANSLATE
+// libMesh includes
+#include <libmesh/replicated_mesh.h>
+#include <libmesh/elem.h>
+#include <libmesh/mesh_generation.h>
+#include <libmesh/equation_systems.h>
+#include <libmesh/nonlinear_implicit_system.h>
+#include <libmesh/boundary_info.h>
+#include <libmesh/dirichlet_boundaries.h>
+#include <libmesh/zero_function.h>
+#include <libmesh/exodusII_io.h>
+
+// Eigen includes
+#include <Eigen/SparseLU>
+
+// BEGIN_TRANSLATE Elasto-plastic Continuum Analysis
+
+namespace MAST {
+namespace Examples {
+namespace Structural {
+namespace Example4 {
+
+class InitExample {
+    
+public:
+    
+    InitExample(libMesh::Parallel::Communicator& comm):
+    q_type    (libMesh::QGAUSS),
+    q_order   (libMesh::FOURTH),
+    fe_order  (libMesh::SECOND),
+    fe_family (libMesh::LAGRANGE),
+    mesh      (new libMesh::ReplicatedMesh(comm)),
+    eq_sys    (new libMesh::EquationSystems(*mesh)),
+    sys       (&eq_sys->add_system<libMesh::NonlinearImplicitSystem>("structural")),
+    p_side_id (1) {
+
+        libMesh::MeshTools::Generation::build_square(*mesh,
+                                                     20, 20,
+                                                     0.0, 10.0,
+                                                     0.0, 10.0,
+                                                     libMesh::QUAD9);
+
+        sys->add_variable("u_x", libMesh::FEType(fe_order, fe_family));
+        sys->add_variable("u_y", libMesh::FEType(fe_order, fe_family));
+
+        sys->get_dof_map().add_dirichlet_boundary
+        (libMesh::DirichletBoundary({3}, {0, 1}, libMesh::ZeroFunction<real_t>()));
+        
+        eq_sys->init();
+
+        mesh->print_info(std::cout);
+        eq_sys->print_info(std::cout);
+    }
+
+    virtual ~InitExample() {
+        
+        delete eq_sys;
+        delete mesh;
+    }
+    
+    libMesh::QuadratureType           q_type;
+    libMesh::Order                    q_order;
+    libMesh::Order                    fe_order;
+    libMesh::FEFamily                 fe_family;
+    libMesh::ReplicatedMesh          *mesh;
+    libMesh::EquationSystems         *eq_sys;
+    libMesh::NonlinearImplicitSystem *sys;
+    uint_t                            p_side_id;
+};
+
+
+
+template <typename TraitsType>
+class Context {
+    
+public:
+
+    Context(InitExample& init):
+    ex_init                      (init),
+    mesh                         (init.mesh),
+    eq_sys                       (init.eq_sys),
+    sys                          (init.sys),
+    elem                         (nullptr),
+    qp                           (-1),
+    fe                           (nullptr),
+    yield                        (nullptr),
+    mp_accessor                  (nullptr),
+    previous_plasticity_accessor (nullptr),
+    mp_sys                       (nullptr),
+    mp_storage                   (nullptr)
+    { }
+    
+    virtual ~Context() { }
+    
+    inline void init_for_qp(uint_t q) {
+        
+        qp = q;
+        
+        mp_accessor->init(*yield,
+                          mp_storage->value(elem->id(), qp).data());
+    }
+    
+    // assembly methods
+    inline uint_t  elem_dim() const {return elem->dim();}
+    inline uint_t  n_nodes() const {return elem->n_nodes();}
+    inline real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
+    inline real_t  qp_location(uint_t i) const { return fe->xyz(qp, i);}
+    inline bool elem_is_quad() const {return (elem->type() == libMesh::QUAD4 ||
+                                              elem->type() == libMesh::QUAD8 ||
+                                              elem->type() == libMesh::QUAD9);}
+    inline bool if_compute_pressure_load_on_side(const uint_t s)
+    { return ex_init.mesh->boundary_info->has_boundary_id(elem, s, ex_init.p_side_id);}
+    
+    
+    InitExample                         &ex_init;
+    libMesh::ReplicatedMesh             *mesh;
+    libMesh::EquationSystems            *eq_sys;
+    libMesh::NonlinearImplicitSystem    *sys;
+    libMesh::ExplicitSystem             *rho_sys;
+    const libMesh::Elem                 *elem;
+    uint_t                               qp;
+    typename TraitsType::fe_shape_t     *fe;
+    typename TraitsType::yield_t        *yield;
+    typename TraitsType::mp_accessor_t  *mp_accessor;
+    typename TraitsType::mp_accessor_t  *previous_plasticity_accessor;
+    typename TraitsType::mp_system_t    *mp_sys;
+    typename TraitsType::mp_storage_t   *mp_storage;
+};
+
+
+
+template <typename BasisScalarType,
+          typename NodalScalarType,
+          typename SolScalarType,
+          uint_t   Dim>
+struct Traits {
+
+    using scalar_t          = typename MAST::DeducedScalarType<typename MAST::DeducedScalarType<BasisScalarType, NodalScalarType>::type, SolScalarType>::type;
+    using traits_t          = MAST::Examples::Structural::Example4::Traits<BasisScalarType, NodalScalarType, SolScalarType, Dim>;
+    using context_t         = Context<traits_t>;
+    using fe_basis_t        = typename MAST::FEBasis::libMeshWrapper::FEBasis<BasisScalarType, Dim>;
+    using fe_shape_t        = typename MAST::FEBasis::Evaluation::FEShapeDerivative<BasisScalarType, NodalScalarType, Dim, Dim, fe_basis_t>;
+    using fe_data_t         = typename MAST::FEBasis::libMeshWrapper::FEData<Dim, fe_basis_t, fe_shape_t>;
+    using fe_side_data_t    = typename MAST::FEBasis::libMeshWrapper::FESideData<Dim, fe_basis_t, fe_shape_t>;
+    using fe_var_t          = typename MAST::FEBasis::FEVarData<BasisScalarType, NodalScalarType, SolScalarType, Dim, Dim, context_t, fe_shape_t>;
+    using modulus_t         = typename MAST::Base::ScalarConstant<SolScalarType>;
+    using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
+    using press_t           = typename MAST::Base::ScalarConstant<SolScalarType>;
+    using area_t            = typename MAST::Base::ScalarConstant<SolScalarType>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t>;
+    using yield_t           = typename MAST::Physics::Elasticity::ElastoPlasticity::vonMisesYieldFunction<scalar_t, prop_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::ElastoPlasticity::StrainEnergy<fe_var_t, yield_t, Dim>;
+    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, Dim>;
+    using element_vector_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
+    using element_matrix_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
+    using assembled_vector_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
+    using assembled_matrix_t = Eigen::SparseMatrix<scalar_t>;
+    using mp_accessor_t     = typename MAST::Physics::Elasticity::ElastoPlasticity::Accessor<scalar_t, yield_t>;
+    using mp_system_t        = MAST::Base::MaterialPointSystem;
+    using mp_storage_t       = MAST::Base::MaterialPointDataStorage<scalar_t>;
+};
+
+
+
+template <typename TraitsType>
+class ElemOps {
+  
+public:
+    
+    using scalar_t = typename TraitsType::scalar_t;
+    using vector_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
+    using matrix_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
+    
+
+    ElemOps(libMesh::Order          q_order,
+            libMesh::QuadratureType q_type,
+            libMesh::Order          fe_order,
+            libMesh::FEFamily       fe_family):
+    E             (nullptr),
+    nu            (nullptr),
+    press         (nullptr),
+    area          (nullptr),
+    _fe_data      (nullptr),
+    _fe_side_data (nullptr),
+    _fe_var       (nullptr),
+    _fe_side_var  (nullptr),
+    _prop         (nullptr),
+    _yield        (nullptr),
+    _energy       (nullptr),
+    _p_load       (nullptr) {
+        
+        _fe_data       = new typename TraitsType::fe_data_t;
+        _fe_data->init(q_order, q_type, fe_order, fe_family);
+        _fe_side_data  = new typename TraitsType::fe_side_data_t;
+        _fe_side_data->init(q_order, q_type, fe_order, fe_family);
+        _fe_var        = new typename TraitsType::fe_var_t;
+        _fe_side_var   = new typename TraitsType::fe_var_t;
+
+        // associate variables with the shape functions
+        _fe_var->set_fe_shape_data(_fe_data->fe_derivative());
+        _fe_side_var->set_fe_shape_data(_fe_side_data->fe_derivative());
+
+        // tell the FE computations which quantities are needed for computation
+        _fe_data->fe_basis().set_compute_dphi_dxi(true);
+        
+        _fe_data->fe_derivative().set_compute_dphi_dx(true);
+        _fe_data->fe_derivative().set_compute_detJxW(true);
+        
+        _fe_side_data->fe_basis().set_compute_dphi_dxi(true);
+        _fe_side_data->fe_derivative().set_compute_normal(true);
+        _fe_side_data->fe_derivative().set_compute_detJxW(true);
+
+        _fe_var->set_compute_du_dx(true);
+        
+        // variables for physics
+        E        = new typename TraitsType::modulus_t(72.e9);
+        nu       = new typename TraitsType::nu_t(0.33);
+        press    = new typename TraitsType::press_t(1.e2);
+        area     = new typename TraitsType::area_t(1.0);
+        _prop    = new typename TraitsType::prop_t;
+        _yield   = new typename TraitsType::yield_t;
+        
+        _prop->set_modulus_and_nu(*E, *nu);
+        _yield->set_limit_stress(1.e6);
+        _yield->set_material(*_prop);
+        _energy   = new typename TraitsType::energy_t;
+        _energy->set_yield_surface(*_yield);
+        _p_load   = new typename TraitsType::press_load_t;
+        _p_load->set_section_area(*area);
+        _p_load->set_pressure(*press);
+        
+        // tell physics kernels about the FE discretization information
+        _energy->set_fe_var_data(*_fe_var);
+        _p_load->set_fe_var_data(*_fe_side_var);
+    }
+    
+    virtual ~ElemOps() {
+        
+        delete _p_load;
+        delete area;
+        delete press;
+        delete _energy;
+        delete _prop;
+        delete nu;
+        delete E;
+        delete _fe_var;
+        delete _fe_side_var;
+        delete _fe_side_data;
+        delete _fe_data;
+    }
+    
+
+    template <typename ContextType, typename AccessorType>
+    inline void compute(ContextType                       &c,
+                        const AccessorType                &v,
+                        typename TraitsType::element_vector_t &res,
+                        typename TraitsType::element_matrix_t *jac) {
+        
+        _fe_data->reinit(c);
+        _fe_var->init(c, v);
+        _energy->compute(c, res, jac);
+        
+        for (uint_t s=0; s<c.elem->n_sides(); s++)
+            if (c.if_compute_pressure_load_on_side(s)) {
+                                
+                _fe_side_data->reinit_for_side(c, s);
+                _fe_side_var->init(c, v);
+                _p_load->compute(c, res, jac);
+            }
+    }
+
+    
+    template <typename ContextType, typename AccessorType, typename ScalarFieldType>
+    inline void derivative(ContextType                       &c,
+                           const ScalarFieldType             &f,
+                           const AccessorType                &v,
+                           typename TraitsType::element_vector_t &res,
+                           typename TraitsType::element_matrix_t *jac) {
+        
+        _fe_data->reinit(c);
+        _fe_var->init(c, v);
+        _energy->derivative(c, f, res, jac);
+        
+        for (uint_t s=0; s<c.elem->n_sides(); s++)
+            if (c.if_compute_pressure_load_on_side(s)) {
+                                
+                _fe_side_data->reinit_for_side(c, s);
+                _fe_side_var->init(c, v);
+                _p_load->derivative(c, f, res, jac);
+            }
+    }
+
+    // parameters
+    typename TraitsType::modulus_t    *E;
+    typename TraitsType::nu_t         *nu;
+    typename TraitsType::press_t      *press;
+    typename TraitsType::area_t       *area;
+    
+private:
+
+    // variables for quadrature and shape function
+    typename TraitsType::fe_data_t         *_fe_data;
+    typename TraitsType::fe_side_data_t    *_fe_side_data;
+    typename TraitsType::fe_var_t          *_fe_var;
+    typename TraitsType::fe_var_t          *_fe_side_var;
+    typename TraitsType::prop_t            *_prop;
+    typename TraitsType::yield_t           *_yield;
+    typename TraitsType::energy_t          *_energy;
+    typename TraitsType::press_load_t      *_p_load;
+};
+
+
+template <typename TraitsType>
+inline void
+compute_residual(Context<TraitsType>                            &c,
+                 ElemOps<TraitsType>                            &e_ops,
+                 const typename TraitsType::assembled_vector_t  &sol,
+                 typename TraitsType::assembled_vector_t        &res) {
+    
+    using scalar_t   = typename TraitsType::scalar_t;
+
+    MAST::Base::Assembly::libMeshWrapper::ResidualAndJacobian<scalar_t, ElemOps<TraitsType>>
+    assembly;
+    
+    assembly.set_elem_ops(e_ops);
+
+    typename TraitsType::assembled_matrix_t
+    *jac = nullptr;
+    
+    res = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    
+    assembly.assemble(c, sol, &res, jac);
+}
+
+
+
+template <typename TraitsType, typename ScalarFieldType>
+inline void
+compute_residual_sensitivity(Context<TraitsType>                            &c,
+                             ElemOps<TraitsType>                            &e_ops,
+                             const ScalarFieldType                          &f,
+                             const typename TraitsType::assembled_vector_t  &sol,
+                             typename TraitsType::assembled_vector_t        &dres) {
+    
+    using scalar_t   = typename TraitsType::scalar_t;
+
+    MAST::Base::Assembly::libMeshWrapper::ResidualAndJacobian<scalar_t, ElemOps<TraitsType>>
+    assembly;
+    
+    assembly.set_elem_ops(e_ops);
+
+    typename TraitsType::assembled_matrix_t
+    *jac = nullptr;
+    
+    dres = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    
+    assembly.assemble(c, f, sol, &dres, jac);
+}
+
+
+
+template <typename TraitsType>
+inline void
+compute_sol(Context<TraitsType>                      &c,
+            ElemOps<TraitsType>                      &e_ops,
+            typename TraitsType::assembled_vector_t  &sol) {
+    
+    using scalar_t   = typename TraitsType::scalar_t;
+
+    MAST::Base::Assembly::libMeshWrapper::ResidualAndJacobian<scalar_t, ElemOps<TraitsType>>
+    assembly;
+    
+    assembly.set_elem_ops(e_ops);
+
+    typename TraitsType::assembled_vector_t
+    res;
+    typename TraitsType::assembled_matrix_t
+    jac;
+    
+    sol = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    res = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    MAST::Numerics::libMeshWrapper::init_sparse_matrix(c.sys->get_dof_map(), jac);
+    
+    assembly.assemble(c, sol, &res, &jac);
+    
+    sol = Eigen::SparseLU<typename TraitsType::assembled_matrix_t>(jac).solve(-res);
+}
+
+
+template <typename TraitsType, typename ScalarFieldType>
+inline void
+compute_sol_sensitivity(Context<TraitsType>                            &c,
+                        ElemOps<TraitsType>                            &e_ops,
+                        const ScalarFieldType                          &f,
+                        const typename TraitsType::assembled_vector_t  &sol,
+                        typename TraitsType::assembled_vector_t        &dsol) {
+    
+    using scalar_t   = typename TraitsType::scalar_t;
+    
+    typename TraitsType::assembled_vector_t
+    *res = nullptr,
+    dres;
+    typename TraitsType::assembled_matrix_t
+    jac,
+    *djac = nullptr;
+    
+    dsol = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    dres = TraitsType::assembled_vector_t::Zero(c.sys->n_dofs());
+    MAST::Numerics::libMeshWrapper::init_sparse_matrix(c.sys->get_dof_map(), jac);
+
+    // assembly of Jacobian matrix
+    {
+        MAST::Base::Assembly::libMeshWrapper::ResidualAndJacobian<scalar_t, ElemOps<TraitsType>>
+        assembly;
+        assembly.set_elem_ops(e_ops);
+        assembly.assemble(c, sol, res, &jac);
+    }
+
+    // assembly of sensitivity RHS
+    {
+        MAST::Base::Assembly::libMeshWrapper::ResidualSensitivity<scalar_t, ElemOps<TraitsType>>
+        sens_assembly;
+        sens_assembly.set_elem_ops(e_ops);
+        sens_assembly.assemble(c, f, sol, &dres, djac);
+    }
+    
+    dsol = Eigen::SparseLU<typename TraitsType::assembled_matrix_t>(jac).solve(-dres);
+}
+
+
+} // namespace Example4
+} // namespace Structural
+} // namespace Examples
+} // namespace MAST
+
+#ifndef MAST_TESTING
+
+int main(int argc, const char** argv) {
+
+    libMesh::LibMeshInit init(argc, argv);
+    
+    using traits_t           = MAST::Examples::Structural::Example4::Traits<real_t, real_t,    real_t, 2>;
+    using traits_complex_t   = MAST::Examples::Structural::Example4::Traits<real_t, real_t, complex_t, 2>;
+
+
+    MAST::Examples::Structural::Example4::InitExample
+    ex_init(init.comm());
+    
+    MAST::Examples::Structural::Example4::Context<traits_t> c(ex_init);
+    MAST::Examples::Structural::Example4::Context<traits_complex_t> c_c(ex_init);
+    
+    MAST::Examples::Structural::Example4::ElemOps<traits_t>
+    e_ops(ex_init.q_order, ex_init.q_type, ex_init.fe_order, ex_init.fe_family);
+    MAST::Examples::Structural::Example4::ElemOps<traits_complex_t>
+    e_ops_c(ex_init.q_order, ex_init.q_type, ex_init.fe_order, ex_init.fe_family);
+
+    typename traits_t::assembled_vector_t
+    sol,
+    dsol;
+
+    typename traits_complex_t::assembled_vector_t
+    sol_c;
+
+    // compute the solution
+    MAST::Examples::Structural::Example4::compute_sol<traits_t>(c, e_ops, sol);
+    
+    // write solution as first time-step
+    libMesh::ExodusII_IO writer(*c.mesh);
+    {
+        for (uint_t i=0; i<sol.size(); i++) c.sys->solution->set(i, sol(i));
+        writer.write_timestep("solution.exo", *c.eq_sys, 1, 1.);
+    }
+
+    // compute the solution sensitivity wrt E
+    (*e_ops_c.E)() += complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol<traits_complex_t>(c_c, e_ops_c, sol_c);
+    (*e_ops_c.E)() -= complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol_sensitivity<traits_t>(c, e_ops, *e_ops.E, sol, dsol);
+    
+    // write solution as first time-step
+    {
+        for (uint_t i=0; i<sol.size(); i++) c.sys->solution->set(i, dsol(i));
+        writer.write_timestep("solution.exo", *c.eq_sys, 2, 2.);
+    }
+    dsol -= sol_c.imag()/ComplexStepDelta;
+    std::cout << dsol.norm() << std::endl;
+
+    // compute the solution sensitivity wrt nu
+    (*e_ops_c.nu)() += complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol<traits_complex_t>(c_c, e_ops_c, sol_c);
+    (*e_ops_c.nu)() -= complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol_sensitivity<traits_t>(c, e_ops, *e_ops.nu, sol, dsol);
+    
+    // write solution as first time-step
+    {
+        for (uint_t i=0; i<sol.size(); i++) c.sys->solution->set(i, dsol(i));
+        writer.write_timestep("solution.exo", *c.eq_sys, 3, 3.);
+    }
+    dsol -= sol_c.imag()/ComplexStepDelta;
+    std::cout << dsol.norm() << std::endl;
+
+    // compute the solution sensitivity wrt p
+    (*e_ops_c.press)() += complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol<traits_complex_t>(c_c, e_ops_c, sol_c);
+    (*e_ops_c.press)() -= complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol_sensitivity<traits_t>(c, e_ops, *e_ops.press, sol, dsol);
+    
+    // write solution as first time-step
+    {
+        for (uint_t i=0; i<sol.size(); i++) c.sys->solution->set(i, dsol(i));
+        writer.write_timestep("solution.exo", *c.eq_sys, 4, 4.);
+    }
+    dsol -= sol_c.imag()/ComplexStepDelta;
+    std::cout << dsol.norm() << std::endl;
+
+    // compute the solution sensitivity wrt section area
+    (*e_ops_c.area)() += complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol<traits_complex_t>(c_c, e_ops_c, sol_c);
+    (*e_ops_c.area)() -= complex_t(0., ComplexStepDelta);
+    MAST::Examples::Structural::Example4::compute_sol_sensitivity<traits_t>(c, e_ops, *e_ops.area, sol, dsol);
+    
+    // write solution as first time-step
+    {
+        for (uint_t i=0; i<sol.size(); i++) c.sys->solution->set(i, dsol(i));
+        writer.write_timestep("solution.exo", *c.eq_sys, 5, 5.);
+    }
+    dsol -= sol_c.imag()/ComplexStepDelta;
+    std::cout << dsol.norm() << std::endl;
+
+    
+    // END_TRANSLATE
+    return 0;
+}
+
+#endif // MAST_TESTING
diff --git a/examples/structural/example_5/example_5.cpp b/examples/structural/example_5/example_5.cpp
index 6c2645e..33d20ac 100644
--- a/examples/structural/example_5/example_5.cpp
+++ b/examples/structural/example_5/example_5.cpp
@@ -203,9 +203,9 @@ struct Traits {
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
     using press_t           = typename ModelType::template pressure_t<scalar_t>;
     using area_t            = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, dim, modulus_t, nu_t, context_t>;
-    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, dim, context_t>;
-    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, dim, context_t>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, dim, modulus_t, nu_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, dim>;
+    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, dim>;
     using element_vector_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
     using element_matrix_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
     using assembled_vector_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp
index 4c33ea0..6ea5018 100644
--- a/examples/structural/example_6/example_6.cpp
+++ b/examples/structural/example_6/example_6.cpp
@@ -19,6 +19,7 @@
 
 // MAST includes
 #include <mast/base/exceptions.hpp>
+#include <mast/base/scalar_constant.hpp>
 #include <mast/util/perf_log.hpp>
 #include <mast/fe/eval/fe_basis_derivatives.hpp>
 #include <mast/fe/libmesh/fe_data.hpp>
@@ -26,7 +27,6 @@
 #include <mast/fe/fe_var_data.hpp>
 #include <mast/fe/scalar_field_wrapper.hpp>
 #include <mast/physics/elasticity/isotropic_stiffness.hpp>
-#include <mast/base/scalar_constant.hpp>
 #include <mast/physics/elasticity/linear_strain_energy.hpp>
 #include <mast/physics/elasticity/pressure_load.hpp>
 #include <mast/optimization/topology/simp/penalized_density.hpp>
@@ -156,10 +156,10 @@ class Context {
     virtual ~Context() { }
     
     // assembly methods
-    uint_t  elem_dim() const {return elem->dim();}
-    uint_t  n_nodes() const {return elem->n_nodes();}
-    real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
-    real_t  qp_location(uint_t i) const { return fe->xyz(qp, i);}
+    inline uint_t  elem_dim() const {return elem->dim();}
+    inline uint_t  n_nodes() const {return elem->n_nodes();}
+    inline real_t  nodal_coord(uint_t nd, uint_t c) const {return elem->point(nd)(c);}
+    inline real_t  qp_location(uint_t i) const { return fe->xyz(qp, i);}
     inline bool elem_is_quad() const {return (elem->type() == libMesh::QUAD4 ||
                                               elem->type() == libMesh::QUAD8 ||
                                               elem->type() == libMesh::QUAD9);}
@@ -203,9 +203,9 @@ struct Traits {
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
     using press_t           = typename ModelType::template pressure_t<scalar_t>;
     using area_t            = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, dim, modulus_t, nu_t, context_t>;
-    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, dim, context_t>;
-    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, dim, context_t>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, dim, modulus_t, nu_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, dim>;
+    using press_load_t      = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, press_t, area_t, dim>;
     using element_vector_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
     using element_matrix_t  = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
     using assembled_vector_t = libMesh::NumericVector<scalar_t>;
diff --git a/include/mast/base/assembly/libmesh/material_point_residual_and_jacobian.hpp b/include/mast/base/assembly/libmesh/material_point_residual_and_jacobian.hpp
new file mode 100644
index 0000000..943c970
--- /dev/null
+++ b/include/mast/base/assembly/libmesh/material_point_residual_and_jacobian.hpp
@@ -0,0 +1,137 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_libmesh_material_point_residual_and_jacobian_h__
+#define __mast_libmesh_material_point_residual_and_jacobian_h__
+
+// MAST includes
+#include <mast/base/mast_data_types.h>
+#include <mast/base/exceptions.hpp>
+#include <mast/base/assembly/libmesh/utility.hpp>
+#include <mast/base/assembly/libmesh/accessor.hpp>
+#include <mast/numerics/utility.hpp>
+
+// libMesh includes
+#include <libmesh/nonlinear_implicit_system.h>
+#include <libmesh/dof_map.h>
+
+namespace MAST {
+namespace Base {
+namespace Assembly {
+namespace libMeshWrapper {
+
+template <typename ScalarType,
+          typename ElemOpsType>
+class MaterialPointResidualAndJacobian {
+
+public:
+    
+    static_assert(std::is_same<ScalarType, typename ElemOpsType::scalar_t>::value,
+                  "Scalar type of assembly and element operations must be same");
+    
+    MaterialPointResidualAndJacobian():
+    _finalize_jac (true),
+    _e_ops        (nullptr)
+    { }
+    
+    virtual ~MaterialPointResidualAndJacobian() { }
+    
+    inline void set_finalize_jac(bool f) { _finalize_jac = f;}
+    
+    inline void set_elem_ops(ElemOpsType& e_ops) { _e_ops = &e_ops; }
+
+    template <typename ContextType, typename VecType, typename MatType>
+    inline void assemble(ContextType   &c,
+                         const VecType &X,
+                         VecType       *R,
+                         MatType       *J) {
+        
+        Assert0( R || J, "Atleast one assembled quantity should be specified.");
+        
+        if (R) MAST::Numerics::Utility::setZero(*R);
+        if (J) MAST::Numerics::Utility::setZero(*J);
+        
+        // iterate over each element, initialize it and get the relevant
+        // analysis quantities
+        typename MAST::Base::Assembly::libMeshWrapper::Accessor<ScalarType, VecType>
+        sol_accessor(*c.sys, X);
+        typename MAST::Base::Assembly::MaterialPointAccessor<ScalarType>
+        mp_accessor(*c.mp_sys, X);
+
+        using elem_vector_t = typename ElemOpsType::vector_t;
+        using elem_matrix_t = typename ElemOpsType::matrix_t;
+        
+        elem_vector_t res_e;
+        elem_matrix_t jac_e;
+        
+        
+        libMesh::MeshBase::const_element_iterator
+        el     = c.mesh->active_local_elements_begin(),
+        end_el = c.mesh->active_local_elements_end();
+        
+        for ( ; el != end_el; ++el) {
+            
+            // set element in the context, which will be used for the initialization routines
+            c.elem = *el;
+            
+            sol_accessor.init(*c.elem);
+            
+            res_e.setZero(sol_accessor.n_dofs());
+            if (J) jac_e.setZero(sol_accessor.n_dofs(), sol_accessor.n_dofs());
+            
+            // perform the element level calculations
+            _e_ops->compute(c,
+                            sol_accessor,
+                            mp_accessor,
+                            res_e,
+                            J?&jac_e:nullptr);
+                        
+            // constrain the quantities to account for hanging dofs,
+            // Dirichlet constraints, etc.
+            if (R && J)
+                MAST::Base::Assembly::libMeshWrapper::constrain_and_add_matrix_and_vector
+                <ScalarType, VecType, MatType, elem_vector_t, elem_matrix_t>
+                (*R, *J, c.sys->get_dof_map(), sol_accessor.dof_indices(), res_e, jac_e);
+            else if (R)
+                MAST::Base::Assembly::libMeshWrapper::constrain_and_add_vector
+                <ScalarType, VecType, elem_vector_t>
+                (*R, c.sys->get_dof_map(), sol_accessor.dof_indices(), res_e);
+            else
+                MAST::Base::Assembly::libMeshWrapper::constrain_and_add_matrix
+                <ScalarType, MatType, elem_matrix_t>
+                (*J, c.sys->get_dof_map(), sol_accessor.dof_indices(), jac_e);
+        }
+
+        // parallel matrix/vector require finalization of communication
+        if (R) MAST::Numerics::Utility::finalize(*R);
+        if (J && _finalize_jac) MAST::Numerics::Utility::finalize(*J);
+    }
+    
+private:
+
+    bool         _finalize_jac;
+    ElemOpsType  *_e_ops;
+};
+
+} // namespace libMeshWrapper
+} // namespace Assembly
+} // namespace Base
+} // namespace MAST
+
+#endif // __mast_libmesh_material_point_residual_and_jacobian_h__
diff --git a/include/mast/base/assembly/libmesh/residual_and_jacobian.hpp b/include/mast/base/assembly/libmesh/residual_and_jacobian.hpp
index 3f39e49..355e1a4 100644
--- a/include/mast/base/assembly/libmesh/residual_and_jacobian.hpp
+++ b/include/mast/base/assembly/libmesh/residual_and_jacobian.hpp
@@ -56,7 +56,7 @@ class ResidualAndJacobian {
     
     inline void set_elem_ops(ElemOpsType& e_ops) { _e_ops = &e_ops; }
 
-    template <typename VecType, typename MatType, typename ContextType>
+    template <typename ContextType, typename VecType, typename MatType>
     inline void assemble(ContextType   &c,
                          const VecType &X,
                          VecType       *R,
diff --git a/include/mast/base/assembly/libmesh/residual_sensitivity.hpp b/include/mast/base/assembly/libmesh/residual_sensitivity.hpp
index 61eb6c7..2ec4590 100644
--- a/include/mast/base/assembly/libmesh/residual_sensitivity.hpp
+++ b/include/mast/base/assembly/libmesh/residual_sensitivity.hpp
@@ -54,7 +54,10 @@ class ResidualSensitivity {
         
     inline void set_elem_ops(ElemOpsType& e_ops) { _e_ops = &e_ops; }
 
-    template <typename VecType, typename MatType, typename ContextType, typename ScalarFieldType>
+    template <typename ContextType,
+              typename VecType,
+              typename MatType,
+              typename ScalarFieldType>
     inline void assemble(ContextType   &c,
                          const ScalarFieldType& f,
                          const VecType &X,
diff --git a/include/mast/base/material_point_data_storage.hpp b/include/mast/base/material_point_data_storage.hpp
new file mode 100644
index 0000000..e66b9be
--- /dev/null
+++ b/include/mast/base/material_point_data_storage.hpp
@@ -0,0 +1,131 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_material_point_data_storage_h__
+#define __mast_material_point_data_storage_h__
+
+
+// libMesh includes
+#include <libmesh/parallel.h>
+
+
+namespace MAST {
+namespace Base {
+
+template <typename ScalarType>
+class MaterialPointDataStorage {
+    
+public:
+    
+    using view_t       = Eigen::Map<typename Eigen::Matrix<ScalarType, Eigen::Dynamic, 1>>;
+    using const_view_t = Eigen::Map<const typename Eigen::Matrix<ScalarType, Eigen::Dynamic, 1>>;
+    
+    
+    MaterialPointDataStorage(libMesh::Parallel::Communicator &comm,
+                             uint_t                          n_comp,
+                             const std::string               &nm):
+    _comm          (comm),
+    _n_comp        (n_comp),
+    _begin_elem_id (-1),
+    _end_elem_id   (-1),
+    _nm            (nm),
+    _data          (nullptr) {
+        
+    }
+    
+    
+    virtual ~MaterialPointDataStorage() {
+        
+        if (_data) delete _data;
+    }
+    
+    
+    inline void init(uint_t n_pts,
+                     uint_t begin_elem_id,
+                     uint_t end_elem_id) {
+        
+        Assert0(!_data, "Data must be cleared before reinitialization");
+        
+        
+        _begin_elem_id  = begin_elem_id;
+        _end_elem_id    = end_elem_id;
+        
+        _data = new ScalarType[n_pts*_n_comp];
+    }
+    
+    
+    inline view_t
+    value(const uint_t e_id, const uint_t qp) {
+
+        Assert2(e_id >= _begin_elem_id,
+                e_id, _begin_elem_id,
+                "Material point id out of bounds for this rank");
+        Assert2(e_id < _end_elem_id,
+                e_id, _end_elem_id,
+                "Material point id out of bounds for this rank");
+
+        uint_t
+        begin_index = (e_id-_begin_elem_id) * _n_comp;
+        
+        return view_t(&(_data[begin_index]) , _n_comp, 1);
+    }
+
+    
+    inline const_view_t
+    value(const uint_t e_id, const uint_t qp) const {
+        
+        Assert2(e_id >= _begin_elem_id,
+                e_id, _begin_elem_id,
+                "Material point id out of bounds for this rank");
+        Assert2(e_id < _end_elem_id,
+                e_id, _end_elem_id,
+                "Material point id out of bounds for this rank");
+
+        uint_t
+        begin_index = (e_id - _begin_elem_id) * _n_comp;
+        
+        return const_view_t(&(_data[begin_index]) , _n_comp, 1);
+    }
+
+    
+private:
+    
+    libMesh::Parallel::Communicator& _comm;
+    
+    const uint_t                     _n_comp;
+    
+    /*!
+     *  first element stored on this processor
+     */
+    uint_t                           _begin_elem_id;
+    
+    /*!
+     *  one past the last element stored on this processor
+     */
+    uint_t                           _end_elem_id;
+    
+    const std::string                _nm;
+    
+    ScalarType                      *_data;
+};
+
+} // namespace Base
+} // namespace MAST
+
+#endif // __mast_material_point_data_storage_h__
diff --git a/include/mast/base/material_point_system.hpp b/include/mast/base/material_point_system.hpp
new file mode 100644
index 0000000..1570e9a
--- /dev/null
+++ b/include/mast/base/material_point_system.hpp
@@ -0,0 +1,119 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_material_point_system_h__
+#define __mast_material_point_system_h__
+
+// MAST includes
+#include <mast/base/mast_data_types.h>
+
+// libMesh includes
+#include <libmesh/mesh_base.h>
+
+namespace MAST {
+namespace Base {
+
+class MaterialPointSystem {
+  
+public:
+    
+    MaterialPointSystem(libMesh::MeshBase& mesh):
+    _mesh         (mesh),
+    _initialized  (false)
+    { }
+    
+    virtual ~MaterialPointSystem() {
+        
+        // delete the material points
+        /*{
+            std::vector<MAST::Base::MaterialPoint*>::iterator
+            it   = _material_points.begin(),
+            end  = _material_points.end();
+            
+            for (; it != end; it++) delete *it;
+        }*/
+        
+        /*// delete the data
+        {
+            std::map<std::string, MAST::Base::MaterialPointDataStorage<ScalarType>*>::iterator
+            it   = _data.begin(),
+            end  = _data.end();
+            
+            for (; it != end; it++) delete it->second;
+        }*/
+    }
+
+    void initialize(uint_t n_qp_per_elem) {
+
+        Assert0(!_initialized, "System already initialization");
+        
+        _n_points_on_rank.resize(_mesh.comm().size(), 0);
+        _begin_point_id.resize(_mesh.comm().size(), 0);
+
+        uint_t
+        prev_end_id = 0;
+        
+        for (uint_t i=0; i<_mesh.comm().size(); i++) {
+            
+            _n_points_on_rank[i] = _mesh.n_elem_on_proc(i) * n_qp_per_elem;
+            _end_point_id[i]     = prev_end_id + _n_points_on_rank[i];
+            prev_end_id          = _end_point_id[i];
+        }
+        
+        for (uint_t i=1; i<_mesh.comm().size(); i++)
+            _begin_point_id[i]   = _end_point_id[i-1];
+    }
+    
+    
+    /*inline MAST::Base::MaterialPointDataStorage<ScalarType>*>&
+    add_variable(const std::string& nm, uint_t n_components)  {
+        
+        Assert0(!initialized, "Variables can be added prior to system initialization");
+
+        
+    }
+    
+    
+    inline MAST::Base::MaterialPointDataStorage<ScalarType>*>&
+    get_variable(const std::string& nm)  {
+        
+    }*/
+
+
+private:
+
+    libMesh::MeshBase                       &_mesh;
+
+    bool                                     _initialized;
+    
+    std::vector<uint_t>                      _n_points_on_rank;
+    
+    std::vector<uint_t>                      _begin_point_id;
+
+    std::vector<uint_t>                      _end_point_id;
+
+    //std::vector<MAST::Base::MaterialPoint*>  _material_points;
+    
+    //std::map<std::string, MAST::Base::MaterialPointDataStorage<ScalarType>*> _data;
+};
+
+} // namespace Base
+} // namespace MAST
+
+#endif // __mast_material_point_system_h__
diff --git a/include/mast/numerics/utility.hpp b/include/mast/numerics/utility.hpp
index 1ddeed8..800ef95 100644
--- a/include/mast/numerics/utility.hpp
+++ b/include/mast/numerics/utility.hpp
@@ -37,6 +37,33 @@ namespace Numerics {
 namespace Utility {
 
 
+template <typename VecType>
+inline typename
+std::enable_if<std::is_same<typename Eigen::internal::traits<VecType>::Scalar,
+                            real_t>::value, real_t>::type
+real_norm(const VecType& v) {
+    return v.norm();
+}
+
+
+template <typename VecType>
+inline typename
+std::enable_if<std::is_same<typename Eigen::internal::traits<VecType>::Scalar,
+                            complex_t>::value, real_t>::type
+real_norm(const VecType& v) {
+    return v.norm();
+}
+
+
+template <typename VecType>
+inline typename
+std::enable_if<std::is_same<typename Eigen::internal::traits<VecType>::Scalar,
+                            adouble_tl_t>::value, real_t>::type
+real_norm(const VecType& v) {
+    return v.norm().getValue();
+}
+
+
 template <typename ValType>
 inline void
 setZero(ValType& m) { m.setZero();}
diff --git a/include/mast/optimization/topology/simp/libmesh/assemble_output_sensitivity.hpp b/include/mast/optimization/topology/simp/libmesh/assemble_output_sensitivity.hpp
index 63af4db..aad937f 100644
--- a/include/mast/optimization/topology/simp/libmesh/assemble_output_sensitivity.hpp
+++ b/include/mast/optimization/topology/simp/libmesh/assemble_output_sensitivity.hpp
@@ -72,9 +72,9 @@ class AssembleOutputSensitivity {
      *  output derivative is defined as a
      * \f[ \frac{dQ}{d\alpha} = \frac{\partial Q}{\partial \alpha} + \lambda^T \frac{\partial R}{\partial \alpha} \f]
      */
-    template <typename Vec1Type,
+    template <typename ContextType,
+              typename Vec1Type,
               typename Vec2Type,
-              typename ContextType,
               typename FilterType>
     inline void assemble(ContextType               &c,
                          const Vec1Type            &X,
diff --git a/include/mast/optimization/topology/simp/libmesh/residual_and_jacobian.hpp b/include/mast/optimization/topology/simp/libmesh/residual_and_jacobian.hpp
index bc71e1d..8c13c01 100644
--- a/include/mast/optimization/topology/simp/libmesh/residual_and_jacobian.hpp
+++ b/include/mast/optimization/topology/simp/libmesh/residual_and_jacobian.hpp
@@ -57,10 +57,10 @@ class ResidualAndJacobian {
     
     inline void set_elem_ops(ElemOpsType& e_ops) { _e_ops = &e_ops; }
 
-    template <typename Vec1Type,
+    template <typename ContextType,
+              typename Vec1Type,
               typename Vec2Type,
-              typename MatType,
-              typename ContextType>
+              typename MatType>
     inline void assemble(ContextType    &c,
                          const Vec1Type &X,
                          const Vec2Type &density,
diff --git a/include/mast/optimization/topology/simp/libmesh/volume.hpp b/include/mast/optimization/topology/simp/libmesh/volume.hpp
index a8ac552..4c12c17 100644
--- a/include/mast/optimization/topology/simp/libmesh/volume.hpp
+++ b/include/mast/optimization/topology/simp/libmesh/volume.hpp
@@ -45,7 +45,7 @@ class Volume {
     Volume() {}
     virtual ~Volume() {}
     
-    template <typename VecType, typename ContextType>
+    template <typename ContextType, typename VecType>
     inline ScalarType compute(ContextType& c,
                               const VecType &density) const {
         
diff --git a/include/mast/physics/elasticity/continuum_stress.hpp b/include/mast/physics/elasticity/continuum_stress.hpp
index 5a03651..4c30756 100644
--- a/include/mast/physics/elasticity/continuum_stress.hpp
+++ b/include/mast/physics/elasticity/continuum_stress.hpp
@@ -31,8 +31,7 @@ namespace LinearContinuum {
 
 template <typename FEVarType,
           typename MaterialPropertyType,
-          uint_t Dim,
-          typename ContextType>
+          uint_t Dim>
 class Stress {
 
 public:
@@ -72,7 +71,7 @@ class Stress {
     }
     
 
-    template <typename AccessorType>
+    template <typename ContextType, typename AccessorType>
     inline void
     compute(ContextType      &c,
             AccessorType     &stress) const {
@@ -112,7 +111,7 @@ class Stress {
     }
     
     
-    template <typename AccessorType, typename ScalarFieldType>
+    template <typename ContextType, typename AccessorType, typename ScalarFieldType>
     inline void derivative(ContextType            &c,
                            const ScalarFieldType  &f,
                            AccessorType           &dstress) const {
@@ -152,7 +151,7 @@ class Stress {
     }
 
     
-    template <typename AccessorType>
+    template <typename ContextType, typename AccessorType>
     inline void adjoint_derivative(ContextType      &c,
                                    AccessorType     &stress_adj) const {
         
diff --git a/include/mast/physics/elasticity/isotropic_stiffness.hpp b/include/mast/physics/elasticity/isotropic_stiffness.hpp
index d1f3369..24e2236 100644
--- a/include/mast/physics/elasticity/isotropic_stiffness.hpp
+++ b/include/mast/physics/elasticity/isotropic_stiffness.hpp
@@ -39,13 +39,15 @@ shear_modulus_derivative(ScalarType  E, ScalarType nu,
 { return dE/2./(1.+nu) - E/2./pow(1.+nu,2) * dnu;}
 
 
-template <typename ScalarType, uint_t Dim, typename ModulusType, typename PoissonType, typename ContextType>
+template <typename ScalarType, uint_t Dim, typename ModulusType, typename PoissonType>
 class IsotropicMaterialStiffness;
 
-template <typename ScalarType, typename ModulusType, typename PoissonType, typename ContextType>
-class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType, ContextType> {
+template <typename ScalarType, typename ModulusType, typename PoissonType>
+class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType> {
 public:
     
+    static const
+    uint_t dim        = 2;
     using E_scalar_t  = typename ModulusType::scalar_t;
     using nu_scalar_t = typename PoissonType::scalar_t;
     using scalar_t    = typename MAST::DeducedScalarType<typename MAST::DeducedScalarType<E_scalar_t, nu_scalar_t>::type, ScalarType>::type;
@@ -68,12 +70,13 @@ class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType, Contex
     
     inline const PoissonType& get_nu() const {return *_nu;}
     
+    template <typename ContextType>
     inline scalar_t G(ContextType& c) const {
         
         return shear_modulus<scalar_t>(_E->value(c), _nu->value(c));
     }
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline scalar_t G_derivative(ContextType& c, const ScalarFieldType& f) const {
         
         return shear_modulus_derivative<scalar_t>(_E->value(c),
@@ -82,6 +85,7 @@ class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType, Contex
                                                   _nu->derivative(c, f));
     }
     
+    template <typename ContextType>
     inline void value(ContextType& c, value_t& m) const {
                 
         Assert0(_E && _nu, "Material values not provided");
@@ -99,7 +103,7 @@ class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType, Contex
     }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(ContextType&           c,
                            const ScalarFieldType& f,
                            value_t&               m) const {
@@ -133,10 +137,12 @@ class IsotropicMaterialStiffness<ScalarType, 2, ModulusType, PoissonType, Contex
 
 
 
-template <typename ScalarType, typename ModulusType, typename PoissonType, typename ContextType>
-class IsotropicMaterialStiffness<ScalarType, 3, ModulusType, PoissonType, ContextType> {
+template <typename ScalarType, typename ModulusType, typename PoissonType>
+class IsotropicMaterialStiffness<ScalarType, 3, ModulusType, PoissonType> {
 public:
     
+    static const
+    uint_t dim        = 3;
     using E_scalar_t  = typename ModulusType::scalar_t;
     using nu_scalar_t = typename PoissonType::scalar_t;
     using scalar_t    = typename MAST::DeducedScalarType<typename MAST::DeducedScalarType<E_scalar_t, nu_scalar_t>::type, ScalarType>::type;
@@ -154,6 +160,7 @@ class IsotropicMaterialStiffness<ScalarType, 3, ModulusType, PoissonType, Contex
         _nu = &nu;
     }
 
+    template <typename ContextType>
     inline void value(const ContextType& c, value_t& m) const {
                 
         Assert0(_E && _nu, "Material values not provided");
@@ -171,7 +178,7 @@ class IsotropicMaterialStiffness<ScalarType, 3, ModulusType, PoissonType, Contex
     }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(const ContextType&     c,
                            const ScalarFieldType& f,
                            value_t&               m) const {
diff --git a/include/mast/physics/elasticity/linear_strain_energy.hpp b/include/mast/physics/elasticity/linear_strain_energy.hpp
index 657a50f..d9664e7 100644
--- a/include/mast/physics/elasticity/linear_strain_energy.hpp
+++ b/include/mast/physics/elasticity/linear_strain_energy.hpp
@@ -31,8 +31,7 @@ namespace LinearContinuum {
 
 template <typename FEVarType,
           typename SectionPropertyType,
-          uint_t Dim,
-          typename ContextType>
+          uint_t Dim>
 class StrainEnergy {
     
 public:
@@ -73,6 +72,7 @@ class StrainEnergy {
         return Dim*_fe_var_data->get_fe_shape_data().n_basis();
     }
     
+    template <typename ContextType>
     inline void compute(ContextType& c,
                         vector_t& res,
                         matrix_t* jac = nullptr) const {
@@ -121,7 +121,7 @@ class StrainEnergy {
         }
     }
 
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(ContextType& c,
                            const ScalarFieldType& f,
                            vector_t& res,
diff --git a/include/mast/physics/elasticity/material_nonlinear_continuum_strain_energy.hpp b/include/mast/physics/elasticity/material_nonlinear_continuum_strain_energy.hpp
new file mode 100644
index 0000000..71dae26
--- /dev/null
+++ b/include/mast/physics/elasticity/material_nonlinear_continuum_strain_energy.hpp
@@ -0,0 +1,189 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_material_nonlinear_continuum_strain_energy_h__
+#define __mast_material_nonlinear_continuum_strain_energy_h__
+
+// MAST includes
+#include <mast/physics/elasticity/linear_elastic_strain_operator.hpp>
+
+namespace MAST {
+namespace Physics {
+namespace Elasticity {
+namespace ElastoPlasticity {
+
+
+template <typename FEVarType,
+          typename YieldSurfaceType,
+          uint_t Dim>
+class StrainEnergy {
+    
+public:
+
+    using scalar_t         = typename FEVarType::scalar_t;
+    using basis_scalar_t   = typename FEVarType::fe_shape_deriv_t::scalar_t;
+    using vector_t         = typename Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
+    using matrix_t         = typename Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
+    using fe_shape_deriv_t = typename FEVarType::fe_shape_deriv_t;
+    static const uint_t
+    n_strain               = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value;
+
+    StrainEnergy():
+    _yield       (nullptr),
+    _fe_var_data (nullptr)
+    { }
+    
+    virtual ~StrainEnergy() { }
+
+    inline void
+    set_yield_surface(const YieldSurfaceType& ys) {
+        
+        Assert0(!_yield, "Yield surface already initialized.");
+        
+        _yield = &ys;
+    }
+
+    inline void set_fe_var_data(const FEVarType& fe_data)
+    {
+        Assert0(!_fe_var_data, "FE data already initialized.");
+        _fe_var_data = &fe_data;
+    }
+
+    inline uint_t n_dofs() const {
+
+        Assert0(_fe_var_data, "FE data not initialized.");
+
+        return Dim*_fe_var_data->get_fe_shape_data().n_basis();
+    }
+    
+    template <typename ContextType>
+    inline void compute(ContextType& c,
+                        vector_t& res,
+                        matrix_t* jac = nullptr) const {
+        
+        Assert0(_fe_var_data, "FE data not initialized.");
+        Assert0(_yield, "Yield surface not initialized");
+        
+        const typename FEVarType::fe_shape_deriv_t
+        &fe = _fe_var_data->get_fe_shape_data();
+        
+        typename Eigen::Matrix<scalar_t, n_strain, 1>
+        epsilon,
+        stress;
+        vector_t
+        vec     = vector_t::Zero(Dim*fe.n_basis());
+        
+        typename YieldSurfaceType::stiff_t
+        mat;
+        
+        matrix_t
+        mat1 = matrix_t::Zero(n_strain, Dim*fe.n_basis()),
+        mat2 = matrix_t::Zero(Dim*fe.n_basis(), Dim*fe.n_basis());
+
+        MAST::Numerics::FEMOperatorMatrix<scalar_t>
+        Bxmat;
+        Bxmat.reinit(n_strain, Dim, fe.n_basis());
+
+        
+        for (uint_t i=0; i<fe.n_q_points(); i++) {
+            
+            c.init_for_qp(i);
+            
+            MAST::Physics::Elasticity::LinearContinuum::strain
+            <scalar_t, scalar_t, FEVarType, Dim>(*_fe_var_data, i, epsilon, Bxmat);
+            
+            // evaluate the stress and tangent stiffness matrix
+            _yield->compute(c, epsilon, *c.mp_accessor, &mat);
+
+            Bxmat.vector_mult_transpose(vec, stress);
+            res += fe.detJxW(i) * vec;
+            
+            if (jac) {
+                
+                Bxmat.left_multiply(mat1, mat);
+                Bxmat.right_multiply_transpose(mat2, mat1);
+                (*jac) += fe.detJxW(i) * mat2;
+            }
+        }
+    }
+
+    template <typename ContextType, typename ScalarFieldType>
+    inline void derivative(ContextType& c,
+                           const ScalarFieldType& f,
+                           vector_t& res,
+                           matrix_t* jac = nullptr) const {
+        
+        Assert0(_fe_var_data, "FE data not initialized.");
+        Assert0(_yield, "Yield surface not initialized");
+        
+        const typename FEVarType::fe_shape_deriv_t
+        &fe = _fe_var_data->get_fe_shape_data();
+
+        typename Eigen::Matrix<scalar_t, n_strain, 1>
+        epsilon,
+        stress;
+        vector_t
+        vec     = vector_t::Zero(Dim*fe.n_basis());
+
+        typename YieldSurfaceType::stiff_t
+        mat;
+        matrix_t
+        mat1 = matrix_t::Zero(n_strain, Dim*fe.n_basis()),
+        mat2 = matrix_t::Zero(Dim*fe.n_basis(), Dim*fe.n_basis());
+
+        MAST::Numerics::FEMOperatorMatrix<scalar_t>
+        Bxmat;
+        Bxmat.reinit(n_strain, Dim, fe.n_basis());
+
+        
+        for (uint_t i=0; i<fe.n_q_points(); i++) {
+            
+            c.init_for_qp(i);
+
+            MAST::Physics::Elasticity::LinearContinuum::strain
+            <scalar_t, scalar_t, FEVarType, Dim>(*_fe_var_data, i, epsilon, Bxmat);
+            
+            //_yield->derivative(c, f, stress, &mat);
+
+            Bxmat.vector_mult_transpose(vec, stress);
+            res += fe.detJxW(i) * vec;
+            
+            if (jac) {
+                
+                Bxmat.left_multiply(mat1, mat);
+                Bxmat.right_multiply_transpose(mat2, mat1);
+                (*jac) += fe.detJxW(i) * mat2;
+            }
+        }
+    }
+
+    
+private:
+    
+    
+    const YieldSurfaceType          *_yield;
+    const FEVarType                 *_fe_var_data;
+};
+
+}  // namespace MaterialNonlinear
+}  // namespace Elasticity
+}  // namespace Physics
+}  // namespace MAST
+
+#endif // __mast_material_nonlinear_continuum_strain_energy_h__
diff --git a/include/mast/physics/elasticity/mindlin_plate_strain_energy.hpp b/include/mast/physics/elasticity/mindlin_plate_strain_energy.hpp
index d206f3b..c81db17 100644
--- a/include/mast/physics/elasticity/mindlin_plate_strain_energy.hpp
+++ b/include/mast/physics/elasticity/mindlin_plate_strain_energy.hpp
@@ -30,8 +30,7 @@ namespace MindlinPlate {
 
 
 template <typename FEVarType,
-          typename SectionPropertyType,
-          typename ContextType>
+          typename SectionPropertyType>
 class StrainEnergy {
     
 public:
@@ -78,6 +77,7 @@ class StrainEnergy {
         return 3*_bending_fe_var_data->get_fe_shape_data().n_basis();
     }
     
+    template <typename ContextType>
     inline void compute(ContextType& c,
                         vector_t& res,
                         matrix_t* jac = nullptr) const {
@@ -175,7 +175,7 @@ class StrainEnergy {
         }
     }
 
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(ContextType& c,
                            const ScalarFieldType& f,
                            vector_t& res,
diff --git a/include/mast/physics/elasticity/plate_bending_section_property.hpp b/include/mast/physics/elasticity/plate_bending_section_property.hpp
index a7e6628..b0098e4 100644
--- a/include/mast/physics/elasticity/plate_bending_section_property.hpp
+++ b/include/mast/physics/elasticity/plate_bending_section_property.hpp
@@ -30,8 +30,7 @@ namespace Elasticity {
 
 template <typename ScalarType,
           typename MaterialType,
-          typename ThicknessType,
-          typename ContextType>
+          typename ThicknessType>
 class PlateBendingSectionProperty {
 public:
     
@@ -59,6 +58,7 @@ class PlateBendingSectionProperty {
     
     inline const ThicknessType& get_thickness() const {return *_th;}
     
+    template <typename ContextType>
     inline void inplane_value(ContextType     &c,
                               inplane_value_t &m) const {
 
@@ -69,7 +69,7 @@ class PlateBendingSectionProperty {
      }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void inplane_derivative(ContextType           &c,
                                    const ScalarFieldType &f,
                                    inplane_value_t       &m) const {
@@ -87,6 +87,7 @@ class PlateBendingSectionProperty {
         m  += dm;
     }
 
+    template <typename ContextType>
     inline void shear_value(ContextType& c,
                             shear_value_t &m) const {
 
@@ -97,7 +98,7 @@ class PlateBendingSectionProperty {
      }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void shear_derivative(ContextType           &c,
                                  const ScalarFieldType &f,
                                  shear_value_t         &m) const {
diff --git a/include/mast/physics/elasticity/pressure_load.hpp b/include/mast/physics/elasticity/pressure_load.hpp
index f869c8b..116614f 100644
--- a/include/mast/physics/elasticity/pressure_load.hpp
+++ b/include/mast/physics/elasticity/pressure_load.hpp
@@ -31,8 +31,7 @@ namespace Elasticity {
 template <typename FEVarType,
           typename PressureFieldType,
           typename SectionAreaType,
-          uint_t Dim,
-          typename ContextType>
+          uint_t Dim>
 class SurfacePressureLoad {
 
 public:
@@ -63,6 +62,7 @@ class SurfacePressureLoad {
         return Dim*_fe_var_data->get_fe_shape_data().n_basis();
     }
 
+    template <typename ContextType>
     inline void compute(ContextType& c,
                         vector_t& res,
                         matrix_t* jac = nullptr) const {
@@ -93,7 +93,7 @@ class SurfacePressureLoad {
     }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(ContextType& c,
                            const ScalarFieldType& f,
                            vector_t& res,
diff --git a/include/mast/physics/elasticity/shell_face_pressure_load.hpp b/include/mast/physics/elasticity/shell_face_pressure_load.hpp
index ffa03ec..c0582b1 100644
--- a/include/mast/physics/elasticity/shell_face_pressure_load.hpp
+++ b/include/mast/physics/elasticity/shell_face_pressure_load.hpp
@@ -29,8 +29,7 @@ namespace Physics {
 namespace Elasticity {
 
 template <typename FEVarType,
-          typename PressureFieldType,
-          typename ContextType>
+          typename PressureFieldType>
 class ShellFacePressureLoad {
 
 public:
@@ -77,6 +76,7 @@ class ShellFacePressureLoad {
         return _fe_var_data->get_fe_shape_data().n_basis();
     }
 
+    template <typename ContextType>
     inline void compute(ContextType& c,
                         vector_t& res,
                         matrix_t* jac = nullptr) const {
@@ -98,7 +98,7 @@ class ShellFacePressureLoad {
     }
     
     
-    template <typename ScalarFieldType>
+    template <typename ContextType, typename ScalarFieldType>
     inline void derivative(ContextType& c,
                            const ScalarFieldType& f,
                            vector_t& res,
diff --git a/include/mast/physics/elasticity/von_mises_stress.hpp b/include/mast/physics/elasticity/von_mises_stress.hpp
index cfe76b8..b668101 100644
--- a/include/mast/physics/elasticity/von_mises_stress.hpp
+++ b/include/mast/physics/elasticity/von_mises_stress.hpp
@@ -1,21 +1,21 @@
 /*
-* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
-* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
-*
-* This library is free software; you can redistribute it and/or
-* modify it under the terms of the GNU Lesser General Public
-* License as published by the Free Software Foundation; either
-* version 2.1 of the License, or (at your option) any later version.
-*
-* This library is distributed in the hope that it will be useful,
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-* Lesser General Public License for more details.
-*
-* You should have received a copy of the GNU Lesser General Public
-* License along with this library; if not, write to the Free Software
-* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+ * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+ * Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
 
 #ifndef __mast_linear_continuum_von_mises_stress_h__
 #define __mast_linear_continuum_von_mises_stress_h__
@@ -29,207 +29,434 @@ namespace Elasticity {
 namespace LinearContinuum {
 
 
-template <typename ScalarType, uint_t Dim>
-using stress_vec_t = Eigen::Matrix<ScalarType,
-                                   MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-                                   1>;
 
 template <typename ScalarType, uint_t Dim>
-using stress_adjoint_mat_t = Eigen::Matrix<ScalarType,
-                                   MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-                                   Eigen::Dynamic>;
-
+struct
+vonMisesStress {};
 
 
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==2, ScalarType>::type
-vonMises_stress(stress_vec_t<ScalarType, Dim>& stress) {
-   
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
+template <typename ScalarType>
+struct
+vonMisesStress<ScalarType, 2> {
     
-    return
-    pow(0.5 * (pow(stress(0)-stress(1),2) +     //(((sigma_xx - sigma_yy)^2    +
-               pow(stress(1),2) +               //  (sigma_yy)^2    +
-               pow(stress(0),2)) +              //  (sigma_xx)^2)/2 +
-        3.0 * (pow(stress(2), 2)), 0.5);        //   3.0 * tau_xy^2)^.5
-}
-
-
-
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==3, ScalarType>::type
-vonMises_stress(stress_vec_t<ScalarType, Dim>& stress) {
-   
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
-    
-    return
-    pow(0.5 * (pow(stress(0)-stress(1),2) +    //(((sigma_xx - sigma_yy)^2    +
-               pow(stress(1)-stress(2),2) +    //  (sigma_yy - sigma_zz)^2    +
-               pow(stress(2)-stress(0),2)) +   //  (sigma_zz - sigma_xx)^2)/2 +
-        3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
-               pow(stress(4), 2) +              //     tau_yz^2 +
-               pow(stress(5), 2)), 0.5);        //     tau_zx^2))^.5
-}
+    template <typename VecType>
+    static inline ScalarType
+    value(const VecType& stress) {
+
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<VecType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        
+        return
+        pow(0.5 * (pow(stress(0)-stress(1),2) +     //(((sigma_xx - sigma_yy)^2    +
+                   pow(stress(1),2) +               //  (sigma_yy)^2    +
+                   pow(stress(0),2)) +              //  (sigma_xx)^2)/2 +
+            3.0 * (pow(stress(2), 2)), 0.5);        //   3.0 * tau_xy^2)^.5
+    }
 
+    
+    
+    template <typename Vec1Type, typename Vec2Type>
+    static inline void
+    derivative(const Vec1Type  &stress,
+               Vec2Type        &dstress) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                dstress.size(),
+                "Incorrect stress vector dimension");
+        
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1),2) +               // (sigma_yy)^2    +
+               pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
+        3.0 * (pow(stress(2), 2));              //  3.0 * tau_xy^2)
+        
+        dstress.setZero();
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (abs(p) > 0.) {
+            
+            dstress(0) =   (stress(0) - stress(1)) + stress(0);
+            dstress(1) = - (stress(0) - stress(1)) + stress(1);
+            dstress(2) = 6. * stress(2);
+            
+            dstress *= 0.5 * pow(p, -0.5);
+        }
+    }
 
+    
+    template <typename VecType, typename MatType>
+    static inline void
+    second_derivative(const VecType   &stress,
+                      MatType         &dstress) {
+        
+        using scalar_t = typename Eigen::internal::traits<VecType>::Scalar;
+        const uint_t
+        ns = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value;
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<VecType>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<MatType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                dstress.rows(),
+                "Incorrect derivative matrix row dimension");
+        Assert1(dstress.cols() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                dstress.cols(),
+                "Incorrect derivative matrix column dimension");
+
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1),2) +               // (sigma_yy)^2    +
+               pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
+        3.0 * (pow(stress(2), 2)),              //  3.0 * tau_xy^2)
+        sqrtp_2 = .5 * pow(p, -0.5);
+        
+        dstress.setZero();
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (abs(p) > 0.) {
+            
+            Eigen::Matrix<scalar_t, ns, 1>
+            ds;
+            ds(0) =   (stress(0) - stress(1)) + stress(0);
+            ds(1) = - (stress(0) - stress(1)) + stress(1);
+            ds(2) = 6. * stress(2);
+            
+            dstress  = -.25 / pow(p, 1.5) * ds * ds.transpose();
+            
+            dstress(0, 0) +=  2. * sqrtp_2;
+            dstress(0, 1) += -1. * sqrtp_2;
+            
+            dstress(1, 0) += -1. * sqrtp_2;
+            dstress(1, 1) +=  2. * sqrtp_2;
+
+            dstress(2, 2) +=  6. * sqrtp_2;
+        }
+    }
 
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==2, ScalarType>::type
-vonMises_stress_derivative(stress_vec_t<ScalarType, Dim>& stress,
-                           stress_vec_t<ScalarType, Dim>& dstress_dp) {
     
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
-    Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            dstress_dp.size(),
-            "Incorrect stress vector dimension");
     
-    ScalarType
-    p =
-    0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
-           pow(stress(1),2) +               // (sigma_yy)^2    +
-           pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
-    3.0 * (pow(stress(2), 2)),              //  3.0 * tau_xy^2)
-    dp = 0.;
+    template <typename Vec1Type, typename Vec2Type>
+    static inline ScalarType
+    derivative_sens(const Vec1Type  &stress,
+                    const Vec2Type  &dstress_dp) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                dstress_dp.size(),
+                "Incorrect stress vector dimension");
+        
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1),2) +               // (sigma_yy)^2    +
+               pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
+        3.0 * (pow(stress(2), 2)),              //  3.0 * tau_xy^2)
+        dp = 0.;
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.)
+            dp =
+            (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
+              (dstress_dp(1)) * (stress(1)) +
+              (dstress_dp(0)) * (stress(0))) +
+             6.0 * (dstress_dp(2) * stress(2))) * 0.5 * pow(p, -0.5);
+        
+        return dp;
+    }
     
-    // if p == 0, then the sensitivity returns nan
-    // Hence, we are avoiding this by setting it to zero whenever p = 0.
-    if (fabs(p) > 0.)
-        dp =
-        (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
-          (dstress_dp(1)) * (stress(1)) +
-          (dstress_dp(0)) * (stress(0))) +
-         6.0 * (dstress_dp(2) * stress(2))) * 0.5 * pow(p, -0.5);
     
-    return dp;
-}
-
-
-
-
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==3, ScalarType>::type
-vonMises_stress_derivative(stress_vec_t<ScalarType, Dim>& stress,
-                           stress_vec_t<ScalarType, Dim>& dstress_dp) {
-    
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
-    Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            dstress_dp.size(),
-            "Incorrect stress vector dimension");
-    
-    ScalarType
-    p =
-    0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
-           pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
-           pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
-    3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
-           pow(stress(4), 2) +              //     tau_yz^2 +
-           pow(stress(5), 2)),              //     tau_zx^2)
-    dp = 0.;
-    
-    // if p == 0, then the sensitivity returns nan
-    // Hence, we are avoiding this by setting it to zero whenever p = 0.
-    if (fabs(p) > 0.)
-        dp =
-        (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
-          (dstress_dp(1) - dstress_dp(2)) * (stress(1) - stress(2)) +
-          (dstress_dp(2) - dstress_dp(0)) * (stress(2) - stress(0))) +
-         6.0 * (dstress_dp(3) * stress(3)+
-                dstress_dp(4) * stress(4)+
-                dstress_dp(5) * stress(5))) * 0.5 * pow(p, -0.5);
-    
-    return dp;
-}
-
-
+    
+    template <typename Vec1Type, typename MatType, typename Vec2Type>
+    static inline void
+    stress_dX(const Vec1Type   &stress,
+              const MatType    &dstress_dX,
+              Vec2Type         &vm_adjoint) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<MatType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<2>::value,
+                dstress_dX.rows(),
+                "Incorrect stress adjoint matrix rows");
+        Assert2(vm_adjoint.size() == dstress_dX.cols(),
+                vm_adjoint.size(),   dstress_dX.cols(),
+                "Incorrect von Mises adjoint vector dimension");
+        
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1),2) +               // (sigma_yy)^2    +
+               pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
+        3.0 * (pow(stress(2), 2));              // 3* (tau_xy^2)
+        
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.)
+            vm_adjoint =
+            (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
+              dstress_dX.row(1) * stress(1) +
+              dstress_dX.row(0) * stress(0)) +
+             6.0 * dstress_dX.row(2) * stress(2)) * 0.5 * pow(p, -0.5);
+    }
+    
+    
+};
 
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==2, void>::type
-vonMises_stress_dX(stress_vec_t<ScalarType, Dim>                &stress,
-                   stress_adjoint_mat_t<ScalarType, Dim>        &dstress_dX,
-                   Eigen::Matrix<ScalarType, Eigen::Dynamic, 1> &vm_adjoint) {
-    
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
-    Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            dstress_dX.rows(),
-            "Incorrect stress adjoint matrix rows");
-    Assert2(vm_adjoint.size() == dstress_dX.cols(),
-            vm_adjoint.size(),   dstress_dX.cols(),
-            "Incorrect von Mises adjoint vector dimension");
-
-    ScalarType
-    p =
-    0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
-           pow(stress(1),2) +               // (sigma_yy)^2    +
-           pow(stress(0),2)) +              // (sigma_xx)^2)/2 +
-    3.0 * (pow(stress(2), 2));              // 3* (tau_xy^2)
-
-    
-    // if p == 0, then the sensitivity returns nan
-    // Hence, we are avoiding this by setting it to zero whenever p = 0.
-    if (fabs(p) > 0.)
-        vm_adjoint =
-        (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
-          dstress_dX.row(1) * stress(1) +
-          dstress_dX.row(0) * stress(0)) +
-         6.0 * dstress_dX.row(2) * stress(2)) * 0.5 * pow(p, -0.5);
-}
 
+template <typename ScalarType>
+struct
+vonMisesStress<ScalarType, 3> {
+    
+    template <typename VecType>
+    static inline ScalarType
+    value(const VecType& stress) {
+
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<VecType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        
+
+        return
+        pow(0.5 * (pow(stress(0)-stress(1),2) +    //(((sigma_xx - sigma_yy)^2    +
+                   pow(stress(1)-stress(2),2) +    //  (sigma_yy - sigma_zz)^2    +
+                   pow(stress(2)-stress(0),2)) +   //  (sigma_zz - sigma_xx)^2)/2 +
+            3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
+                   pow(stress(4), 2) +              //     tau_yz^2 +
+                   pow(stress(5), 2)), 0.5);        //     tau_zx^2))^.5
+    }
+    
+    
+    /*!
+     * computes the derivative of
+     */
+    template <typename Vec1Type, typename Vec2Type>
+    static inline void
+    derivative(const Vec1Type  &stress,
+               Vec2Type        &dstress) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                dstress.size(),
+                "Incorrect stress vector dimension");
+
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
+               pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
+        3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
+               pow(stress(4), 2) +              //     tau_yz^2 +
+               pow(stress(5), 2));              //     tau_zx^2)
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.) {
+            
+            dstress(0) =   (stress(0) - stress(1)) - (stress(2) - stress(0));
+            dstress(1) = - (stress(0) - stress(1)) + (stress(1) - stress(2));
+            dstress(2) = - (stress(1) - stress(2)) + (stress(2) - stress(0));
+            dstress(3) = 6. * stress(3);
+            dstress(4) = 6. * stress(4);
+            dstress(5) = 6. * stress(5);
+            
+            dstress *= 0.5 * pow(p, -0.5);
+        }
+    }
 
+    
+    
+    template <typename VecType, typename MatType>
+    static inline void
+    second_derivative(const VecType   &stress,
+                      MatType         &dstress) {
+        
+        using scalar_t = typename Eigen::internal::traits<VecType>::Scalar;
+        const uint_t
+        ns = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value;
+
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<VecType>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<MatType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                dstress.rows(),
+                "Incorrect derivative matrix row dimension");
+        Assert1(dstress.cols() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                dstress.cols(),
+                "Incorrect derivative matrix column dimension");
+
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
+               pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
+        3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
+               pow(stress(4), 2) +              //     tau_yz^2 +
+               pow(stress(5), 2)),              //     tau_zx^2)
+        sqrtp_2 = .5 * pow(p, -0.5);
+
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.) {
+
+            Eigen::Matrix<scalar_t, ns, 1>
+            ds;
+            ds(0) =   (stress(0) - stress(1)) - (stress(2) - stress(0));
+            ds(1) = - (stress(0) - stress(1)) + (stress(1) - stress(2));
+            ds(2) = - (stress(1) - stress(2)) + (stress(2) - stress(0));
+            ds(3) = 6. * stress(3);
+            ds(4) = 6. * stress(4);
+            ds(5) = 6. * stress(5);
+
+            dstress  = -.25 / pow(p, 1.5) * ds * ds.transpose();
+            
+            dstress(0, 0) +=   2. * sqrtp_2;
+            dstress(0, 1) +=  -1. * sqrtp_2;
+            dstress(0, 2) +=  -1. * sqrtp_2;
+
+            dstress(1, 0) +=  -1. * sqrtp_2;
+            dstress(1, 1) +=   2. * sqrtp_2;
+            dstress(1, 2) +=  -1. * sqrtp_2;
+
+            dstress(2, 0) +=  -1. * sqrtp_2;
+            dstress(2, 1) +=  -1. * sqrtp_2;
+            dstress(2, 2) +=   2. * sqrtp_2;
+
+            dstress(3, 3) += 6. * sqrtp_2;
+            dstress(4, 4) += 6. * sqrtp_2;
+            dstress(5, 5) += 6. * sqrtp_2;
+        }
+    }
 
-template <typename ScalarType, uint_t Dim>
-inline
-typename std::enable_if<Dim==3, void>::type
-vonMises_stress_dX(stress_vec_t<ScalarType, Dim>                &stress,
-                   stress_adjoint_mat_t<ScalarType, Dim>        &dstress_dX,
-                   Eigen::Matrix<ScalarType, Eigen::Dynamic, 1> &vm_adjoint) {
-    
-    Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            stress.size(),
-            "Incorrect stress vector dimension");
-    Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
-            dstress_dX.rows(),
-            "Incorrect stress adjoint matrix rows");
-    Assert2(vm_adjoint.size() == dstress_dX.cols(),
-            vm_adjoint.size(),   dstress_dX.cols(),
-            "Incorrect von Mises adjoint vector dimension");
-
-    ScalarType
-    p =
-    0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
-           pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
-           pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
-    3.0 * (pow(stress(3), 2) +              // 3* (tau_xx^2 +
-           pow(stress(4), 2) +              //     tau_yy^2 +
-           pow(stress(5), 2));              //     tau_zz^2)
-
-    
-    // if p == 0, then the sensitivity returns nan
-    // Hence, we are avoiding this by setting it to zero whenever p = 0.
-    if (fabs(p) > 0.)
-        vm_adjoint =
-        (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
-          (dstress_dX.row(1) - dstress_dX.row(2)) * (stress(1) - stress(2)) +
-          (dstress_dX.row(2) - dstress_dX.row(0)) * (stress(2) - stress(0))) +
-         6.0 * (dstress_dX.row(3) * stress(3)+
-                dstress_dX.row(4) * stress(4)+
-                dstress_dX.row(5) * stress(5))) * 0.5 * pow(p, -0.5);
-}
+    
+    template <typename Vec1Type, typename Vec2Type>
+    static inline ScalarType
+    derivative_sens(const Vec1Type  &stress,
+                    const Vec2Type  &dstress_dp) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                dstress_dp.size(),
+                "Incorrect stress vector dimension");
+
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
+               pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
+        3.0 * (pow(stress(3), 2) +              // 3* (tau_xy^2 +
+               pow(stress(4), 2) +              //     tau_yz^2 +
+               pow(stress(5), 2)),              //     tau_zx^2)
+        dp = 0.;
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.)
+            dp =
+            (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
+              (dstress_dp(1) - dstress_dp(2)) * (stress(1) - stress(2)) +
+              (dstress_dp(2) - dstress_dp(0)) * (stress(2) - stress(0))) +
+             6.0 * (dstress_dp(3) * stress(3)+
+                    dstress_dp(4) * stress(4)+
+                    dstress_dp(5) * stress(5))) * 0.5 * pow(p, -0.5);
+        
+        return dp;
+    }
 
+    
+    template <typename Vec1Type, typename MatType, typename Vec2Type>
+    static inline void
+    stress_dX(const Vec1Type   &stress,
+              const MatType    &dstress_dX,
+              Vec2Type         &vm_adjoint) {
+        
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec1Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<Vec2Type>::Scalar>::value,
+                      "Scalar type must be same");
+        static_assert(std::is_same<ScalarType, typename Eigen::internal::traits<MatType>::Scalar>::value,
+                      "Scalar type must be same");
+        Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                stress.size(),
+                "Incorrect stress vector dimension");
+        Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<3>::value,
+                dstress_dX.rows(),
+                "Incorrect stress adjoint matrix rows");
+        Assert2(vm_adjoint.size() == dstress_dX.cols(),
+                vm_adjoint.size(),   dstress_dX.cols(),
+                "Incorrect von Mises adjoint vector dimension");
+
+        ScalarType
+        p =
+        0.5 * (pow(stress(0)-stress(1),2) +    //((sigma_xx - sigma_yy)^2    +
+               pow(stress(1)-stress(2),2) +    // (sigma_yy - sigma_zz)^2    +
+               pow(stress(2)-stress(0),2)) +   // (sigma_zz - sigma_xx)^2)/2 +
+        3.0 * (pow(stress(3), 2) +              // 3* (tau_xx^2 +
+               pow(stress(4), 2) +              //     tau_yy^2 +
+               pow(stress(5), 2));              //     tau_zz^2)
+        
+        
+        // if p == 0, then the sensitivity returns nan
+        // Hence, we are avoiding this by setting it to zero whenever p = 0.
+        if (fabs(p) > 0.)
+            vm_adjoint =
+            (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
+              (dstress_dX.row(1) - dstress_dX.row(2)) * (stress(1) - stress(2)) +
+              (dstress_dX.row(2) - dstress_dX.row(0)) * (stress(2) - stress(0))) +
+             6.0 * (dstress_dX.row(3) * stress(3)+
+                    dstress_dX.row(4) * stress(4)+
+                    dstress_dX.row(5) * stress(5))) * 0.5 * pow(p, -0.5);
+    }
+};
 
 }  // namespace LinearContinuum
 }  // namespace Elasticity
diff --git a/include/mast/physics/elasticity/von_mises_yield_criterion.hpp b/include/mast/physics/elasticity/von_mises_yield_criterion.hpp
new file mode 100644
index 0000000..a3186c7
--- /dev/null
+++ b/include/mast/physics/elasticity/von_mises_yield_criterion.hpp
@@ -0,0 +1,502 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_elastoplasticity_von_mises_yield_criterion_h__
+#define __mast_elastoplasticity_von_mises_yield_criterion_h__
+
+// MAST includes
+#include <mast/base/mast_data_types.h>
+#include <mast/physics/elasticity/linear_elastic_strain_operator.hpp>
+#include <mast/physics/elasticity/von_mises_stress.hpp>
+#include <mast/numerics/utility.hpp>
+
+
+namespace MAST {
+namespace Physics {
+namespace Elasticity {
+namespace ElastoPlasticity {
+
+template <typename ScalarType,
+          typename YieldCriterionType>
+class Accessor {
+  
+public:
+    
+    using stress_t       =
+    Eigen::Map<Eigen::Matrix<ScalarType, YieldCriterionType::n_strain, 1>>;
+    using const_stress_t =
+    Eigen::Map<const Eigen::Matrix<ScalarType, YieldCriterionType::n_strain, 1>>;
+
+    Accessor():
+    _n_dofs  (-1),
+    _data    (nullptr),
+    _yield   (nullptr) { }
+    
+    
+    virtual ~Accessor() { }
+
+
+    /*!
+     * initializes the data in this accessor with a specified data independent of a system
+     */
+    inline void init(YieldCriterionType   &yield,
+                     ScalarType           *data) {
+        
+        _yield  = &yield;
+        _data   = data;
+        _n_dofs = yield.n_variables();
+    }
+    
+    
+    inline stress_t stress() {
+        
+        return stress_t(_data, YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline const_stress_t stress() const {
+        
+        return const_stress_t(_data, YieldCriterionType::n_strain, 1);
+    }
+
+
+    inline stress_t strain() {
+        
+        return stress_t(&(_data[YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline const_stress_t strain() const {
+        
+        return const_stress_t(&(_data[YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+
+    inline stress_t plastic_strain() {
+        
+        return stress_t(&(_data[2*YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline const_stress_t plastic_strain() const {
+        
+        return const_stress_t(&(_data[2*YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline stress_t back_stress() {
+        
+        Assert0(_yield->if_kinematic_hardening,
+                "Kinematic hardening not enabled for yield criterion");
+        
+        return stress_t(&(_data[3*YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline const_stress_t back_stress() const {
+        
+        Assert0(_yield->if_kinematic_hardening,
+                "Kinematic hardening not enabled for yield criterion");
+        
+        return const_stress_t(&(_data[3*YieldCriterionType::n_strain]), YieldCriterionType::n_strain, 1);
+    }
+
+    
+    inline ScalarType consistency_parameter() const {
+        
+        return _data[_n_dofs-1];
+    }
+
+    inline ScalarType& consistency_parameter() {
+        
+        return _data[_n_dofs-1];
+    }
+
+    
+private:
+   
+    uint_t               _n_dofs;
+    ScalarType          *_data;
+    YieldCriterionType  *_yield;
+};
+
+
+
+template <typename ScalarType,
+          typename MaterialType>
+class vonMisesYieldFunction {
+    
+public:
+    
+    using scalar_t   = typename MaterialType::scalar_t;
+    static_assert(std::is_same<scalar_t, typename MaterialType::scalar_t>::value,
+                  "Scalar type should be same for yield function and material stiffness");
+    using material_t = MaterialType;
+    static const
+    uint_t dim       = MaterialType::dim;
+    static const
+    uint_t n_strain  = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<dim>::value;
+    using stiff_t    = typename MaterialType::value_t;
+    
+    
+    vonMisesYieldFunction():
+    _if_kinematic_hardening (false),
+    _vm_lim                 (0.),
+    _material               (nullptr) {
+        
+    }
+    
+    
+    virtual ~vonMisesYieldFunction() { }
+    
+    
+    
+    inline void set_limit_stress(const real_t v) {
+        
+        _vm_lim = v;
+    }
+
+
+    inline void set_material(material_t  &m) {
+        
+        _material = &m;
+    }
+
+    
+    inline bool if_kinematic_hardening() const {
+        
+        return _if_kinematic_hardening;
+    }
+
+    
+    inline uint_t n_variables() const {
+
+        uint_t
+        n  = 3*n_strain+1; // stress, strain, plastic strain and consistency parameter
+        
+        // if kinematic hardening is enabled then a backstress is included
+        n += _if_kinematic_hardening?n_strain:0;
+            
+        return n;
+    }
+
+
+    
+    template <typename ContextType,
+              typename VecType,
+              typename AccessorType>
+    inline void compute(ContextType     &c,
+                        const VecType   &strain,
+                        AccessorType    &internal,
+                        stiff_t         *stiff = nullptr) const {
+        
+        Assert2(strain.size() == n_strain,
+                strain.size(), n_strain,
+                "Incompatible strain vector dimension");
+        Assert2(strain.size() == n_strain,
+                strain.size(), n_strain,
+                "Incompatible strain vector dimension");
+
+        
+        Eigen::Matrix<scalar_t, n_strain, 1>
+        strain_e  = Eigen::Matrix<scalar_t, n_strain, 1>::Zero(n_strain);
+        
+        // elastic part of the strain
+        strain_e = strain - internal.plastic_strain();
+        
+        stiff_t
+        m_stiff = stiff_t::Zero();
+        
+        // material stiffness matrix
+        _material->value(c, m_stiff);
+        
+        // store the strain in the accessor
+        internal.strain() = strain;
+        
+        // stress due to elastic strain
+        internal.stress() = m_stiff * strain_e;
+        
+        // evaluate the value of yield criterion
+        scalar_t
+        yield = this->yield_function(c, internal);
+        
+        // if the yield function is violated then an update needs to be computed
+        if (_active_yield_function(yield)) {
+            
+            Eigen::Matrix<scalar_t, n_strain+1, 1>
+            dx   = Eigen::Matrix<scalar_t, n_strain+1, 1>::Zero(),
+            res  = Eigen::Matrix<scalar_t, n_strain+1, 1>::Zero();
+            
+            Eigen::Matrix<scalar_t, n_strain+1, n_strain+1>
+            jac    = Eigen::Matrix<scalar_t, n_strain+1, n_strain+1>::Zero();
+
+            bool
+            terminate = false;
+
+            real_t
+            res_norm        = 0.,
+            dx_norm         = 0.,
+            res_scaled_norm = 0.,
+            stress_norm     = 0.,
+            tol             = 1.e-10;
+            
+            uint_t
+            iter     = 0,
+            max_it   = 20;
+            
+            while (!terminate) {
+                
+                this->return_mapping_residual_and_jacobian(c, strain, internal, res, &jac);
+                
+                // check the norm and if we need to continue iteration
+                res_scaled_norm  = res_norm = MAST::Numerics::Utility::real_norm(res);
+                stress_norm      = MAST::Numerics::Utility::real_norm(internal.stress());
+                
+                if (stress_norm > 0.)
+                    res_scaled_norm /= stress_norm;
+                
+                std::cout
+                << "Iter: " << std::setw(5) << iter
+                << " || res ||_2 = "
+                << std::setw(20) << res_norm
+                << ", || res ||_2/|| sigma ||_2 = "
+                << std::setw(20) << res_scaled_norm
+                << ", || dx ||_2 = "
+                << std::setw(20) << dx_norm
+                << std::endl;
+
+                if (!terminate) {
+                    
+                    ////////////////////////////
+                    // update to the solution
+                    ////////////////////////////
+                    dx = Eigen::FullPivLU<Eigen::Matrix<scalar_t, n_strain+1, n_strain+1>>(jac).solve(res);
+                    
+                    internal.stress() -= dx.topRows(n_strain);
+                    internal.consistency_parameter() -= dx(n_strain);
+                    
+                    // increment the iteration
+                    iter++;
+                    
+                    dx_norm = MAST::Numerics::Utility::real_norm(dx);
+                    
+                    if (dx_norm < tol) {
+                        
+                        terminate = true;
+                        std::cout
+                        << "Terminating Return-Mapping due to converged solution update"
+                        << std::endl;
+                    }
+                }
+                
+                if (iter == max_it) {
+                    
+                    terminate = true;
+                    std::cout
+                    << "Terminating Return-Mapping due to maximum iteration"
+                    << std::endl;
+                }
+                else if (res_norm < tol) {
+                    
+                    terminate = true;
+                    std::cout
+                    << "Terminating Return-Mapping due to converged residual"
+                    << std::endl;
+                }
+                else if (res_scaled_norm < tol) {
+                    
+                    terminate = true;
+                    std::cout
+                    << "Terminating Return-Mapping due to converged scaled residual"
+                    << std::endl;
+                }
+            }
+        }
+        
+        // set the plastic strain value based on the increment
+        Eigen::Matrix<scalar_t, n_strain, 1>
+        n    = Eigen::Matrix<scalar_t, n_strain, 1>::Zero();
+        
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<scalar_t, dim>::
+        derivative(internal.stress(), n);
+        
+        internal.plastic_strain()  = c.previous_plasticity_accessor->plastic_strain();
+        internal.plastic_strain() += (internal.consistency_parameter() -
+                                      c.previous_plasticity_accessor->consistency_parameter()) * n;
+
+        
+        // compute the stiffness tangent if the matrix was provided
+        if (stiff) {
+            
+            if (!_active_yield_function(yield))
+                // for elastic material response the standard material stiffness is acceptable
+                *stiff = m_stiff;
+            else
+                // for plastic material the tangent stiffness is to be computed using the current
+                // solution data
+                this->compute_tangent_stiffness(c, internal, *stiff);
+        }
+    }
+    
+    
+    /*!
+     * computes \f[ f(\sigma) = \sqrt(\sigma_{vm}) - \bar{\sigma}  \f], where \f$ \sigma_{vm} \f$ is the
+     * von-Mises stress and \f$ \bar{sigma} \f$ is the limit value.
+     */
+    template <typename ContextType,
+              typename AccessorType>
+    inline scalar_t yield_function(ContextType   &c,
+                                   AccessorType  &internal) const {
+
+        /*Eigen::Matrix<scalar_t, n_strain, 1>
+        s;
+        
+        for (uint_t i=0; i<n_strain; i++)
+            s(i) = stress(i);// - internal(i);
+         */
+        
+        scalar_t
+        v =  MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<scalar_t, dim>::
+        value(internal.stress());
+        
+        v -= _vm_lim;
+        
+        return v;
+    }
+    
+    
+
+    /*!
+     * computes \f$  C^{tan} = C - \frac{ C \partial_\sigma f C n }{ \partial_\sigma f C n } \f$
+     */
+    template <typename ContextType,
+              typename AccessorType>
+    inline void compute_tangent_stiffness(ContextType         &c,
+                                          const AccessorType  &internal,
+                                          stiff_t             &stiff) const {
+
+        Assert2(internal.stress().size() == n_strain,
+                internal.stress().size(), n_strain,
+                "Incompatible vector size");
+        Assert2(stiff.rows() == n_strain,
+                stiff.rows(), n_strain,
+                "Incompatible matrix row dimension");
+        Assert2(stiff.cols() == n_strain,
+                stiff.cols(), n_strain,
+                "Incompatible matrix column dimension");
+
+
+        _material->value(c, stiff);
+
+        Eigen::Matrix<scalar_t, n_strain, 1>
+        n   = Eigen::Matrix<scalar_t, n_strain, 1>::Zero(n_strain),
+        Cn  = Eigen::Matrix<scalar_t, n_strain, 1>::Zero(n_strain);
+        
+        // computation of \partial f/partial \sigma
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<scalar_t, dim>::
+        derivative(internal.stress(), n);
+        
+        Cn  = stiff * n;
+        
+        stiff -= (Cn*Cn.transpose())/(n.dot(Cn));
+    }
+    
+    
+    template <typename ContextType,
+              typename AccessorType,
+              typename Vec1Type,
+              typename Vec2Type,
+              typename MatType>
+    inline void
+    return_mapping_residual_and_jacobian(ContextType    &c,
+                                         const Vec1Type &strain,
+                                         AccessorType   &accessor,
+                                         Vec2Type       &res,
+                                         MatType        *jac = nullptr) const {
+                
+        Eigen::Matrix<scalar_t, n_strain, 1>
+        n    = Eigen::Matrix<scalar_t, n_strain, 1>::Zero();
+        Eigen::Matrix<scalar_t, n_strain, n_strain>
+        dnds = Eigen::Matrix<scalar_t, n_strain, n_strain>::Zero();
+
+        stiff_t       m_stiff;
+
+        // material stiffness matrix
+        _material->value(c, m_stiff);
+
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<scalar_t, dim>::
+        derivative(accessor.stress(), n);
+        
+        // derivative of normal, needed for Jacobian
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<scalar_t, dim>::
+        second_derivative(accessor.stress(), dnds);
+
+        ////////////////////////////
+        // residual: stress update
+        ////////////////////////////
+        res.topRows(n_strain) =
+        accessor.stress() - c.previous_plasticity_accessor->stress() -
+        m_stiff * (strain - c.previous_plasticity_accessor->strain() -
+                   (accessor.consistency_parameter() -
+                    c.previous_plasticity_accessor->consistency_parameter()) * n);
+        
+        // yield criterion
+        res(n_strain) = yield_function(c, accessor);
+        
+        if (jac) {
+            
+            ////////////////////////////
+            // Jacobian
+            ////////////////////////////
+            jac->topLeftCorner(n_strain, n_strain) =
+            (accessor.consistency_parameter() -
+             c.previous_plasticity_accessor->consistency_parameter()) * m_stiff * dnds;
+            for (uint_t i=0; i<n_strain; i++) (*jac)(i,i) += 1.;
+            
+            jac->topRightCorner(n_strain, 1) = m_stiff * n;
+            
+            jac->bottomLeftCorner(1, n_strain) = n.transpose();
+        }
+    }
+    
+private:
+    
+    template <typename ValType>
+    inline bool _active_yield_function(const ValType &v) const {
+        
+        return v >= 0;
+    }
+
+    inline bool _active_yield_function(const complex_t &v) const {
+        
+        return v.real() >= 0;
+    }
+
+    
+    bool             _if_kinematic_hardening;
+    real_t           _vm_lim;
+    MaterialType    *_material;
+};
+
+} // namespace MAST
+} // namespace Physics
+} // namespace Elasticity
+} // namespace ElastoPlasticity
+
+#endif // __mast_elastoplasticity_von_mises_yield_criterion_h__
diff --git a/include/mast/solvers/eigen/nonlinear_solver.hpp b/include/mast/solvers/eigen/nonlinear_solver.hpp
new file mode 100644
index 0000000..a3360e9
--- /dev/null
+++ b/include/mast/solvers/eigen/nonlinear_solver.hpp
@@ -0,0 +1,208 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef __mast_numerics_nonlinear_solver_h__
+#define __mast_numerics_nonlinear_solver_h__
+
+
+// MAST includes
+#include <mast/base/parameter_data.hpp>
+
+
+namespace MAST {
+namespace Numerics {
+namespace EigenWrapper {
+
+template <typename VecType,
+          typename MatType,
+          typename ComputeKernelType,
+          typename LinearSolverType = Eigen::FullPivLU<MatType>,
+          typename LineSearchType = void>
+class NonlinearSolver {
+    
+public:
+    
+    static_assert(std::is_same<VecType::Scalar, MatType::Scalar>::value,
+                  "Vector and Matrix must use same scalar types");
+    using scalar_t        = typename VecType::Scalar;
+    using linear_solver_t = LinearSolverType;
+    
+    uint_t       output_indent_level;
+    uint_t       verbose_level;
+    uint_t       max_iterations;
+    uint_t       update_jacobian_frequency;
+    real_t       residual_atol;
+    real_t       residual_rtol;
+    real_t       dx_atol;
+    
+    
+    NonlinearSolver():
+    _linear_solver   (nullptr),
+    _line_search     (nullptr),
+    _compute         (nullptr),
+    _output_indent_level       (0),
+    _verbose_level             (1),
+    _max_iterations            (20),
+    _update_jacobian_frequency (1),
+    residual_atol              (1.e-10),
+    residual_rtol              (1.e-8),
+    dx_atol                    (1.e-10)
+    { }
+    
+    virtual ~NonlinearSolver() {
+        
+        delete _parameter_data;
+    }
+    
+    
+    inline void set_linear_solver(LinearSolverType &solver) {
+        
+        _linear_solver = &solver;
+    }
+
+    inline void set_line_search(LineSearchType &ls) {
+        
+        _line_search = &ls;
+    }
+
+    inline void set_compute_kernel(ComputeKernelType &compute) {
+        
+        _compute = &compute;
+    }
+
+    
+    inline MAST::Base::ParameterData&
+    get_parameter_data() {
+        return _parameter_data;
+    }
+
+    /*!
+     *  solve starting with the initial guess
+     */
+    template <typename ContextType>
+    inline void solve(ContextType c, VecType& x) {
+        
+        Assert0(_linear_solver, "Linear solver not set");
+        Assert0(_compute, "Residual and Jacobian compute object not set");
+        
+        bool
+        terminate = false;
+        
+        uint_t
+        n_dofs          = c->n_dofs(),
+        it              = 0;
+        
+        real_t
+        dx_norm  = 0.,
+        res0     = 0.,  // norm at first iteration
+        res_norm = 0.;  // norm at current iteration
+
+        VecType
+        res = VecType::Zero(n_dofs),
+        dx  = VecType::Zero(n_dofs),
+        x1  = VecType::Zero(n_dofs);
+        
+
+        while (!terminate) {
+            
+            if (it%jac_update_freq == 0) {
+                
+                _compute->compute(c, x, res, &jac);
+                _linear_solver->compute(jac);
+            }
+            else
+                _compute->compute(c, x, res, nullptr);
+
+            // first check for residual convergence
+            res_norm = res.norm();
+            
+            if (verbose > 0)
+                std::cout
+                << std::setw(output_indent_level) << " "
+                << "Iter: " << std::setw(5) << it << std::endl
+                << std::setw(output_indent_level+2) << " "
+                << "|| res ||_2 = "
+                << std::setw(20) << res_norm << std::endl;
+            
+            // check for convergence of residual, otherwise solve for next step
+            if (_check_residual_convergence(it, res0, res_norm)) {
+                
+                terminate = true;
+                
+                if (verbose > 0)
+                    std::cout
+                    << std::setw(output_indent_level+2) << " "
+                    << "Terminating due to residual norm convergence." << std::endl;
+            }
+            else {
+                
+                // solve for the solution update, which is also the
+                // search direction for the line search
+                dx = _linear_solver->solve(res);
+                
+                if (_line_search) {
+                    
+                    dx *= -1.;
+                    _line_search->compute(c, x0, dx, x1);
+                    dx = x1-x0;
+                }
+                else
+                    // use the standard NR udpate if no linear search is specified
+                    x1 = x0 - dx;
+
+                dx_norm = dx.norm();
+
+                if (verbose > 0)
+                    std::cout
+                    << std::setw(output_indent_level+2) << " "
+                    << "|| dx ||_2 = "
+                    << std::setw(20) << dx_norm << std::endl;
+
+                // check for the convergence of solution update
+                if (dx_norm <= dx_atol) {
+                    
+                    terminate = true;
+                    
+                    if (verbose > 0)
+                        std::cout
+                        << std::setw(output_indent_level+2) << " "
+                        << "Terminating due to solution update convergence." << std::endl;
+                }
+            }
+            
+            it++;
+        }
+    }
+    
+    
+protected:
+    
+    LinearSolverType          *_linear_solver;
+
+    LineSearchType            *_line_search;
+    
+    ComputeKernelType         *_compute;
+};
+
+} // namespace EigenWrapper
+} // namespace Numerics
+} // namespace MAST
+
+
+#endif // __mast_numerics_nonlinear_solver_h__
diff --git a/tests/physics/elasticity/CMakeLists.txt b/tests/physics/elasticity/CMakeLists.txt
index 7a981f7..0b7711e 100644
--- a/tests/physics/elasticity/CMakeLists.txt
+++ b/tests/physics/elasticity/CMakeLists.txt
@@ -84,16 +84,28 @@ set_tests_properties(vonMisesStressComplexStep
         FIXTURES_SETUP     vonMisesStressComplexStep)
 
 
+
+
 #von Mises stress: adol-c sensitivity
 if (MAST_ENABLE_ADOLC EQUAL 1)
     target_sources(mast_catch_tests
                     PRIVATE
-                    ${CMAKE_CURRENT_LIST_DIR}/von_mises_stress_adolc_sensitivity.cpp)
+                    ${CMAKE_CURRENT_LIST_DIR}/von_mises_stress_adolc_sensitivity.cpp
+                    ${CMAKE_CURRENT_LIST_DIR}/von_mises_yield_criterion.cpp)
+                    
     add_test(NAME vonMisesStressAdolC
              COMMAND $<TARGET_FILE:mast_catch_tests> -w NoTests "von_mises_stress_adolc")
     set_tests_properties(vonMisesStressAdolC
             PROPERTIES
             LABELS "SEQ"
             FIXTURES_SETUP     vonMisesStressAdolC)
+
+    #von Mises yield criterion
+    add_test(NAME vonMisesYieldCriterion
+             COMMAND $<TARGET_FILE:mast_catch_tests> -w NoTests "von_mises_yield_criterion")
+    set_tests_properties(vonMisesYieldCriterion
+            PROPERTIES
+            LABELS "SEQ"
+            FIXTURES_SETUP     vonMisesYieldCriterion)
 endif()
 
diff --git a/tests/physics/elasticity/flat_shell_face_pressure_load_complex_step_sensitivity.cpp b/tests/physics/elasticity/flat_shell_face_pressure_load_complex_step_sensitivity.cpp
index c48d994..f604fb8 100644
--- a/tests/physics/elasticity/flat_shell_face_pressure_load_complex_step_sensitivity.cpp
+++ b/tests/physics/elasticity/flat_shell_face_pressure_load_complex_step_sensitivity.cpp
@@ -69,7 +69,7 @@ struct Traits {
     using fe_data_t    = typename MAST::FEBasis::libMeshWrapper::FEData<Dim, fe_basis_t, fe_shape_t>;
     using fe_var_t     = typename MAST::FEBasis::FEVarData<BasisScalarType, NodalScalarType, SolScalarType, 1, Dim, Context, fe_shape_t>;
     using pressure_t   = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using press_load_t = typename MAST::Physics::Elasticity::ShellFacePressureLoad<fe_var_t, pressure_t, Context>;
+    using press_load_t = typename MAST::Physics::Elasticity::ShellFacePressureLoad<fe_var_t, pressure_t>;
 };
 
 
diff --git a/tests/physics/elasticity/isotropic_stiffness.cpp b/tests/physics/elasticity/isotropic_stiffness.cpp
index b99ee0a..6c12920 100644
--- a/tests/physics/elasticity/isotropic_stiffness.cpp
+++ b/tests/physics/elasticity/isotropic_stiffness.cpp
@@ -43,7 +43,7 @@ struct Traits {
     using scalar_t          = ScalarType;
     using modulus_t         = typename MAST::Base::ScalarConstant<ScalarType>;
     using nu_t              = typename MAST::Base::ScalarConstant<ScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<ScalarType, Dim, modulus_t, nu_t, Context>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<ScalarType, Dim, modulus_t, nu_t>;
 };
 
 
diff --git a/tests/physics/elasticity/linear_strain_energy.cpp b/tests/physics/elasticity/linear_strain_energy.cpp
index a584a53..7e3e974 100644
--- a/tests/physics/elasticity/linear_strain_energy.cpp
+++ b/tests/physics/elasticity/linear_strain_energy.cpp
@@ -68,8 +68,8 @@ struct Traits {
     using fe_var_t          = typename MAST::FEBasis::FEVarData<BasisScalarType, NodalScalarType, SolScalarType, Dim, Dim, Context, fe_shape_t>;
     using modulus_t         = typename MAST::Base::ScalarConstant<SolScalarType>;
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t, Context>;
-    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, Dim, Context>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::LinearContinuum::StrainEnergy<fe_var_t, prop_t, Dim>;
 };
 
 
diff --git a/tests/physics/elasticity/mindlin_plate_strain_energy_complex_step_sensitivity.cpp b/tests/physics/elasticity/mindlin_plate_strain_energy_complex_step_sensitivity.cpp
index 61108da..bac431c 100644
--- a/tests/physics/elasticity/mindlin_plate_strain_energy_complex_step_sensitivity.cpp
+++ b/tests/physics/elasticity/mindlin_plate_strain_energy_complex_step_sensitivity.cpp
@@ -71,9 +71,9 @@ struct Traits {
     using modulus_t         = typename MAST::Base::ScalarConstant<scalar_t>;
     using nu_t              = typename MAST::Base::ScalarConstant<scalar_t>;
     using thickness_t       = typename MAST::Base::ScalarConstant<scalar_t>;
-    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<scalar_t, 2, modulus_t, nu_t, Context>;
-    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<scalar_t, material_t, thickness_t, Context>;
-    using energy_t          = typename MAST::Physics::Elasticity::MindlinPlate::StrainEnergy<fe_var_t, section_t, Context>;
+    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<scalar_t, 2, modulus_t, nu_t>;
+    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<scalar_t, material_t, thickness_t>;
+    using energy_t          = typename MAST::Physics::Elasticity::MindlinPlate::StrainEnergy<fe_var_t, section_t>;
 };
 
 
diff --git a/tests/physics/elasticity/plate_bending_section_property_complex_step_sensitivity.cpp b/tests/physics/elasticity/plate_bending_section_property_complex_step_sensitivity.cpp
index 908e5c4..41e13cd 100644
--- a/tests/physics/elasticity/plate_bending_section_property_complex_step_sensitivity.cpp
+++ b/tests/physics/elasticity/plate_bending_section_property_complex_step_sensitivity.cpp
@@ -44,8 +44,8 @@ struct Traits {
     using modulus_t         = typename MAST::Base::ScalarConstant<ScalarType>;
     using nu_t              = typename MAST::Base::ScalarConstant<ScalarType>;
     using thickness_t       = typename MAST::Base::ScalarConstant<ScalarType>;
-    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<ScalarType, 2, modulus_t, nu_t, Context>;
-    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<ScalarType, material_t, thickness_t, Context>;
+    using material_t        = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<ScalarType, 2, modulus_t, nu_t>;
+    using section_t         = typename MAST::Physics::Elasticity::PlateBendingSectionProperty<ScalarType, material_t, thickness_t>;
 
 };
 
diff --git a/tests/physics/elasticity/pressure_load.cpp b/tests/physics/elasticity/pressure_load.cpp
index 8956432..b1591ec 100644
--- a/tests/physics/elasticity/pressure_load.cpp
+++ b/tests/physics/elasticity/pressure_load.cpp
@@ -71,7 +71,7 @@ struct Traits {
     using fe_var_t     = typename MAST::FEBasis::FEVarData<BasisScalarType, NodalScalarType, SolScalarType, Dim, Dim, Context, fe_shape_t>;
     using pressure_t   = typename MAST::Base::ScalarConstant<SolScalarType>;
     using area_t       = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using press_load_t = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, pressure_t, area_t, Dim, Context>;
+    using press_load_t = typename MAST::Physics::Elasticity::SurfacePressureLoad<fe_var_t, pressure_t, area_t, Dim>;
 };
 
 
diff --git a/tests/physics/elasticity/stress_evaluation.cpp b/tests/physics/elasticity/stress_evaluation.cpp
index 51417f6..1a97472 100644
--- a/tests/physics/elasticity/stress_evaluation.cpp
+++ b/tests/physics/elasticity/stress_evaluation.cpp
@@ -68,8 +68,8 @@ struct Traits {
     using fe_var_t          = typename MAST::FEBasis::FEVarData<BasisScalarType, NodalScalarType, SolScalarType, Dim, Dim, Context, fe_shape_t>;
     using modulus_t         = typename MAST::Base::ScalarConstant<SolScalarType>;
     using nu_t              = typename MAST::Base::ScalarConstant<SolScalarType>;
-    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t, Context>;
-    using stress_t          = typename MAST::Physics::Elasticity::LinearContinuum::Stress<fe_var_t, prop_t, Dim, Context>;
+    using prop_t            = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<SolScalarType, Dim, modulus_t, nu_t>;
+    using stress_t          = typename MAST::Physics::Elasticity::LinearContinuum::Stress<fe_var_t, prop_t, Dim>;
     using stress_vec_t      = typename Eigen::Matrix<scalar_t, stress_t::n_strain, 1>;
     using stress_adj_mat_t  = typename Eigen::Matrix<scalar_t, stress_t::n_strain, Eigen::Dynamic>;
     using stress_storage_t  = typename Eigen::Matrix<scalar_t, stress_t::n_strain, Eigen::Dynamic>;
diff --git a/tests/physics/elasticity/von_mises_stress_adolc_sensitivity.cpp b/tests/physics/elasticity/von_mises_stress_adolc_sensitivity.cpp
index 80d5e70..e654561 100644
--- a/tests/physics/elasticity/von_mises_stress_adolc_sensitivity.cpp
+++ b/tests/physics/elasticity/von_mises_stress_adolc_sensitivity.cpp
@@ -33,13 +33,6 @@ namespace Elasticity {
 namespace vonMisesStress {
 namespace AdolC {
 
-template <typename ScalarType, uint_t Dim>
-using stress_vec_t = MAST::Physics::Elasticity::LinearContinuum::stress_vec_t<ScalarType, Dim>;
-
-template <typename ScalarType, uint_t Dim>
-using stress_adjoint_mat_t = MAST::Physics::Elasticity::LinearContinuum::stress_adjoint_mat_t<ScalarType, Dim>;
-
-
 
 template <uint_t   Dim>
 inline void test_von_mises_stress_sensitivity()  {
@@ -49,30 +42,48 @@ inline void test_von_mises_stress_sensitivity()  {
     n_strain = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value;
     
     
-    stress_vec_t<real_t, Dim>
-    stress   = stress_vec_t<real_t, Dim>::Random(),
-    dstress  = stress_vec_t<real_t, Dim>::Random();
+    Eigen::Matrix<real_t, n_strain, 1>
+    stress   = Eigen::Matrix<real_t, n_strain, 1>::Random(),
+    dstress  = Eigen::Matrix<real_t, n_strain, 1>::Random(),
+    dsvm     = Eigen::Matrix<real_t, n_strain, 1>::Zero(),
+    dsvm_ad  = Eigen::Matrix<real_t, n_strain, 1>::Zero();
     
-    stress_adjoint_mat_t<real_t, Dim>
-    stress_adj_mat = stress_adjoint_mat_t<real_t, Dim>::Random(n_strain, n_basis);
+
+    Eigen::Matrix<real_t, n_strain, n_strain>
+    d2stress       = Eigen::Matrix<real_t, n_strain, n_strain>::Zero(),
+    d2stress_ad    = Eigen::Matrix<real_t, n_strain, n_strain>::Zero();
+
+    Eigen::Matrix<real_t, n_strain, Eigen::Dynamic>
+    stress_adj_mat = Eigen::Matrix<real_t, n_strain, Eigen::Dynamic>::Random(n_strain, n_basis);
 
     Eigen::Matrix<real_t, Eigen::Dynamic, 1>
     vm_adj    = Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Zero(n_basis),
-    vm_adj_cs = Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Zero(n_basis);
+    vm_adj_ad = Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Zero(n_basis);
     
     real_t
-    vm  = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<real_t, Dim>(stress),
-    dvm = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_derivative<real_t, Dim>(stress, dstress),
-    dvm_cs = 0.;
+    vm  = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    value(stress),
+    dvm = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    derivative_sens(stress, dstress),
+    dvm_ad = 0.;
+
+    // stress adjoint
+    MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    stress_dX(stress, stress_adj_mat, vm_adj);
     
-    MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_dX<real_t, Dim>(stress,
-                                                                                stress_adj_mat,
-                                                                                vm_adj);
+    // dsigma_vm/dsigma
+    MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    derivative(stress, dsvm);
+
+    // d2sigma_vm/dsigma2
+    MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    second_derivative(stress, d2stress);
     
     {
         // the number of directions for which we compute the sensitivity is = n_strain
-        stress_vec_t<adouble_tl_t, Dim>
-        stress_ad;
+        Eigen::Matrix<adouble_tl_t, n_strain, 1>
+        stress_ad,
+        ds_ad;
         
         for (uint_t i=0; i<n_strain; i++) {
             
@@ -81,16 +92,48 @@ inline void test_von_mises_stress_sensitivity()  {
         }
         
         adouble_tl_t
-        vm_ad = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<adouble_tl_t, Dim>(stress_ad);
+        vm_ad = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<adouble_tl_t, Dim>::
+        value(stress_ad);
 
-        dvm_cs = *vm_ad.getADValue();
+        dvm_ad = *vm_ad.getADValue();
     }
+
+    
+    // first and second derivative of sigma_vm wrt stress vector
+    for (uint_t i=0; i<n_strain; i++) {
+        
+        // the number of directions for which we compute the sensitivity is = n_strain
+        Eigen::Matrix<adouble_tl_t, n_strain, 1>
+        stress_ad,
+        ds_ad;
+        
+        real_t
+        v = 1.;
+        
+        for (uint_t j=0; j<n_strain; j++) stress_ad(j) = stress(j);
+        
+        // derivative wrt ith sigma value
+        stress_ad(i).setADValue(&v);
+
+        adouble_tl_t
+        vm_ad = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<adouble_tl_t, Dim>::
+        value(stress_ad);
+
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<adouble_tl_t, Dim>::
+        derivative(stress_ad, ds_ad);
+        
+        dsvm_ad(i) = *vm_ad.getADValue();
+
+        for (uint_t j=0; j<n_strain; j++)
+            d2stress_ad(j,i) = *ds_ad(j).getADValue();
+    }
+    
     
     // adjoint
     {
         adtl::setNumDir(n_basis);
 
-        stress_vec_t<adouble_tl_t, Dim>
+        Eigen::Matrix<adouble_tl_t, n_strain, 1>
         stress_ad;
 
         // the adjoint can be computed in adol-c traceless vector mode
@@ -106,16 +149,23 @@ inline void test_von_mises_stress_sensitivity()  {
         }
         
         adouble_tl_t
-        vm_ad = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<adouble_tl_t, Dim>(stress_ad);
+        vm_ad = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<adouble_tl_t, Dim>::
+        value(stress_ad);
         
         for (uint_t j=0; j<n_basis; j++)
-            vm_adj_cs(j) = vm_ad.getADValue(j);
+            vm_adj_ad(j) = vm_ad.getADValue(j);
     }
     
-    CHECK(dvm == Catch::Detail::Approx(dvm_cs));
+    CHECK(dvm == Catch::Detail::Approx(dvm_ad));
 
     CHECK_THAT(MAST::Test::eigen_matrix_to_std_vector(vm_adj),
-               Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(vm_adj_cs)));
+               Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(vm_adj_ad)));
+
+    CHECK_THAT(MAST::Test::eigen_matrix_to_std_vector(dsvm),
+               Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(dsvm_ad)));
+
+    CHECK_THAT(MAST::Test::eigen_matrix_to_std_vector(d2stress),
+               Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(d2stress_ad)));
 }
 
 
diff --git a/tests/physics/elasticity/von_mises_stress_complex_step_sensitivity.cpp b/tests/physics/elasticity/von_mises_stress_complex_step_sensitivity.cpp
index 094e797..f14f719 100644
--- a/tests/physics/elasticity/von_mises_stress_complex_step_sensitivity.cpp
+++ b/tests/physics/elasticity/von_mises_stress_complex_step_sensitivity.cpp
@@ -33,13 +33,6 @@ namespace Elasticity {
 namespace vonMisesStress {
 namespace ComplexStep {
 
-template <typename ScalarType, uint_t Dim>
-using stress_vec_t = MAST::Physics::Elasticity::LinearContinuum::stress_vec_t<ScalarType, Dim>;
-
-template <typename ScalarType, uint_t Dim>
-using stress_adjoint_mat_t = MAST::Physics::Elasticity::LinearContinuum::stress_adjoint_mat_t<ScalarType, Dim>;
-
-
 
 template <uint_t   Dim>
 inline void test_von_mises_stress_sensitivity()  {
@@ -49,42 +42,43 @@ inline void test_von_mises_stress_sensitivity()  {
     n_strain = MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value;
     
     
-    stress_vec_t<real_t, Dim>
-    stress   = stress_vec_t<real_t, Dim>::Random(),
-    dstress  = stress_vec_t<real_t, Dim>::Random();
+    Eigen::Matrix<real_t, n_strain, 1>
+    stress   = Eigen::Matrix<real_t, n_strain, 1>::Random(),
+    dstress  = Eigen::Matrix<real_t, n_strain, 1>::Random();
     
-    stress_adjoint_mat_t<real_t, Dim>
-    stress_adj_mat = stress_adjoint_mat_t<real_t, Dim>::Random(n_strain, n_basis);
+    Eigen::Matrix<real_t, n_strain, Eigen::Dynamic>
+    stress_adj_mat = Eigen::Matrix<real_t, n_strain, Eigen::Dynamic>::Random(n_strain, n_basis);
 
     Eigen::Matrix<real_t, Eigen::Dynamic, 1>
     vm_adj    = Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Zero(n_basis),
     vm_adj_cs = Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Zero(n_basis);
     
     real_t
-    vm  = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<real_t, Dim>(stress),
-    dvm = MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_derivative<real_t, Dim>(stress, dstress),
+    vm  = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    value(stress),
+    dvm = MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    derivative_sens(stress, dstress),
     dvm_cs = 0.;
     
-    MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_dX<real_t, Dim>(stress,
-                                                                                stress_adj_mat,
-                                                                                vm_adj);
+    MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<real_t, Dim>::
+    stress_dX(stress, stress_adj_mat, vm_adj);
     
     // extremely small perturbations do not work for this routine that involves
     // square and square-root operations
     for (uint_t i=0; i<n_strain; i++) {
         
-        stress_vec_t<complex_t, Dim>
+        Eigen::Matrix<complex_t, n_strain, 1>
         stress_c   = stress.template cast<complex_t>();
         
         stress_c(i) += complex_t(0., sqrt(ComplexStepDelta));
         
         // linearized contribution of ith stress component
         dvm_cs += dstress(i)/sqrt(ComplexStepDelta) *
-        MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<complex_t, Dim>(stress_c).imag();
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<complex_t, Dim>::value(stress_c).imag();
 
         // linearized contribution of ith stress component
         vm_adj_cs += stress_adj_mat.row(i)/sqrt(ComplexStepDelta) *
-        MAST::Physics::Elasticity::LinearContinuum::vonMises_stress<complex_t, Dim>(stress_c).imag();
+        MAST::Physics::Elasticity::LinearContinuum::vonMisesStress<complex_t, Dim>::value(stress_c).imag();
     }
     
     CHECK(dvm == Catch::Detail::Approx(dvm_cs));
diff --git a/tests/physics/elasticity/von_mises_yield_criterion.cpp b/tests/physics/elasticity/von_mises_yield_criterion.cpp
new file mode 100644
index 0000000..5a4ed12
--- /dev/null
+++ b/tests/physics/elasticity/von_mises_yield_criterion.cpp
@@ -0,0 +1,324 @@
+/*
+* MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
+* Copyright (C) 2013-2020  Manav Bhatia and MAST authors
+*
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* This library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this library; if not, write to the Free Software
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+// Catch includes
+#include "catch.hpp"
+
+// MAST includes
+#include <mast/physics/elasticity/von_mises_yield_criterion.hpp>
+#include <mast/physics/elasticity/isotropic_stiffness.hpp>
+#include <mast/base/scalar_constant.hpp>
+
+// Test includes
+#include <test_helpers.h>
+
+namespace MAST {
+namespace Test {
+namespace Physics {
+namespace Elasticity {
+namespace vonMisesYieldCriterion {
+
+template <typename AccessorType>
+struct Context {
+    
+    Context():
+    previous_plasticity_accessor (nullptr)
+    { }
+    
+    
+    AccessorType *previous_plasticity_accessor;
+};
+
+
+template <typename ScalarType, uint_t Dim>
+struct Traits {
+
+    using scalar_t  = ScalarType;
+    using modulus_t = typename MAST::Base::ScalarConstant<ScalarType>;
+    using nu_t      = typename MAST::Base::ScalarConstant<ScalarType>;
+    using prop_t    = typename MAST::Physics::Elasticity::IsotropicMaterialStiffness<ScalarType, Dim, modulus_t, nu_t>;
+    using yield_t   = MAST::Physics::Elasticity::ElastoPlasticity::vonMisesYieldFunction<ScalarType, prop_t>;
+    using accessor_t= typename MAST::Physics::Elasticity::ElastoPlasticity::Accessor<scalar_t, yield_t>;
+    static const uint_t
+    n_strain        = yield_t::n_strain;
+};
+
+
+template <typename Traits>
+struct Prop {
+    
+    Prop():
+    E    (new typename Traits::modulus_t(72.e9)),
+    nu   (new typename Traits::nu_t(0.33)),
+    prop (new typename Traits::prop_t) {
+        
+        prop->set_modulus_and_nu(*E, *nu);
+    }
+    
+    virtual ~Prop() { }
+    
+    std::unique_ptr<typename Traits::modulus_t>  E;
+    std::unique_ptr<typename Traits::nu_t>       nu;
+    std::unique_ptr<typename Traits::prop_t>     prop;
+};
+
+
+
+inline void copy_value(real_t                                    factor,
+                       Eigen::Matrix<real_t, Eigen::Dynamic, 1>& from,
+                       Eigen::Matrix<real_t, Eigen::Dynamic, 1>& to) {
+    
+    to = factor * from;
+}
+
+
+// this removes the AD-derivative data from copy
+inline void copy_value(real_t                                          factor,
+                       Eigen::Matrix<adouble_tl_t, Eigen::Dynamic, 1>& from,
+                       Eigen::Matrix<adouble_tl_t, Eigen::Dynamic, 1>& to) {
+    
+    for (uint_t i=0; i<from.size(); i++) to(i) = factor * from(i).getValue();
+}
+
+
+
+template <typename Traits>
+inline void von_mises_yield_criterion_jacobian
+(Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain,   1> &strain,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain+1, 1> &x,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain+1, 1> &res,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain+1, Traits::n_strain+1> *jac) {
+    
+    using scalar_t         = typename Traits::scalar_t;
+    
+    Context<typename Traits::accessor_t> c;
+    Prop<Traits> p;
+    
+    typename Traits::yield_t
+    yield;
+
+    const uint_t
+    n_strain = Traits::n_strain,
+    n_dofs   = yield.n_variables();
+    
+
+    Eigen::Matrix<scalar_t, n_strain, n_strain>
+    stiff = Eigen::Matrix<scalar_t, n_strain, n_strain>::Zero();
+
+    Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>
+    v0 = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Zero(n_dofs),
+    v1 = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Zero(n_dofs);
+
+    v1.topRows(n_strain) = x.topRows(n_strain); // copy stress
+    v1(n_dofs-1)         = x(n_strain);         // copy consistency parameter
+    
+    copy_value(0.8, v1, v0);
+    
+    MAST::Physics::Elasticity::ElastoPlasticity::Accessor<scalar_t, typename Traits::yield_t>
+    accessor0,
+    accessor1;
+    
+    accessor0.init(yield, v0.data());
+    accessor1.init(yield, v1.data());
+        
+    c.previous_plasticity_accessor = &accessor0;
+    
+    yield.set_material(*p.prop);
+    yield.set_limit_stress(5.e6);
+    yield.return_mapping_residual_and_jacobian(c, strain, accessor1, res, jac);
+}
+
+
+
+template <typename Traits>
+inline void von_mises_yield_criterion_tangent_stiffness
+(Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain,   1> &strain,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain+1, 1> &x,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain, 1> &stress,
+ Eigen::Matrix<typename Traits::scalar_t, Traits::n_strain, Traits::n_strain> *Cmat) {
+    
+    using scalar_t         = typename Traits::scalar_t;
+    
+    Context<typename Traits::accessor_t> c;
+    Prop<Traits> p;
+    
+    typename Traits::yield_t
+    yield;
+
+    const uint_t
+    n_strain = Traits::n_strain,
+    n_dofs   = yield.n_variables();
+    
+
+    Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>
+    v0 = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Zero(n_dofs),
+    v1 = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Zero(n_dofs);
+
+    v1.topRows(n_strain) = x.topRows(n_strain);  // copy stress
+    v1(n_dofs-1)         = x(n_strain);          // copy consistency parameter
+    
+    copy_value(0.8, v1, v0);
+    
+    MAST::Physics::Elasticity::ElastoPlasticity::Accessor<scalar_t, typename Traits::yield_t>
+    accessor0,
+    accessor1;
+    
+    accessor0.init(yield, v0.data());
+    accessor1.init(yield, v1.data());
+        
+    c.previous_plasticity_accessor = &accessor0;
+
+    yield.set_material(*p.prop);
+    yield.set_limit_stress(5.e6);
+    yield.compute(c, strain, accessor1, Cmat);
+
+    /*// we use this as the starting point for the
+    for (uint_t i=0; i<v0.size(); i++)
+    std::cout << std::setw(20) << v0(i) << std::setw(20) << v1(i) << std::endl;
+    copy_value(1., v1, v0);
+    yield.compute(c, strain, accessor1, Cmat);
+
+    for (uint_t i=0; i<v0.size(); i++)
+    std::cout << std::setw(20) << v0(i) << std::setw(20) << v1(i) << std::endl;*/
+    // copy the stress value back
+    stress = accessor1.stress();
+}
+
+
+
+inline void
+test_von_mises_yield_criterion_jacobian() {
+    
+    using traits_t         = Traits<real_t, 2>;
+    
+    const uint_t
+    n_strain = traits_t::yield_t::n_strain;
+
+    Eigen::Matrix<real_t, n_strain, 1>
+    strain = 1.e-4 * Eigen::Matrix<real_t, n_strain, 1>::Random(),
+    stress;
+
+    Eigen::Matrix<real_t, n_strain, n_strain>
+    Cmat    = Eigen::Matrix<real_t, n_strain, n_strain>::Zero(),
+    Cmat_ad = Eigen::Matrix<real_t, n_strain, n_strain>::Zero();
+
+    // the return-mapping solver will stack variables in x as
+    // {stress, consistency_param}
+    Eigen::Matrix<real_t, n_strain+1, 1>
+    x      = Eigen::Matrix<real_t, n_strain+1, 1>::Random(),
+    res    = Eigen::Matrix<real_t, n_strain+1, 1>::Zero();
+    // scale the stresses by factor
+    x.topRows(n_strain) *= 1.e2;
+    
+    Eigen::Matrix<real_t, n_strain+1, n_strain+1>
+    jac    = Eigen::Matrix<real_t, n_strain+1, n_strain+1>::Zero(),
+    jac_ad = Eigen::Matrix<real_t, n_strain+1, n_strain+1>::Zero();
+    
+    von_mises_yield_criterion_jacobian<traits_t>(strain, x, res, &jac);
+    von_mises_yield_criterion_tangent_stiffness<traits_t>(strain, x, stress, &Cmat);
+    
+    // compute the quantities using adouble_tl_t to verify the linearization
+    {
+        adtl::setNumDir(n_strain+1);
+        using traits_ad_t = Traits<adouble_tl_t, 2>;
+        
+        Eigen::Matrix<adouble_tl_t, n_strain, 1>
+        strain_ad = Eigen::Matrix<adouble_tl_t, n_strain, 1>::Zero();
+                
+        Eigen::Matrix<adouble_tl_t, n_strain+1, 1>
+        x_ad   = Eigen::Matrix<adouble_tl_t, n_strain+1, 1>::Zero(),
+        res_   = Eigen::Matrix<adouble_tl_t, n_strain+1, 1>::Zero();
+        
+        Eigen::Matrix<adouble_tl_t, n_strain+1, n_strain+1>
+        *jac_ = nullptr;
+        
+        // the adjoint can be computed in adol-c traceless vector mode
+        // with ndof components.
+        for (uint_t i=0; i<n_strain; i++) strain_ad(i) = strain(i);
+        
+        for (uint_t i=0; i<n_strain+1; i++) {
+            
+            x_ad(i) = x(i);
+            x_ad(i).setADValue(i, 1.);
+        }
+        
+        von_mises_yield_criterion_jacobian<traits_ad_t>(strain_ad, x_ad, res_, jac_);
+
+        for (uint_t i=0; i<n_strain+1; i++)
+            for (uint_t j=0; j<n_strain+1; j++)
+                jac_ad(i,j) = res_(i).getADValue(j);
+    }
+
+    
+    // compute the quantities using adouble_tl_t to verify the tangent stiffness
+    {
+        adtl::setNumDir(n_strain);
+        using traits_ad_t = Traits<adouble_tl_t, 2>;
+        
+        Eigen::Matrix<adouble_tl_t, n_strain, 1>
+        strain_ad = Eigen::Matrix<adouble_tl_t, n_strain, 1>::Zero(),
+        stress_;
+                
+        Eigen::Matrix<adouble_tl_t, n_strain, n_strain>
+        *Cmat_ = nullptr;
+
+        Eigen::Matrix<adouble_tl_t, n_strain+1, 1>
+        x_ad   = Eigen::Matrix<adouble_tl_t, n_strain+1, 1>::Zero();
+        
+        Eigen::Matrix<adouble_tl_t, n_strain+1, n_strain+1>
+        *jac_ = nullptr;
+        
+        // the adjoint can be computed in adol-c traceless vector mode
+        // with ndof components.
+        for (uint_t i=0; i<n_strain; i++) {
+          
+            strain_ad(i) = strain(i);
+            strain_ad(i).setADValue(i, 1.);
+        }
+        
+        //for (uint_t i=0; i<n_strain+1; i++) x_ad(i) = x(i);
+        
+        von_mises_yield_criterion_tangent_stiffness<traits_ad_t>(strain_ad, x_ad, stress_, Cmat_);
+
+        for (uint_t i=0; i<n_strain; i++)
+            for (uint_t j=0; j<n_strain; j++)
+                Cmat_ad(i,j) = stress_(i).getADValue(j);
+    }
+
+    
+    CHECK_THAT(MAST::Test::eigen_matrix_to_std_vector(jac),
+               Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(jac_ad)));
+
+    //std::cout << Cmat << std::endl << Cmat_ad << std::endl;
+    //CHECK_THAT(MAST::Test::eigen_matrix_to_std_vector(Cmat),
+    //           Catch::Approx(MAST::Test::eigen_matrix_to_std_vector(Cmat_ad)));
+}
+
+
+TEST_CASE("von_mises_yield_criterion",
+          "[Physics][Elasticity][Elastoplasticity][vonMisesPlasticity]") {
+    
+    vonMisesYieldCriterion::test_von_mises_yield_criterion_jacobian();
+}
+
+} // namespace MAST
+} // namespace Test
+} // namespace Physics
+} // namespace Elasticity
+} // namespace vonMisesYieldCriterion