diff --git a/cardmech/src/elasticity/assemble/DGElasticity.cpp b/cardmech/src/elasticity/assemble/DGElasticity.cpp
index b49274155d3c7f08ee9832e42ff7c536b634cd72..526b1602066da85220603a7a639ae0fe42a8911e 100644
--- a/cardmech/src/elasticity/assemble/DGElasticity.cpp
+++ b/cardmech/src/elasticity/assemble/DGElasticity.cpp
@@ -729,16 +729,30 @@ Vector DGElasticity::Projection(const Vector &u, const Vector &reference) const
     }
     auto r = projection.find_row(c());
     int childCount = recursiveCells[dlevel].size();
-    DGVectorFieldElement E(u, c);
-    int nV = E.shape_size() * E.Dim();
-    for (const auto &child: recursiveCells[dlevel]) {
-      auto fineRow = reference.find_row(child->Center());
-      for(int i=0;i<nV;++i) {
-        projection(r, i) += reference(fineRow, i) / childCount;
+    DGVectorFieldElement E(projection, c);
+
+    for(int i = 0;i< E.NodalPoints(); ++i){
+      Point z = E.NodalPoint(i);
+      bool found=false;
+      for (const auto &child: recursiveCells[dlevel]) {
+        DGVectorFieldElement CE(reference, reference.find_cell(child->Center()));
+        for(int j = 0;j< CE.NodalPoints(); ++j){
+          if(CE.NodalPoint(j) == z){
+            auto rf = reference.find_row(child->Center());
+            for(int k=0;k<E.Dim();++k) {
+              projection(r, k*E.shape_size() + i) = reference(rf, k*CE.shape_size() + j);
+            }
+            found = true;
+          }
+          if(found)
+            break;
+        }
+        if(found)
+          break;
       }
     }
   }
-    return projection;
+  return projection;
 }
 
 
@@ -815,7 +829,7 @@ double DGElasticity::EnergyError(const Vector &u) const {
       double w = E.QWeight(q);
       const Point &z = E.QPoint(q);
       Tensor DD = (E.VectorGradient(q, u) - eProblem.DeformationGradient(timeStep, z));
-      err += w * abs(cellMat.IsoDerivative(sym(One + DD), DD));
+      err += w * cellMat.IsoDerivative(sym(One + DD), DD);
       err += w * abs(cellMat.AnisoDerivative(One + DD, Q, DD));
     }
   }
@@ -833,7 +847,7 @@ double DGElasticity::EnergyNorm(const Vector &u) const {
         for (int q = 0; q < E.nQ(); ++q) {
             double w = E.QWeight(q);
             Tensor DD = E.VectorGradient(q, u);
-            err += w * abs( cellMat.TotalDerivative(One + DD, One + DD, Q , DD));
+          err += w * abs(cellMat.TotalDerivative(One + DD, One+DD, Q, DD));
         }
     }
     return sqrt(PPM->Sum(err));
diff --git a/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp b/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
index 0ea525b560e917e14459f6c60deda3dba16eed90..47f205f54a01fc73b1493c938ef607d726cd8e71 100644
--- a/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
+++ b/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
@@ -50,7 +50,7 @@ public:
                                                                             isPlanar(planar) {}
 
   string MeshName() const override {
-    return (isPlanar ? "UnitSquare" : "UnitCube");
+    return (isPlanar ? "UnitBlock2DQuad" : "UnitBlock3DQuad");
   }
 
   VectorField Deformation(const Point &z) const {
diff --git a/tools/jobsystem/run_dg_benchmarks.py b/tools/jobsystem/run_dg_benchmarks.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391