Index: include/systems/fem_context.h
===================================================================
--- include/systems/fem_context.h	(revision 5024)
+++ include/systems/fem_context.h	(working copy)
@@ -78,44 +78,44 @@
    * Returns the value of the solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Number interior_value(unsigned int var, unsigned int qp);
+  Number interior_value(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the value of the solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Number side_value(unsigned int var, unsigned int qp);
+  Number side_value(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the value of the solution variable \p var at the physical
    * point \p p on the current element
    */
-  Number point_value(unsigned int var, const Point &p);
+  Number point_value(unsigned int var, const Point &p) const;
 
   /**
    * Returns the gradient of the solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Gradient interior_gradient(unsigned int var, unsigned int qp);
+  Gradient interior_gradient(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the gradient of the solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Gradient side_gradient(unsigned int var, unsigned int qp);
+  Gradient side_gradient(unsigned int var, unsigned int qp) const;
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
   /**
    * Returns the hessian of the solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Tensor interior_hessian(unsigned int var, unsigned int qp);
+  Tensor interior_hessian(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the hessian of the solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Tensor side_hessian(unsigned int var, unsigned int qp);
+  Tensor side_hessian(unsigned int var, unsigned int qp) const;
 
 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
 
@@ -123,44 +123,44 @@
    * Returns the value of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Number fixed_interior_value(unsigned int var, unsigned int qp);
+  Number fixed_interior_value(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the value of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Number fixed_side_value(unsigned int var, unsigned int qp);
+  Number fixed_side_value(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the value of the fixed_solution variable \p var at the physical
    * point \p p on the current element
    */
-  Number fixed_point_value(unsigned int var, const Point &p);
+  Number fixed_point_value(unsigned int var, const Point &p) const;
 
   /**
    * Returns the gradient of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Gradient fixed_interior_gradient(unsigned int var, unsigned int qp);
+  Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the gradient of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Gradient fixed_side_gradient(unsigned int var, unsigned int qp);
+  Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const;
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
   /**
    * Returns the hessian of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element interior
    */
-  Tensor fixed_interior_hessian(unsigned int var, unsigned int qp);
+  Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const;
 
   /**
    * Returns the hessian of the fixed_solution variable \p var at the quadrature
    * point \p qp on the current element side
    */
-  Tensor fixed_side_hessian(unsigned int var, unsigned int qp);
+  Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const;
 
 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
 
Index: src/systems/fem_context.C
===================================================================
--- src/systems/fem_context.C	(revision 5024)
+++ src/systems/fem_context.C	(working copy)
@@ -140,7 +140,7 @@
 
 
 
-Number FEMContext::interior_value(unsigned int var, unsigned int qp)
+Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -166,7 +166,7 @@
 
 
 
-Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp)
+Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -193,7 +193,7 @@
 
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
-Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp)
+Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -220,7 +220,7 @@
 
 
 
-Number FEMContext::side_value(unsigned int var, unsigned int qp)
+Number FEMContext::side_value(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -246,7 +246,7 @@
 
 
 
-Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp)
+Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -273,7 +273,7 @@
 
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
-Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp)
+Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -300,7 +300,7 @@
 
 
 
-Number FEMContext::point_value(unsigned int var, const Point &p)
+Number FEMContext::point_value(unsigned int var, const Point &p) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -326,7 +326,7 @@
 
 
 
-Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp)
+Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -352,7 +352,7 @@
 
 
 
-Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp)
+Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -379,7 +379,7 @@
 
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
-Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp)
+Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -406,7 +406,7 @@
 
 
 
-Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp)
+Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -432,7 +432,7 @@
 
 
 
-Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp)
+Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -459,7 +459,7 @@
 
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
-Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp)
+Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);
@@ -486,7 +486,7 @@
 
 
 
-Number FEMContext::fixed_point_value(unsigned int var, const Point &p)
+Number FEMContext::fixed_point_value(unsigned int var, const Point &p) const
 {
   // Get local-to-global dof index lookup
   libmesh_assert (dof_indices.size() > var);