diff --git a/scripts/lc-builds/toss3_gcc8.1.0-debug-mpi-cxx.sh b/scripts/lc-builds/toss3_gcc8.1.0-debug-mpi-cxx.sh
new file mode 100755
index 000000000..3f9400d73
--- /dev/null
+++ b/scripts/lc-builds/toss3_gcc8.1.0-debug-mpi-cxx.sh
@@ -0,0 +1,25 @@
+#!/usr/bin/env bash
+
+SCRIPT_PATH=${0%/*}
+. "$SCRIPT_PATH/utils/parse-args.sh"
+
+# Inherit build directory name from script name
+BUILD_SUFFIX="lc_$(TMP=${BASH_SOURCE##*/}; echo ${TMP%.*})"
+
+rm -rf ${BUILD_SUFFIX} 2>/dev/null
+mkdir -p ${BUILD_SUFFIX}/install
+mkdir -p ${BUILD_SUFFIX}/build && cd ${BUILD_SUFFIX}/build
+
+module load cmake/3.14.5
+module load gcc/8.1.0
+
+cmake \
+  ${SRC_DIR} \
+  -DCMAKE_BUILD_TYPE=Debug \
+  -DCMAKE_CXX_COMPILER=/usr/tce/packages/gcc/gcc-8.1.0/bin/g++ \
+  -C ${HOST_CONFIGS_DIR}/lc-builds/toss3/gcc8.1.0_tpl.cmake \
+  -DENABLE_OPENMP=On \
+  -DENABLE_MPI=On \
+  -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR \
+  -DENABLE_STATIC_CXXONLY=On \
+  $CMAKE_ARGS \
diff --git a/src/CRKSPH/SolidCRKSPHHydroBase.cc b/src/CRKSPH/SolidCRKSPHHydroBase.cc
index 705e8d27d..d3b56ec55 100644
--- a/src/CRKSPH/SolidCRKSPHHydroBase.cc
+++ b/src/CRKSPH/SolidCRKSPHHydroBase.cc
@@ -35,7 +35,6 @@
 #include "Neighbor/ConnectivityMap.hh"
 #include "Utilities/timingUtilities.hh"
 #include "Utilities/safeInv.hh"
-#include "Utilities/ThreePointDamagedNodeCoupling.hh"
 #include "SolidMaterial/SolidEquationOfState.hh"
 
 #include "SolidCRKSPHHydroBase.hh"
diff --git a/src/Damage/CMakeLists.txt b/src/Damage/CMakeLists.txt
index e7e96b8fc..f1b0bd67b 100644
--- a/src/Damage/CMakeLists.txt
+++ b/src/Damage/CMakeLists.txt
@@ -16,6 +16,9 @@ set(Damage_inst
     JohnsonCookFailureStrainPolicy
     JohnsonCookDamagePolicy
     DamagedPressurePolicy
+    PairMaxDamageNodeCoupling
+    ThreePointDamagedNodeCoupling
+    DamageGradientNodeCoupling
    )
 
 
@@ -51,6 +54,10 @@ set(Damage_headers
     oneMinusDamage.hh
     weibullFlawDistributionBenzAsphaug.hh
     weibullFlawDistributionOwen.hh
+    DamageGradientNodeCoupling.hh
+    PairMaxDamageNodeCoupling.hh
+    MaxDamageNodeCoupling.hh
+    ThreePointDamagedNodeCoupling.hh
     )
 
 spheral_install_python_files(
diff --git a/src/Utilities/DamageGradientNodeCoupling.cc b/src/Damage/DamageGradientNodeCoupling.cc
similarity index 98%
rename from src/Utilities/DamageGradientNodeCoupling.cc
rename to src/Damage/DamageGradientNodeCoupling.cc
index 4fc69f5b9..38ce07074 100644
--- a/src/Utilities/DamageGradientNodeCoupling.cc
+++ b/src/Damage/DamageGradientNodeCoupling.cc
@@ -10,7 +10,7 @@
 //
 // Created by JMO, Fri Jul 31 14:46:25 PDT 2015
 //----------------------------------------------------------------------------//
-#include "Utilities/DamageGradientNodeCoupling.hh"
+#include "Damage/DamageGradientNodeCoupling.hh"
 #include "Utilities/pointDistances.hh"
 #include "Utilities/DBC.hh"
 #include "Utilities/Timer.hh"
diff --git a/src/Utilities/DamageGradientNodeCoupling.hh b/src/Damage/DamageGradientNodeCoupling.hh
similarity index 100%
rename from src/Utilities/DamageGradientNodeCoupling.hh
rename to src/Damage/DamageGradientNodeCoupling.hh
diff --git a/src/Utilities/DamageGradientNodeCouplingInst.cc.py b/src/Damage/DamageGradientNodeCouplingInst.cc.py
similarity index 85%
rename from src/Utilities/DamageGradientNodeCouplingInst.cc.py
rename to src/Damage/DamageGradientNodeCouplingInst.cc.py
index 250ba43d0..6f3b0a646 100644
--- a/src/Utilities/DamageGradientNodeCouplingInst.cc.py
+++ b/src/Damage/DamageGradientNodeCouplingInst.cc.py
@@ -2,7 +2,7 @@
 //------------------------------------------------------------------------------
 // Explicit instantiation.
 //------------------------------------------------------------------------------
-#include "Utilities/DamageGradientNodeCoupling.cc"
+#include "Damage/DamageGradientNodeCoupling.cc"
 
 namespace Spheral {
   template class DamageGradientNodeCoupling< Dim< %(ndim)s > >;
diff --git a/src/Damage/DamageModel.cc b/src/Damage/DamageModel.cc
index f4320659f..3f7aa9994 100644
--- a/src/Damage/DamageModel.cc
+++ b/src/Damage/DamageModel.cc
@@ -17,13 +17,13 @@
 #include "Boundary/Boundary.hh"
 #include "Kernel/TableKernel.hh"
 #include "Neighbor/ConnectivityMap.hh"
-#include "Utilities/NodeCoupling.hh"
-#include "Utilities/DamagedNodeCoupling.hh"
-#include "Utilities/DamageGradientNodeCoupling.hh"
-#include "Utilities/ThreePointDamagedNodeCoupling.hh"
 #include "Utilities/GeometricUtilities.hh"
 #include "Utilities/safeInv.hh"
 #include "Utilities/Timer.hh"
+#include "Utilities/NodeCoupling.hh"
+#include "Damage/PairMaxDamageNodeCoupling.hh"
+#include "Damage/DamageGradientNodeCoupling.hh"
+#include "Damage/ThreePointDamagedNodeCoupling.hh"
 
 #include "boost/shared_ptr.hpp"
 
@@ -150,7 +150,7 @@ initialize(const Scalar /*time*/,
     break;
 
   case DamageCouplingAlgorithm::PairMaxDamage:
-    mNodeCouplingPtr = std::make_shared<DamagedNodeCoupling<Dimension>>(state, pairs);
+    mNodeCouplingPtr = std::make_shared<PairMaxDamageNodeCoupling<Dimension>>(state, pairs);
     break;
 
   case DamageCouplingAlgorithm::DamageGradient:
diff --git a/src/Damage/DamageModel.hh b/src/Damage/DamageModel.hh
index c8aabf47a..8b81e9ae5 100644
--- a/src/Damage/DamageModel.hh
+++ b/src/Damage/DamageModel.hh
@@ -42,6 +42,7 @@ enum class DamageCouplingAlgorithm {
   PairMaxDamage = 1,
   DamageGradient = 2,
   ThreePointDamage = 3,
+  TensorPairMaxDamage = 4,
 };
 
 template<typename Dimension>
diff --git a/src/Utilities/MaxDamageNodeCoupling.hh b/src/Damage/MaxDamageNodeCoupling.hh
similarity index 100%
rename from src/Utilities/MaxDamageNodeCoupling.hh
rename to src/Damage/MaxDamageNodeCoupling.hh
diff --git a/src/Utilities/DamagedNodeCoupling.cc b/src/Damage/PairMaxDamageNodeCoupling.cc
similarity index 93%
rename from src/Utilities/DamagedNodeCoupling.cc
rename to src/Damage/PairMaxDamageNodeCoupling.cc
index bbc8f070a..61f780757 100644
--- a/src/Utilities/DamagedNodeCoupling.cc
+++ b/src/Damage/PairMaxDamageNodeCoupling.cc
@@ -1,12 +1,12 @@
 //---------------------------------Spheral++----------------------------------//
-// DamagedNodeCoupling
+// PairMaxDamageNodeCoupling
 //
 // A functor class encapsulating how we couple solid nodes in the presence of
 // multiple materials and damage.
 //
 // Created by JMO, Fri Jul 31 14:46:25 PDT 2015
 //----------------------------------------------------------------------------//
-#include "Utilities/DamagedNodeCoupling.hh"
+#include "Damage/PairMaxDamageNodeCoupling.hh"
 #include "Utilities/DBC.hh"
 #include "Field/FieldList.hh"
 #include "Hydro/HydroFieldNames.hh"
@@ -49,8 +49,8 @@ double scalarDamage(const FieldListType& damage,
 // Constructor
 //------------------------------------------------------------------------------
 template<typename Dimension>
-DamagedNodeCoupling<Dimension>::
-DamagedNodeCoupling(const State<Dimension>& state,
+PairMaxDamageNodeCoupling<Dimension>::
+PairMaxDamageNodeCoupling(const State<Dimension>& state,
                     NodePairList& pairs):
   NodeCoupling() {
   const auto n = pairs.size();
diff --git a/src/Utilities/DamagedNodeCoupling.hh b/src/Damage/PairMaxDamageNodeCoupling.hh
similarity index 76%
rename from src/Utilities/DamagedNodeCoupling.hh
rename to src/Damage/PairMaxDamageNodeCoupling.hh
index e8d6a3db5..57f645e46 100644
--- a/src/Utilities/DamagedNodeCoupling.hh
+++ b/src/Damage/PairMaxDamageNodeCoupling.hh
@@ -1,5 +1,5 @@
 //---------------------------------Spheral++----------------------------------//
-// DamagedNodeCoupling
+// PairMaxDamageNodeCoupling
 //
 // A functor class encapsulating how we couple solid nodes in the presence of
 // multiple materials and damage.
@@ -9,8 +9,8 @@
 // Created by JMO: Fri Jul 31 14:46:25 PDT 2015
 //        updated: Fri Feb  5 12:48:20 PST 2021 JMO
 //----------------------------------------------------------------------------//
-#ifndef __Spheral_DamagedNodeCoupling__
-#define __Spheral_DamagedNodeCoupling__
+#ifndef __Spheral_PairMaxDamageNodeCoupling__
+#define __Spheral_PairMaxDamageNodeCoupling__
 
 #include "Utilities/NodeCoupling.hh"
 #include "DataBase/State.hh"
@@ -18,7 +18,7 @@
 namespace Spheral {
 
 template<typename Dimension>
-class DamagedNodeCoupling: public NodeCoupling {
+class PairMaxDamageNodeCoupling: public NodeCoupling {
 public:
   typedef typename Dimension::Scalar Scalar;
   typedef typename Dimension::Vector Vector;
@@ -26,12 +26,12 @@ public:
   typedef typename Dimension::SymTensor SymTensor;
 
   // Constructor.
-  DamagedNodeCoupling(const State<Dimension>& state,
+  PairMaxDamageNodeCoupling(const State<Dimension>& state,
                       NodePairList& pairs);
 
 private:
   // Forbidden methods.
-  DamagedNodeCoupling();
+  PairMaxDamageNodeCoupling();
 };
 
 }
diff --git a/src/Utilities/DamagedNodeCouplingInst.cc.py b/src/Damage/PairMaxDamageNodeCouplingInst.cc.py
similarity index 67%
rename from src/Utilities/DamagedNodeCouplingInst.cc.py
rename to src/Damage/PairMaxDamageNodeCouplingInst.cc.py
index 77b1f5b6a..13b75157f 100644
--- a/src/Utilities/DamagedNodeCouplingInst.cc.py
+++ b/src/Damage/PairMaxDamageNodeCouplingInst.cc.py
@@ -2,9 +2,9 @@
 //------------------------------------------------------------------------------
 // Explicit instantiation.
 //------------------------------------------------------------------------------
-#include "Utilities/DamagedNodeCoupling.cc"
+#include "Damage/PairMaxDamageNodeCoupling.cc"
 
 namespace Spheral {
-  template class DamagedNodeCoupling< Dim< %(ndim)s > >;
+  template class PairMaxDamageNodeCoupling< Dim< %(ndim)s > >;
 }
 """
diff --git a/src/Utilities/ThreePointDamagedNodeCoupling.cc b/src/Damage/ThreePointDamagedNodeCoupling.cc
similarity index 99%
rename from src/Utilities/ThreePointDamagedNodeCoupling.cc
rename to src/Damage/ThreePointDamagedNodeCoupling.cc
index 1edca1d02..47595f1e5 100644
--- a/src/Utilities/ThreePointDamagedNodeCoupling.cc
+++ b/src/Damage/ThreePointDamagedNodeCoupling.cc
@@ -13,7 +13,7 @@
 //
 // Created by JMO, Fri Jan  1 15:18:38 PST 2021
 //----------------------------------------------------------------------------//
-#include "Utilities/ThreePointDamagedNodeCoupling.hh"
+#include "Damage/ThreePointDamagedNodeCoupling.hh"
 #include "Utilities/pointDistances.hh"
 #include "Utilities/DBC.hh"
 #include "Utilities/Timer.hh"
diff --git a/src/Utilities/ThreePointDamagedNodeCoupling.hh b/src/Damage/ThreePointDamagedNodeCoupling.hh
similarity index 100%
rename from src/Utilities/ThreePointDamagedNodeCoupling.hh
rename to src/Damage/ThreePointDamagedNodeCoupling.hh
diff --git a/src/Utilities/ThreePointDamagedNodeCouplingInst.cc.py b/src/Damage/ThreePointDamagedNodeCouplingInst.cc.py
similarity index 84%
rename from src/Utilities/ThreePointDamagedNodeCouplingInst.cc.py
rename to src/Damage/ThreePointDamagedNodeCouplingInst.cc.py
index e7fb3b557..863db8b40 100644
--- a/src/Utilities/ThreePointDamagedNodeCouplingInst.cc.py
+++ b/src/Damage/ThreePointDamagedNodeCouplingInst.cc.py
@@ -2,7 +2,7 @@
 //------------------------------------------------------------------------------
 // Explicit instantiation.
 //------------------------------------------------------------------------------
-#include "Utilities/ThreePointDamagedNodeCoupling.cc"
+#include "Damage/ThreePointDamagedNodeCoupling.cc"
 
 namespace Spheral {
   template class ThreePointDamagedNodeCoupling< Dim< %(ndim)s > >;
diff --git a/src/Pybind11Wraps/Damage/DamageGradientNodeCoupling.py b/src/Pybind11Wraps/Damage/DamageGradientNodeCoupling.py
new file mode 100644
index 000000000..00442d0b7
--- /dev/null
+++ b/src/Pybind11Wraps/Damage/DamageGradientNodeCoupling.py
@@ -0,0 +1,15 @@
+#-------------------------------------------------------------------------------
+# DamageGradientNodeCoupling
+#-------------------------------------------------------------------------------
+from PYB11Generator import *
+import NodeCoupling
+
+@PYB11template("Dimension")
+class DamageGradientNodeCoupling(NodeCoupling.NodeCoupling):
+    """A functor class encapsulating how we couple solid nodes in the presence of
+multiple materials and damage.
+
+This one attempts to mock up the shielding effect of ThreePointDamagedNodeCoupling
+by using local damage gradient to estimate when nodes are separated by
+regions of greater damage (or fractures)."""
+
diff --git a/src/Pybind11Wraps/Damage/DamageMOD.py b/src/Pybind11Wraps/Damage/DamageMOD.py
index ed901795c..bf8cdbb98 100644
--- a/src/Pybind11Wraps/Damage/DamageMOD.py
+++ b/src/Pybind11Wraps/Damage/DamageMOD.py
@@ -13,6 +13,9 @@
 from TensorDamageModel import *
 from ProbabilisticDamageModel import *
 from JohnsonCookDamage import *
+from PairMaxDamageNodeCoupling import *
+from DamageGradientNodeCoupling import *
+from ThreePointDamagedNodeCoupling import *
 
 #-------------------------------------------------------------------------------
 # Includes
@@ -26,6 +29,9 @@
                   '"Damage/weibullFlawDistributionBenzAsphaug.hh"',
                   '"Damage/weibullFlawDistributionOwen.hh"',
                   '"Damage/computeFragmentField.hh"',
+                  '"Damage/PairMaxDamageNodeCoupling.hh"',
+                  '"Damage/ThreePointDamagedNodeCoupling.hh"',
+                  '"Damage/DamageGradientNodeCoupling.hh"',
                   '"FileIO/FileIO.hh"']
 
 #-------------------------------------------------------------------------------
@@ -44,7 +50,8 @@
 DamageCouplingAlgorithm = PYB11enum(("DirectDamage",
                                      "PairMaxDamage",
                                      "DamageGradient",
-                                     "ThreePointDamage"), export_values=True)
+                                     "ThreePointDamage",
+                                     "TensorPairMaxDamage"), export_values=True)
 
 # Note the following enums are deprecated
 EffectiveDamageAlgorithm = PYB11enum(("CopyDamage",
@@ -107,6 +114,10 @@ def computeFragmentField(nodeList = "const NodeList<%(Dimension)s>&",
 ProbabilisticDamageModel%(ndim)id = PYB11TemplateClass(ProbabilisticDamageModel, template_parameters="%(Dimension)s")
 JohnsonCookDamage%(ndim)id = PYB11TemplateClass(JohnsonCookDamage, template_parameters="%(Dimension)s")
 
+PairMaxDamageNodeCoupling%(ndim)id = PYB11TemplateClass(PairMaxDamageNodeCoupling, template_parameters="%(Dimension)s")
+ThreePointDamagedNodeCoupling%(ndim)id = PYB11TemplateClass(ThreePointDamagedNodeCoupling, template_parameters="%(Dimension)s")
+DamageGradientNodeCoupling%(ndim)id = PYB11TemplateClass(DamageGradientNodeCoupling, template_parameters="%(Dimension)s")
+
 weibullFlawDistributionBenzAsphaug%(ndim)id = PYB11TemplateFunction(weibullFlawDistributionBenzAsphaug, template_parameters="%(Dimension)s")
 weibullFlawDistributionOwen%(ndim)id = PYB11TemplateFunction(weibullFlawDistributionOwen, template_parameters="%(Dimension)s")
 computeFragmentField%(ndim)id = PYB11TemplateFunction(computeFragmentField, template_parameters="%(Dimension)s")
diff --git a/src/Pybind11Wraps/Damage/PairMaxDamageNodeCoupling.py b/src/Pybind11Wraps/Damage/PairMaxDamageNodeCoupling.py
new file mode 100644
index 000000000..990282bdb
--- /dev/null
+++ b/src/Pybind11Wraps/Damage/PairMaxDamageNodeCoupling.py
@@ -0,0 +1,25 @@
+#-------------------------------------------------------------------------------
+# PairMaxDamageNodeCoupling
+#-------------------------------------------------------------------------------
+from PYB11Generator import *
+import NodeCoupling
+
+@PYB11template("Dimension")
+class PairMaxDamageNodeCoupling(NodeCoupling.NodeCoupling):
+    """A functor class encapsulating how we couple solid nodes in the presence of
+multiple materials and damage.
+
+This form simply directly damages each pair based on their mutual damage."""
+
+    PYB11typedefs = """
+  typedef typename %(Dimension)s::Scalar Scalar;
+  typedef typename %(Dimension)s::Vector Vector;
+  typedef typename %(Dimension)s::Tensor Tensor;
+  typedef typename %(Dimension)s::SymTensor SymTensor;
+"""
+
+    def pyinit(self,
+               state = "const State<%(Dimension)s>&",
+               pairs = "NodePairList&"):
+        "Constructor"
+
diff --git a/src/Pybind11Wraps/Damage/ThreePointDamagedNodeCoupling.py b/src/Pybind11Wraps/Damage/ThreePointDamagedNodeCoupling.py
new file mode 100644
index 000000000..b1dc3f6c2
--- /dev/null
+++ b/src/Pybind11Wraps/Damage/ThreePointDamagedNodeCoupling.py
@@ -0,0 +1,26 @@
+#-------------------------------------------------------------------------------
+# ThreePointDamagedNodeCoupling
+#-------------------------------------------------------------------------------
+from PYB11Generator import *
+import NodeCoupling
+
+@PYB11template("Dimension")
+class ThreePointDamagedNodeCoupling(NodeCoupling.NodeCoupling):
+    """A functor class encapsulating how we couple solid nodes in the presence of
+multiple materials and damage.
+
+This for uses the "three point" formalism, which allows damaged points to
+cut communication between pairs that talk across them."""
+
+    PYB11typedefs = """
+  typedef typename %(Dimension)s::Scalar Scalar;
+  typedef typename %(Dimension)s::Vector Vector;
+  typedef typename %(Dimension)s::Tensor Tensor;
+  typedef typename %(Dimension)s::SymTensor SymTensor;
+"""
+
+    def pyinit(self,
+               state = "const State<%(Dimension)s>&",
+               W = "const TableKernel<%(Dimension)s>&",
+               pairs = "NodePairList&"):
+        "Constructor"
diff --git a/src/Pybind11Wraps/Utilities/NodeCoupling.py b/src/Pybind11Wraps/Utilities/NodeCoupling.py
index 8e5137134..b651a7976 100644
--- a/src/Pybind11Wraps/Utilities/NodeCoupling.py
+++ b/src/Pybind11Wraps/Utilities/NodeCoupling.py
@@ -18,60 +18,3 @@ def __call__(self,
         "Functional method to override for coupling (nodeListi, i) <-> (nodeListj, j)"
         return "double"
 
-#-------------------------------------------------------------------------------
-# DamagedNodeCoupling
-#-------------------------------------------------------------------------------
-@PYB11template("Dimension")
-class DamagedNodeCoupling(NodeCoupling):
-    """A functor class encapsulating how we couple solid nodes in the presence of
-multiple materials and damage.
-
-This form simply directly damages each pair based on their mutual damage."""
-
-    PYB11typedefs = """
-  typedef typename %(Dimension)s::Scalar Scalar;
-  typedef typename %(Dimension)s::Vector Vector;
-  typedef typename %(Dimension)s::Tensor Tensor;
-  typedef typename %(Dimension)s::SymTensor SymTensor;
-"""
-
-    def pyinit(self,
-               state = "const State<%(Dimension)s>&",
-               pairs = "NodePairList&"):
-        "Constructor"
-
-#-------------------------------------------------------------------------------
-# DamageGradientNodeCoupling
-#-------------------------------------------------------------------------------
-@PYB11template("Dimension")
-class DamageGradientNodeCoupling(NodeCoupling):
-    """A functor class encapsulating how we couple solid nodes in the presence of
-multiple materials and damage.
-
-This one attempts to mock up the shielding effect of ThreePointDamagedNodeCoupling
-by using local damage gradient to estimate when nodes are separated by
-regions of greater damage (or fractures)."""
-
-#-------------------------------------------------------------------------------
-# ThreePointDamagedNodeCoupling
-#-------------------------------------------------------------------------------
-@PYB11template("Dimension")
-class ThreePointDamagedNodeCoupling(NodeCoupling):
-    """A functor class encapsulating how we couple solid nodes in the presence of
-multiple materials and damage.
-
-This for uses the "three point" formalism, which allows damaged points to
-cut communication between pairs that talk across them."""
-
-    PYB11typedefs = """
-  typedef typename %(Dimension)s::Scalar Scalar;
-  typedef typename %(Dimension)s::Vector Vector;
-  typedef typename %(Dimension)s::Tensor Tensor;
-  typedef typename %(Dimension)s::SymTensor SymTensor;
-"""
-
-    def pyinit(self,
-               state = "const State<%(Dimension)s>&",
-               W = "const TableKernel<%(Dimension)s>&",
-               pairs = "NodePairList&"):
-        "Constructor"
diff --git a/src/Pybind11Wraps/Utilities/UtilitiesMOD.py b/src/Pybind11Wraps/Utilities/UtilitiesMOD.py
index 7a3f55a2a..ebcd32632 100644
--- a/src/Pybind11Wraps/Utilities/UtilitiesMOD.py
+++ b/src/Pybind11Wraps/Utilities/UtilitiesMOD.py
@@ -48,9 +48,6 @@
                   '"Utilities/Timer.hh"',
                   '"Utilities/DomainNode.hh"',
                   '"Utilities/NodeCoupling.hh"',
-                  '"Utilities/DamagedNodeCoupling.hh"',
-                  '"Utilities/ThreePointDamagedNodeCoupling.hh"',
-                  '"Utilities/DamageGradientNodeCoupling.hh"',
                   '"Utilities/uniform_random.hh"',
                   '<algorithm>']
 
@@ -258,10 +255,6 @@ def computeShepardsInterpolation(fieldList = "const FieldList<%(Dimension)s, %(D
 # Instantiate stuff for the dimensions Spheral is building
 for ndim in dims:
     exec('''
-DamagedNodeCoupling%(ndim)id = PYB11TemplateClass(DamagedNodeCoupling, template_parameters="%(Dimension)s")
-ThreePointDamagedNodeCoupling%(ndim)id = PYB11TemplateClass(ThreePointDamagedNodeCoupling, template_parameters="%(Dimension)s")
-DamageGradientNodeCoupling%(ndim)id = PYB11TemplateClass(DamageGradientNodeCoupling, template_parameters="%(Dimension)s")
-
 # Functors
 VectorScalarFunctor%(ndim)id = PYB11TemplateClass(SpheralFunctor, template_parameters=("%(Vector)s", "double"))
 VectorVectorFunctor%(ndim)id = PYB11TemplateClass(SpheralFunctor, template_parameters=("%(Vector)s", "%(Vector)s"))
diff --git a/src/Utilities/CMakeLists.txt b/src/Utilities/CMakeLists.txt
index e44a02429..0892afeb7 100644
--- a/src/Utilities/CMakeLists.txt
+++ b/src/Utilities/CMakeLists.txt
@@ -12,9 +12,6 @@ set(Utilities_inst
     medianPosition
     computeShepardsInterpolation
     overlayRemapFields
-    DamagedNodeCoupling
-    ThreePointDamagedNodeCoupling
-    DamageGradientNodeCoupling
    )
 
 set(Utilities_sources
@@ -46,8 +43,6 @@ set(Utilities_headers
     Bessel.hh
     CounterClockwiseComparator.hh
     DBC.hh
-    DamageGradientNodeCoupling.hh
-    DamagedNodeCoupling.hh
     DataTypeTraits.hh
     DomainNode.hh
     DomainNodeInline.hh
@@ -56,7 +51,6 @@ set(Utilities_headers
     Functors.hh
     GeometricUtilities.hh
     KeyTraits.hh
-    MaxDamageNodeCoupling.hh
     NodeCoupling.hh
     NodeListHashMap.hh
     OpenMP_wrapper.hh
@@ -71,7 +65,6 @@ set(Utilities_headers
     SafeIndexMap.hh
     SpheralFunctions.hh
     SurfaceNodeCoupling.hh
-    ThreePointDamagedNodeCoupling.hh
     Timer.hh
     Tree.hh
     TreeInline.hh
diff --git a/src/Utilities/SpheralTimers.cc b/src/Utilities/SpheralTimers.cc
index a5066ea2c..ccec22071 100644
--- a/src/Utilities/SpheralTimers.cc
+++ b/src/Utilities/SpheralTimers.cc
@@ -7,10 +7,10 @@
 #include "Timer.hh"
 #include <list>
 
-// Must initialize the static list defined in Timer.hh
-#ifdef TIMER
-std::list<Timer*> Timer::TimerList(0); 
-#endif
+// // Must initialize the static list defined in Timer.hh
+// #ifdef TIMER
+// std::list<Timer*> Timer::TimerList(0); 
+// #endif
 
 //------------------------------------------------------------------------------
 // Root Timers
diff --git a/src/Utilities/Timer.cc b/src/Utilities/Timer.cc
index 2c4e0e59d..7bf235704 100644
--- a/src/Utilities/Timer.cc
+++ b/src/Utilities/Timer.cc
@@ -19,10 +19,6 @@ using std::abs;
 
 #include "Timer.hh"
 
-// Must initialize the static list defined in Timer.hh
-// list<Timer*> Timer::TimerList(0); 
-
-
 static int ID_counter = 0;
 
 #ifdef PAPI
@@ -74,7 +70,7 @@ Timer::Timer(const string& name):
   //cout << " Root Timer construction for = " << name << endl;
 
   // Make a list of pointers of all Managed Timers.
-  TimerList.push_back(this);
+  TimerRegistrar::timerList().push_back(this);
   
   setup();
 }
@@ -85,7 +81,7 @@ Timer::Timer(const string& name, Timer& parent):
 
   //cout << "   Timer construction for = " << name << endl;
 
-  TimerList.push_back(this);
+  TimerRegistrar::timerList().push_back(this);
 
   setup();
 }
@@ -96,7 +92,7 @@ Timer::Timer(const string& name, Timer& parent, bool):
 
   //cout << "   Diagnostic Timer construction for = " << name << endl;
 
-  TimerList.push_back(this);
+  TimerRegistrar::timerList().push_back(this);
 
   setup();
   diagnostic = true;
@@ -240,7 +236,7 @@ void Timer::TimerSummary( const std::string fname,
   number_procs=1;
 #endif
 
-  // printf("Rank %d, Procs %d, num timers %d\n", rank, number_procs, TimerList.size());
+  // printf("Rank %d, Procs %d, num timers %d\n", rank, number_procs, TimerRegistrar::timerList().size());
 
   list<Timer*>::iterator tli;  
   
@@ -248,7 +244,7 @@ void Timer::TimerSummary( const std::string fname,
 
   // go thru all Timers in static list and obtain the avg/min/max of
   // all processors onto rank0 -- obviously trivial for serial runs.
-  for(tli=TimerList.begin() ; tli != TimerList.end() ; ++tli) {
+  for(tli=TimerRegistrar::timerList().begin() ; tli != TimerRegistrar::timerList().end() ; ++tli) {
       
     // printf(" Timer: %-16s Parent: %-16s diagnostic = %d\n",
     //    (*tli)->Name().c_str(), (*tli)->Parent.Name().c_str(), 
@@ -290,10 +286,10 @@ void Timer::TimerSummary( const std::string fname,
   list<Timer*> ParentList;
   list<Timer*>::iterator pti; 
 
-  //cout << " size of TimerList = " << TimerList.size() << endl;
+  //cout << " size of TimerList = " << TimerRegistrar::timerList().size() << endl;
 
   // for each timer in entire Timer List
-  for(tli=TimerList.begin() ; tli != TimerList.end() ; tli++) {
+  for(tli=TimerRegistrar::timerList().begin() ; tli != TimerRegistrar::timerList().end() ; tli++) {
 
     //printf(" Timer #%3d: %-16s Parent: %-16s diagnostic = %d\n",
     //   (*tli)->ID,
@@ -378,7 +374,7 @@ void Timer::TimerSummary( const std::string fname,
 
       writeSingleLineSeparator(OUT);
 
-      for(tli=TimerList.begin() ; tli != TimerList.end() ; tli++) {
+      for(tli=TimerRegistrar::timerList().begin() ; tli != TimerRegistrar::timerList().end() ; tli++) {
 
         if ((*tli)->diagnostic) continue;
 
@@ -425,7 +421,7 @@ void Timer::TimerSummary( const std::string fname,
   double largestTime = -1.0;
   Timer *DiagnosticParent=NULL;
     
-  for(tli=TimerList.begin() ; tli != TimerList.end() ; tli++) {    
+  for(tli=TimerRegistrar::timerList().begin() ; tli != TimerRegistrar::timerList().end() ; tli++) {    
     if ((*tli)->diagnostic == true) Ndiag++;
     if ((*tli)->avgWC > largestTime) {
       DiagnosticParent = *tli;
@@ -466,7 +462,7 @@ void Timer::TimerSummary( const std::string fname,
     
     writeSingleLineSeparator(OUT);  
     
-    for(tli=TimerList.begin() ; tli != TimerList.end() ; tli++) {
+    for(tli=TimerRegistrar::timerList().begin() ; tli != TimerRegistrar::timerList().end() ; tli++) {
 
       if ((*tli)->diagnostic == false) continue;
 
@@ -824,9 +820,4 @@ static void PAPI_cleanup(int *i) {
 
 #endif // on PAPI
 
-#else
-//#include "Timer.H"
-// Must initialize the static list defined in Timer.hh
-//list<Timer*> Timer::TimerList(0); 
-
 #endif  // TIMER
diff --git a/src/Utilities/Timer.hh b/src/Utilities/Timer.hh
index 2fcf639e4..c4533f6c4 100644
--- a/src/Utilities/Timer.hh
+++ b/src/Utilities/Timer.hh
@@ -12,6 +12,23 @@
 #include <list>
 #include <string>
 
+//------------------------------------------------------------------------------
+// A registrar to hold onto the static list of Timers
+//------------------------------------------------------------------------------
+class Timer;
+class TimerRegistrar {
+public:
+  static std::list<Timer*>& timerList() {
+    if (mTimerListPtr == nullptr) mTimerListPtr = new std::list<Timer*>;
+    return *mTimerListPtr;
+  }
+private:
+  static std::list<Timer*>* mTimerListPtr;
+};
+
+//------------------------------------------------------------------------------
+// Timer
+//------------------------------------------------------------------------------
 // if TIMER is not defined, then we just have a stub class.
 #ifdef TIMER
 
@@ -63,8 +80,6 @@ public:
   
   long int Count() { return count; }
 
-  static std::list<Timer*> TimerList;
-
   static void TimerSummary(const std::string fname = "time.table",
                            const bool printAllTimers = false);
   
@@ -101,6 +116,9 @@ private:
 #else
 
 
+//------------------------------------------------------------------------------
+// Timer (stub)
+//------------------------------------------------------------------------------
 // stub Timer class
 #include <string>
 #include <iostream>
@@ -118,8 +136,6 @@ public:
   Timer(const std::string&, Timer&, bool) {}
   ~Timer() {}
 
-  static std::list<Timer*> TimerList;
-
   inline void setup(){}
   inline void start(){}
   inline void stop(){}