Grab bag of changes for Blast, incl. Axom update

src/CRKSPH/CRKSPH.hh

@@ -8,6 +8,7 @@
 
 #include "CRKSPH/CRKSPHBase.hh"
 #include "Geometry/CellFaceFlag.hh"
+#include "Neighbor/PairwiseField.hh"
 #include "RK/RKCorrectionParams.hh"
 
 #include <memory>
@@ -20,7 +21,6 @@ template<typename Dimension> class TableKernel;
 template<typename Dimension> class DataBase;
 template<typename Dimension, typename Value> class Field;
 template<typename Dimension, typename Value> class FieldList;
-template<typename Dimension, typename Value, size_t numElements> class PairwiseField;
 class FileIO;
 }
 

src/CRKSPH/CRKSPHRZ.hh

@@ -12,6 +12,7 @@
 
 #include "CRKSPH/CRKSPHBase.hh"
 #include "Geometry/Dimension.hh"
+#include "Neighbor/PairwiseField.hh"
 
 namespace Spheral {
 
@@ -22,7 +23,6 @@ template<typename Dimension> class TableKernel;
 template<typename Dimension> class DataBase;
 template<typename Dimension, typename Value> class Field;
 template<typename Dimension, typename Value> class FieldList;
-template<typename Dimension, typename Value, size_t numElements> class PairwiseField;
 class FileIO;
 
 class CRKSPHRZ: public CRKSPHBase<Dim<2>> {

src/FSISPH/SolidFSISPH.cc

@@ -459,6 +459,8 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
   TIME_BEGIN("SolidFSISPHpreStepInitialize");
   if (mApplySelectDensitySum){
   switch(this->densityUpdate()){
+   case FSIMassDensityMethod::FSIConsistentSumMassDensity:
+     // fallthrough intended
    case FSIMassDensityMethod::FSISumMassDensity:
     {
       const auto& W = this->kernel();
@@ -467,7 +469,8 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
       const auto  mass = state.fields(HydroFieldNames::mass, 0.0);
       const auto  H = state.fields(HydroFieldNames::H, SymTensor::zero);
             auto  massDensity = state.fields(HydroFieldNames::massDensity, 0.0);
-      computeFSISPHSumMassDensity(connectivityMap, W, mSumDensityNodeLists, position, mass, H, massDensity);
+      const auto consistentSum = (this->densityUpdate() == FSIMassDensityMethod::FSIConsistentSumMassDensity);
+      computeFSISPHSumMassDensity(connectivityMap, W, mSumDensityNodeLists, position, mass, H, consistentSum, massDensity);
       for (auto boundaryItr = this->boundaryBegin(); boundaryItr < this->boundaryEnd(); ++boundaryItr) (*boundaryItr)->applyFieldListGhostBoundary(massDensity);
       for (auto boundaryItr = this->boundaryBegin(); boundaryItr < this->boundaryEnd(); ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary();
      break;

src/FSISPH/SolidFSISPH.hh

@@ -9,6 +9,7 @@
 #define __Spheral_SolidFSISPH_hh__
 
 #include "Physics/GenericHydro.hh"
+#include "Neighbor/PairwiseField.hh"
 #include "Utilities/SpheralMessage.hh"
 
 #include <string>
@@ -33,6 +34,7 @@ enum class FSIMassDensityMethod {
   FSISumMassDensity = 0,
   PressureCorrectSumMassDensity = 1,
   HWeightedSumMassDensity = 2,
+  FSIConsistentSumMassDensity = 3,
 };
 
 template<typename Dimension> class State;
@@ -43,7 +45,6 @@ template<typename Dimension> class TableKernel;
 template<typename Dimension> class DataBase;
 template<typename Dimension, typename Value> class Field;
 template<typename Dimension, typename Value> class FieldList;
-template<typename Dimension, typename Value, size_t numElements> class PairwiseField;
 class FileIO;
 
 template<typename Dimension>

src/FSISPH/SolidFSISPHEvaluateDerivatives.cc

@@ -560,8 +560,8 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
         if (freeParticle) {
           DvDti += mj*deltaDvDt;
           DvDtj -= mi*deltaDvDt;
-        } 
-      
+        }
+
         // Velocity Gradient
         //-----------------------------------------------------------
         // construct our interface velocity 

src/FSISPH/computeFSISPHSumMassDensity.cc

@@ -14,6 +14,7 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
                             const FieldList<Dimension, typename Dimension::Vector>& position,
                             const FieldList<Dimension, typename Dimension::Scalar>& mass,
                             const FieldList<Dimension, typename Dimension::SymTensor>& H,
+                            const bool consistentSum,
                             FieldList<Dimension, typename Dimension::Scalar>& massDensity) {
 
   // Pre-conditions.
@@ -61,14 +62,15 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
       const auto& ri = position(nodeListi, i);
       const auto& rj = position(nodeListj, j);
       const auto  rij = ri - rj;
+      const auto normalSum = (nodeListi == nodeListj) || consistentSum;
 
       if(sumDensityNodeLists[nodeListi]==1){
         const auto& Hi = H(nodeListi, i);
         const auto  Hdeti = Hi.Determinant();
         const auto  etai = (Hi*rij).magnitude();
         const auto  Wi = W.kernelValue(etai, Hdeti);
 
-        massDensity_thread(nodeListi, i) += (nodeListi == nodeListj ? mj : mi)*Wi; 
+        massDensity_thread(nodeListi, i) += (normalSum ? mj : mi)*Wi; 
       } 
 
       if(sumDensityNodeLists[nodeListj]==1){
@@ -77,7 +79,7 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
         const auto  etaj = (Hj*rij).magnitude();
         const auto  Wj = W.kernelValue(etaj, Hdetj);
 
-        massDensity_thread(nodeListj, j) += (nodeListi == nodeListj ? mi : mj)*Wj;
+        massDensity_thread(nodeListj, j) += (normalSum ? mi : mj)*Wj;
       }
     }
 

src/FSISPH/computeFSISPHSumMassDensity.hh

@@ -22,9 +22,10 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
                             const FieldList<Dimension, typename Dimension::Vector>& position,
                             const FieldList<Dimension, typename Dimension::Scalar>& mass,
                             const FieldList<Dimension, typename Dimension::SymTensor>& H,
+                            const bool consistentSum,
                             FieldList<Dimension, typename Dimension::Scalar>& massDensity);
 
 }
 
 
- #endif
+ #endif

src/FSISPH/computeFSISPHSumMassDensityInst.cc.py

@@ -14,6 +14,7 @@
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>&,
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>&,
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>&,
+                                            const bool consistentSum,
                                                   FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>&);
 }
 """

src/GSPH/GenericRiemannHydro.hh

@@ -13,6 +13,7 @@
 
 #include "Physics/GenericHydro.hh"
 #include "Physics/Physics.hh"
+#include "Neighbor/PairwiseField.hh"
 
 namespace Spheral {
 
@@ -39,7 +40,6 @@ template<typename Dimension> class RiemannSolverBase;
 template<typename Dimension> class DataBase;
 template<typename Dimension, typename Value> class Field;
 template<typename Dimension, typename Value> class FieldList;
-template<typename Dimension, typename Value, size_t numElements> class PairwiseField;
 class FileIO;
 
 template<typename Dimension>

src/GSPH/MFV.hh

@@ -36,6 +36,7 @@
 #include <memory>
 
 #include "GSPH/GenericRiemannHydro.hh"
+#include "Neighbor/PairwiseField.hh"
 
 namespace Spheral {
 
@@ -55,7 +56,6 @@ template<typename Dimension> class RiemannSolverBase;
 template<typename Dimension> class DataBase;
 template<typename Dimension, typename Value> class Field;
 template<typename Dimension, typename Value> class FieldList;
-template<typename Dimension, typename Value, size_t numElements> class PairwiseField;
 class FileIO;
 
 template<typename Dimension>

src/Geometry/GeomPolyhedron.cc

@@ -57,7 +57,6 @@ GeomPolyhedron():
   mCentroid(),
   mRinterior2(-1.0),
   mConvex(true),
-  mSurfaceMeshPtr(nullptr),
   mSurfaceMeshQueryPtr(nullptr),
   mSignedDistancePtr(nullptr) {
   if (mDevnull == NULL) mDevnull = fopen("/dev/null", "w");
@@ -78,7 +77,6 @@ GeomPolyhedron(const vector<GeomPolyhedron::Vector>& points):
   mCentroid(),
   mRinterior2(-1.0),
   mConvex(true),
-  mSurfaceMeshPtr(nullptr),
   mSurfaceMeshQueryPtr(nullptr),
   mSignedDistancePtr(nullptr) {
   TIME_BEGIN("Polyhedron_construct1");
@@ -255,7 +253,6 @@ GeomPolyhedron(const vector<GeomPolyhedron::Vector>& points,
   mCentroid(),
   mRinterior2(-1.0),
   mConvex(false),
-  mSurfaceMeshPtr(nullptr),
   mSurfaceMeshQueryPtr(nullptr),
   mSignedDistancePtr(nullptr) {
   TIME_BEGIN("Polyhedron_construct2");
@@ -290,7 +287,6 @@ GeomPolyhedron(const GeomPolyhedron& rhs):
   mCentroid(rhs.mCentroid),
   mRinterior2(rhs.mRinterior2),
   mConvex(rhs.mConvex),
-  mSurfaceMeshPtr(nullptr),
   mSurfaceMeshQueryPtr(nullptr),
   mSignedDistancePtr(nullptr) {
   for (Facet& facet: mFacets) facet.mVerticesPtr = &mVertices;
@@ -329,7 +325,6 @@ operator=(const GeomPolyhedron& rhs) {
 //------------------------------------------------------------------------------
 GeomPolyhedron::
 ~GeomPolyhedron() {
-  if (mSurfaceMeshPtr != nullptr) delete mSurfaceMeshPtr;
   if (mSurfaceMeshQueryPtr != nullptr) delete mSurfaceMeshQueryPtr;
   if (mSignedDistancePtr != nullptr) delete mSignedDistancePtr;
 }
@@ -347,7 +342,7 @@ contains(const GeomPolyhedron::Vector& point,
 
     // Experimental version using Axom
     using AxPoint = axom::quest::InOutOctree<3>::SpacePt;
-    if (mSurfaceMeshPtr == nullptr) this->buildAxomData();
+    if (!mSurfaceMeshPtr) this->buildAxomData();
     const auto inside = mSurfaceMeshQueryPtr->within(AxPoint(&const_cast<Vector&>(point)[0]));
     if (not inside and countBoundary) {
       return this->distance(point) < tol;
@@ -747,7 +742,7 @@ distance(const GeomPolyhedron::Vector& p,
 
     // Experimental version using Axom
     using AxPoint = axom::quest::InOutOctree<3>::SpacePt;
-    if (mSurfaceMeshPtr == nullptr) this->buildAxomData();
+    if (!mSurfaceMeshPtr) this->buildAxomData();
     return std::abs(mSignedDistancePtr->computeDistance(AxPoint(&const_cast<Vector&>(p)[0])));
 
   } else {
@@ -932,10 +927,9 @@ setBoundingBox() {
   TIME_END("Polyhedron_BB_R2");
 
   // Clear any existing Axom information, so it's reconstructed if needed
-  if (mSurfaceMeshPtr != nullptr) delete mSurfaceMeshPtr;
   if (mSurfaceMeshQueryPtr != nullptr) delete mSurfaceMeshQueryPtr;
   if (mSignedDistancePtr != nullptr) delete mSignedDistancePtr;
-  mSurfaceMeshPtr = nullptr;
+  mSurfaceMeshPtr.reset();
   mSurfaceMeshQueryPtr = nullptr;
   mSignedDistancePtr = nullptr;
   TIME_END("Polyhedron_BB");
@@ -979,10 +973,10 @@ buildAxomData() const {
   bb.addPoint(AxPoint(&xmin[0]));
   bb.addPoint(AxPoint(&xmax[0]));
   axom::mint::write_vtk(meshPtr, "blago.vtk");
-  mSurfaceMeshPtr = meshPtr;
+  mSurfaceMeshPtr = std::shared_ptr<axom::quest::InOutOctree<3>::SurfaceMesh>(meshPtr);
   mSurfaceMeshQueryPtr = new AxOctree(bb, mSurfaceMeshPtr);
   mSurfaceMeshQueryPtr->generateIndex();
-  mSignedDistancePtr = new AxDistance(mSurfaceMeshPtr,
+  mSignedDistancePtr = new AxDistance(mSurfaceMeshPtr.get(),
                                       true);               // is_watertight
 }
 

src/Geometry/GeomPolyhedron.hh

@@ -154,7 +154,7 @@ private:
   Vector mXmin, mXmax, mCentroid;
   double mRinterior2;
   bool mConvex;
-  mutable axom::quest::InOutOctree<3>::SurfaceMesh* mSurfaceMeshPtr;
+  mutable std::shared_ptr<axom::quest::InOutOctree<3>::SurfaceMesh> mSurfaceMeshPtr;
   mutable axom::quest::InOutOctree<3>* mSurfaceMeshQueryPtr;
   mutable axom::quest::SignedDistance<3>* mSignedDistancePtr;
 

src/Neighbor/ConnectivityMap.cc

@@ -8,6 +8,7 @@
 #include "ConnectivityMap.hh"
 #include "NodeList/NodeList.hh"
 #include "Neighbor/Neighbor.hh"
+#include "Neighbor/NodePairList.hh"
 #include "DataBase/DataBase.hh"
 #include "Field/FieldList.hh"
 #include "Boundary/Boundary.hh"
@@ -134,7 +135,8 @@ ConnectivityMap():
   mNodeTraversalIndices(),
   mKeys(FieldStorageType::CopyFields),
   mCouplingPtr(std::make_shared<NodeCoupling>()),
-  mIntersectionConnectivity() {
+  mIntersectionConnectivity(),
+  mExcludePairs([](int, int, int, int) { return false; }){
 }
 
 //------------------------------------------------------------------------------
@@ -286,6 +288,7 @@ patchConnectivity(const FieldList<Dimension, int>& flags,
       culledPairs.insert(culledPairs.end(), culledPairs_thread.begin(), culledPairs_thread.end());
     }
   }
+
   mNodePairListPtr = culledPairListPtr;
 
   // Sort the NodePairList in order to enforce domain decomposition independence.
@@ -930,9 +933,14 @@ computeConnectivity() {
 
                     // We don't include self-interactions.
                     if ((iNodeList != jNodeList) or (i != j)) {
-                      neighbors[jNodeList].push_back(j);
-                      if (calculatePairInteraction(iNodeList, i, jNodeList, j, firstGhostNodej)) nodePairs_private.push_back(NodePairIdxType(i, iNodeList, j, jNodeList));
-                      if (domainDecompIndependent) keys[jNodeList].push_back(pair<int, Key>(j, mKeys(jNodeList, j)));
+
+                      // Exclude specified pairs.
+                      if (!mExcludePairs(iNodeList, i, jNodeList, j))
+                      {
+                        neighbors[jNodeList].push_back(j);
+                        if (calculatePairInteraction(iNodeList, i, jNodeList, j, firstGhostNodej)) nodePairs_private.push_back(NodePairIdxType(i, iNodeList, j, jNodeList));
+                        if (domainDecompIndependent) keys[jNodeList].push_back(pair<int, Key>(j, mKeys(jNodeList, j)));
+                      }
                     }
                   }
                 }
@@ -1139,6 +1147,8 @@ computeConnectivity() {
   //   }
   // }
 
+  mNodePairListPtr->makeUnique();
+
   // Post conditions.
   BEGIN_CONTRACT_SCOPE
   // Make sure that the correct number of nodes have been completed.
@@ -1158,4 +1168,53 @@ computeConnectivity() {
   TIME_END("ConnectivityMap_computeConnectivity");
 }
 
+//------------------------------------------------------------------------------
+// Remove given pairs from the connectivity in place
+//------------------------------------------------------------------------------
+template<typename Dimension>
+void
+ConnectivityMap<Dimension>::
+removeConnectivity(std::function<bool(int, int, int, int)> excludePairs)
+{
+  const auto domainDecompIndependent = NodeListRegistrar<Dimension>::instance().domainDecompositionIndependent();
+  const auto numNodeLists = mNodeLists.size();
+  for (auto iNodeList = 0u; iNodeList != numNodeLists; ++iNodeList) {
+    const auto numNodes = ((domainDecompIndependent or mBuildGhostConnectivity or mBuildOverlapConnectivity) ? 
+                           mNodeLists[iNodeList]->numNodes() :
+                           mNodeLists[iNodeList]->numInternalNodes());
+    for (auto i = 0u; i < numNodes; ++i) {
+      auto& neighbors = mConnectivity[i];
+      CHECK(neighbors.size() == numNodeLists);
+      for (auto jNodeList = 0u; jNodeList < numNodeLists; ++jNodeList) {
+        auto& neighborsj = neighbors[jNodeList];
+        auto l = 0u;
+        for (auto k = 0u; k < neighborsj.size(); ++k)
+        {
+          const auto j = neighborsj[k];
+          if (!excludePairs(iNodeList, i, jNodeList, j))
+          {
+            neighborsj[l++] = neighborsj[k];
+          }
+        }
+        neighborsj.resize(l);
+      }
+    }
+  }
+
+  auto& nodePairList = *mNodePairListPtr;
+  const auto npairs = nodePairList.size();
+  auto l = 0u;
+  for (auto k = 0u; k < npairs; ++k) {
+    const auto iNodeList = nodePairList[k].i_list;
+    const auto jNodeList = nodePairList[k].j_list;
+    const auto i = nodePairList[k].i_node;
+    const auto j = nodePairList[k].j_node;
+    if (!excludePairs(iNodeList, i, jNodeList, j))
+    {
+      nodePairList[l++] = nodePairList[k];
+    }
+  }
+  nodePairList.resize(l);
+}
+
 }