Skip to content

Grab bag of changes for Blast, incl. Axom update #333

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
brbass wants to merge 14 commits into
base: develop
Choose a base branch
from
Draft

Grab bag of changes for Blast, incl. Axom update #333

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
2500c87
Added optional node pair exclusion
Feb 13, 2024
375600d
Test functions for ensuring that node pairs are unique.
Feb 21, 2024
4220f38
Switched location of culling.
Mar 27, 2024
b990b4d
Added function to integrator for culling
Oct 28, 2024
b8b73a6
Revert "Added function to integrator for culling"
Oct 29, 2024
679d504
Ignore connectivity in FSISPH if f_couple=0
Dec 4, 2024
02119b0
Tools for density update
Feb 18, 2025
5a962ac
Merge branch 'develop' into feature/cull-pairs
May 16, 2025
4c8518f
Removed couple check
May 19, 2025
fe24007
Axom changes
May 19, 2025
c22bde7
Implicit delete wasn't working
May 19, 2025
93464f0
Wrong bracket
May 21, 2025
b087c08
Formatting
May 21, 2025
90b7fe9
Requested changes
May 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions src/CRKSPH/CRKSPH.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,15 @@ CRKSPH(DataBase<Dimension>& dataBase,
mPairAccelerationsPtr() {
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
template<typename Dimension>
CRKSPH<Dimension>::
~CRKSPH(){
// Needs to be here due to implicit PairwiseField delete
}

Copy link
Copy Markdown
Contributor Author

@brbass brbass May 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Had to add this to all the hydro packages. Otherwise, the build fails since the PairwiseField destructor isn't available in the header file.

//------------------------------------------------------------------------------
// Register the state we need/are going to evolve.
//------------------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion src/CRKSPH/CRKSPH.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ public:
CRKSPH& operator=(const CRKSPH&) = delete;

// Destructor.
virtual ~CRKSPH() = default;
virtual ~CRKSPH();

// Register the state Hydro expects to use and evolve.
virtual
Expand Down
8 changes: 8 additions & 0 deletions src/CRKSPH/CRKSPHRZ.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,14 @@ CRKSPHRZ(DataBase<Dimension>& dataBase,
mPairAccelerationsPtr() {
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
CRKSPHRZ::
~CRKSPHRZ(){
// Needs to be here due to implicit PairwiseField delete
}

//------------------------------------------------------------------------------
// On problem start up some dependent state needs to be calculated
//------------------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion src/CRKSPH/CRKSPHRZ.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ public:
CRKSPHRZ& operator=(const CRKSPHRZ&) = delete;

// Destructor.
virtual ~CRKSPHRZ() = default;
virtual ~CRKSPHRZ();

// A second optional method to be called on startup, after Physics::initializeProblemStartup has
// been called.
Expand Down
9 changes: 9 additions & 0 deletions src/CRKSPH/SolidCRKSPH.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,15 @@ SolidCRKSPH(DataBase<Dimension>& dataBase,
mPlasticStrain0 = dataBase.newSolidFieldList(0.0, SolidFieldNames::plasticStrain + "0");
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
template<typename Dimension>
SolidCRKSPH<Dimension>::
~SolidCRKSPH(){
// Needs to be here due to implicit PairwiseField delete
}

//------------------------------------------------------------------------------
// On problem start up, we need to initialize our internal data.
//------------------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion src/CRKSPH/SolidCRKSPH.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ public:
SolidCRKSPH& operator=(const SolidCRKSPH&) = delete;

// Destructor.
virtual ~SolidCRKSPH() = default;
virtual ~SolidCRKSPH();

// Tasks we do once on problem startup.
virtual
Expand Down
8 changes: 8 additions & 0 deletions src/CRKSPH/SolidCRKSPHRZ.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,14 @@ SolidCRKSPHRZ(DataBase<Dimension>& dataBase,
mSelfAccelerations(FieldStorageType::CopyFields) {
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
SolidCRKSPHRZ::
~SolidCRKSPHRZ(){
// Needs to be here due to implicit PairwiseField delete
}

//------------------------------------------------------------------------------
// On problem start up, we need to initialize our internal data.
//------------------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion src/CRKSPH/SolidCRKSPHRZ.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ public:
SolidCRKSPHRZ& operator=(const SolidCRKSPHRZ&) = delete;

// Destructor.
virtual ~SolidCRKSPHRZ() = default;
virtual ~SolidCRKSPHRZ();

// A second optional method to be called on startup, after Physics::initializeProblemStartup has
// been called.
Expand Down
14 changes: 13 additions & 1 deletion src/FSISPH/SolidFSISPH.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,6 +241,15 @@ SolidFSISPH(DataBase<Dimension>& dataBase,
mInterfaceAngles = dataBase.newFluidFieldList(0.0, FSIFieldNames::interfaceAngles);
}

//------------------------------------------------------------------------------
// Destructor.
//------------------------------------------------------------------------------
template<typename Dimension>
SolidFSISPH<Dimension>::
~SolidFSISPH(){
// Needs to be here due to implicit PairwiseField delete
}

//------------------------------------------------------------------------------
// On problem start up, we need to initialize our internal data.
//------------------------------------------------------------------------------
Expand Down Expand Up @@ -459,6 +468,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();
Expand All @@ -467,7 +478,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;
Expand Down
3 changes: 2 additions & 1 deletion src/FSISPH/SolidFSISPH.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ enum class FSIMassDensityMethod {
FSISumMassDensity = 0,
PressureCorrectSumMassDensity = 1,
HWeightedSumMassDensity = 2,
FSIConsistentSumMassDensity = 3,
};

template<typename Dimension> class State;
Expand Down Expand Up @@ -93,7 +94,7 @@ public:
SolidFSISPH(const SolidFSISPH&) = delete;
SolidFSISPH& operator=(const SolidFSISPH&) = delete;

virtual ~SolidFSISPH() = default;
virtual ~SolidFSISPH();

// A second optional method to be called on startup, after Physics::initializeProblemStartup has
// been called.
Expand Down
4 changes: 2 additions & 2 deletions src/FSISPH/SolidFSISPHEvaluateDerivatives.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -560,8 +560,8 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
if (freeParticle) {
DvDti += mj*deltaDvDt;
DvDtj -= mi*deltaDvDt;
}
}

// Velocity Gradient
//-----------------------------------------------------------
// construct our interface velocity
Expand Down
6 changes: 4 additions & 2 deletions src/FSISPH/computeFSISPHSumMassDensity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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){
Expand All @@ -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;
}
}

Expand Down
3 changes: 2 additions & 1 deletion src/FSISPH/computeFSISPHSumMassDensity.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions src/FSISPH/computeFSISPHSumMassDensityInst.cc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>&);
}
"""
16 changes: 5 additions & 11 deletions src/Geometry/GeomPolyhedron.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -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");
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand All @@ -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;
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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
}

Expand Down
2 changes: 1 addition & 1 deletion src/Geometry/GeomPolyhedron.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
Loading