diff --git a/CMakeLists.txt b/CMakeLists.txt index e456bb3b..773b9757 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -160,7 +160,7 @@ endif() set(TRIBOL_DEBUG_DEFINE_STRING "${TRIBOL_DEBUG_DEFINE_STRING}" CACHE STRING "" FORCE) # Define AXOM_DEBUG for SLIC macros -string(REPLACE "TRIBOL" "AXOM" AXOM_DEBUG_DEFINE_STRING TRIBOL_DEBUG_DEFINE_STRING) +string(REPLACE "TRIBOL" "AXOM" AXOM_DEBUG_DEFINE_STRING "${TRIBOL_DEBUG_DEFINE_STRING}") set(AXOM_DEBUG_DEFINE_STRING "${AXOM_DEBUG_DEFINE_STRING}" CACHE STRING "" FORCE) #------------------------------------------------------------------------------ diff --git a/README.md b/README.md index fafbc263..a64ff734 100644 --- a/README.md +++ b/README.md @@ -71,7 +71,7 @@ to build Tribol. The Tribol contact physics library requires: - CMake 3.14 or higher -- C++14 compiler +- C++17 compiler - MPI - mfem - axom diff --git a/scripts/spack/packages/mfem/package.py b/scripts/spack/packages/mfem/package.py index 0e808959..aadaba4f 100644 --- a/scripts/spack/packages/mfem/package.py +++ b/scripts/spack/packages/mfem/package.py @@ -11,7 +11,7 @@ class Mfem(BuiltinMfem): # Note: Make sure this sha coincides with the git submodule # Note: We add a number to the end of the real version number to indicate that we have # moved forward past the release. Increment the last number when updating the commit sha. - version("4.8.0.1", commit="d9c1c34fdfaf3f7a9f56dfc82f7c083082a36fca") + version("4.9.0", commit="d9d6526cc1749980a2ba1da16e2c1ca1e07d82ec") variant('asan', default=False, description='Add Address Sanitizer flags') diff --git a/scripts/spack/packages/tribol/package.py b/scripts/spack/packages/tribol/package.py index 78153b1c..514fcada 100644 --- a/scripts/spack/packages/tribol/package.py +++ b/scripts/spack/packages/tribol/package.py @@ -81,6 +81,7 @@ class Tribol(CachedCMakePackage, CudaPackage, ROCmPackage): # Other libraries depends_on("mfem@4.7.0.2:+lapack") + depends_on("mfem@4.9.0:+lapack", when="+enzyme") depends_on("axom@0.9:") depends_on("raja@2024.02.0:", when="+raja") @@ -95,6 +96,8 @@ class Tribol(CachedCMakePackage, CudaPackage, ROCmPackage): depends_on("mfem+metis+mpi", when="+redecomp") depends_on("mfem+asan", when="+asan") + # Tribol uses MFEM's enzyme header + depends_on("mfem+enzyme", when="+enzyme") with when("+profiling"): depends_on("caliper+mpi") diff --git a/src/examples/examples_common.hpp b/src/examples/examples_common.hpp index 4c6f5d81..bd4b75d4 100644 --- a/src/examples/examples_common.hpp +++ b/src/examples/examples_common.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef EXAMPLES_COMMON_HPP_ -#define EXAMPLES_COMMON_HPP_ +#ifndef SRC_EXAMPLES_EXAMPLES_COMMON_HPP_ +#define SRC_EXAMPLES_EXAMPLES_COMMON_HPP_ // axom includes #include "axom/core.hpp" @@ -425,4 +425,4 @@ void initialize_logger( tribol::CommT problem_comm ) */ void finalize_logger() { slic::finalize(); } -#endif /* EXAMPLES_COMMON_HPP_ */ +#endif /* SRC_EXAMPLES_EXAMPLES_COMMON_HPP_ */ diff --git a/src/examples/mortar_lm_patch_test.cpp b/src/examples/mortar_lm_patch_test.cpp index 9fa6cd16..384622a0 100644 --- a/src/examples/mortar_lm_patch_test.cpp +++ b/src/examples/mortar_lm_patch_test.cpp @@ -8,13 +8,7 @@ // Tribol includes #include "tribol/interface/tribol.hpp" #include "tribol/utils/TestUtils.hpp" -#include "tribol/utils/Math.hpp" #include "tribol/common/Parameters.hpp" -#include "tribol/mesh/MethodCouplingData.hpp" -#include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/mesh/InterfacePairs.hpp" -#include "tribol/mesh/MeshData.hpp" -#include "tribol/geom/GeomUtilities.hpp" // Axom includes #include "axom/slic.hpp" @@ -25,9 +19,7 @@ #endif // C/C++ includes -#include // for std::sin() and std::cos() -#include // for std::string and operators -#include // for std::ostringstream +#include // for std::string and operators //---------------------------------------------------------------------------------- // Example command line arguments for running this example. This test creates two @@ -196,10 +188,8 @@ int main( int argc, char** argv ) // instantiate mfem vector for RHS contributions int rhs_size = mesh.dim * mesh.numTotalNodes + // equilibrium equations mesh.numNonmortarSurfaceNodes; // gap equations - RealT b[rhs_size]; - tribol::initRealArray( &b[0], rhs_size, 0. ); - mfem::Vector rhs( &b[0], rhs_size ); - rhs = 0.; // initialize + mfem::Vector rhs( rhs_size ); + rhs = 0.0; // initialize SLIC_INFO( "Finalized initial oversized sparse matrix and created rhs vector." ); @@ -230,7 +220,7 @@ int main( int argc, char** argv ) ///////////////////////////////////////// // populate RHS with gap contributions // ///////////////////////////////////////// - mesh.getGapEvals( &b[0] ); + mesh.getGapEvals( rhs.GetData() ); SLIC_INFO( "Populated RHS gap contributions." ); diff --git a/src/redecomp/transfer/GridFnTransfer.hpp b/src/redecomp/transfer/GridFnTransfer.hpp index cad1cd51..0763b8f3 100644 --- a/src/redecomp/transfer/GridFnTransfer.hpp +++ b/src/redecomp/transfer/GridFnTransfer.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_REDECOMP_GRIDFNTRANSFER_HPP_ -#define SRC_REDECOMP_GRIDFNTRANSFER_HPP_ +#ifndef SRC_REDECOMP_TRANSFER_GRIDFNTRANSFER_HPP_ +#define SRC_REDECOMP_TRANSFER_GRIDFNTRANSFER_HPP_ #include "mfem.hpp" @@ -43,4 +43,4 @@ class GridFnTransfer { } // end namespace redecomp -#endif /* SRC_REDECOMP_GRIDFNTRANSFER_HPP_ */ +#endif /* SRC_REDECOMP_TRANSFER_GRIDFNTRANSFER_HPP_ */ diff --git a/src/redecomp/transfer/MatrixTransfer.hpp b/src/redecomp/transfer/MatrixTransfer.hpp index 7e4f03de..484699ce 100644 --- a/src/redecomp/transfer/MatrixTransfer.hpp +++ b/src/redecomp/transfer/MatrixTransfer.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_REDECOMP_MATRIXTRANSFER_HPP_ -#define SRC_REDECOMP_MATRIXTRANSFER_HPP_ +#ifndef SRC_REDECOMP_TRANSFER_MATRIXTRANSFER_HPP_ +#define SRC_REDECOMP_TRANSFER_MATRIXTRANSFER_HPP_ #include "mfem.hpp" @@ -202,4 +202,4 @@ class MatrixTransfer { } // end namespace redecomp -#endif /* SRC_REDECOMP_MATRIXTRANSFER_HPP_ */ +#endif /* SRC_REDECOMP_TRANSFER_MATRIXTRANSFER_HPP_ */ diff --git a/src/redecomp/transfer/SparseMatrixTransfer.hpp b/src/redecomp/transfer/SparseMatrixTransfer.hpp index 8406a708..9090d502 100644 --- a/src/redecomp/transfer/SparseMatrixTransfer.hpp +++ b/src/redecomp/transfer/SparseMatrixTransfer.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_REDECOMP_SPARSEMATRIXTRANSFER_HPP_ -#define SRC_REDECOMP_SPARSEMATRIXTRANSFER_HPP_ +#ifndef SRC_REDECOMP_TRANSFER_SPARSEMATRIXTRANSFER_HPP_ +#define SRC_REDECOMP_TRANSFER_SPARSEMATRIXTRANSFER_HPP_ #include @@ -113,4 +113,4 @@ class SparseMatrixTransfer { } // end namespace redecomp -#endif /* SRC_REDECOMP_SPARSEMATRIXTRANSFER_HPP_ */ +#endif /* SRC_REDECOMP_TRANSFER_SPARSEMATRIXTRANSFER_HPP_ */ diff --git a/src/redecomp/transfer/TransferByElements.hpp b/src/redecomp/transfer/TransferByElements.hpp index f6548777..d1017b62 100644 --- a/src/redecomp/transfer/TransferByElements.hpp +++ b/src/redecomp/transfer/TransferByElements.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_REDECOMP_TRANSFERBYELEMENTS_HPP_ -#define SRC_REDECOMP_TRANSFERBYELEMENTS_HPP_ +#ifndef SRC_REDECOMP_TRANSFER_TRANSFERBYELEMENTS_HPP_ +#define SRC_REDECOMP_TRANSFER_TRANSFERBYELEMENTS_HPP_ #include "redecomp/transfer/GridFnTransfer.hpp" @@ -47,4 +47,4 @@ class TransferByElements : public GridFnTransfer { } // end namespace redecomp -#endif /* SRC_REDECOMP_TRANSFERBYELEMENTS_HPP_ */ +#endif /* SRC_REDECOMP_TRANSFER_TRANSFERBYELEMENTS_HPP_ */ diff --git a/src/redecomp/transfer/TransferByNodes.hpp b/src/redecomp/transfer/TransferByNodes.hpp index dedf91c8..417096c7 100644 --- a/src/redecomp/transfer/TransferByNodes.hpp +++ b/src/redecomp/transfer/TransferByNodes.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_REDECOMP_TRANSFERBYNODES_HPP_ -#define SRC_REDECOMP_TRANSFERBYNODES_HPP_ +#ifndef SRC_REDECOMP_TRANSFER_TRANSFERBYNODES_HPP_ +#define SRC_REDECOMP_TRANSFER_TRANSFERBYNODES_HPP_ #include "mfem.hpp" @@ -120,4 +120,4 @@ class TransferByNodes : public GridFnTransfer { } // end namespace redecomp -#endif /* SRC_REDECOMP_TRANSFERBYNODES_HPP_ */ +#endif /* SRC_REDECOMP_TRANSFER_TRANSFERBYNODES_HPP_ */ diff --git a/src/redecomp/utils/BisecTree.hpp b/src/redecomp/utils/BisecTree.hpp index 2c0bb8d4..b7b78db5 100644 --- a/src/redecomp/utils/BisecTree.hpp +++ b/src/redecomp/utils/BisecTree.hpp @@ -331,4 +331,4 @@ class BisecTree { } // end namespace redecomp -#endif /* SRC_REDECOMP_BISECTREE_HPP_ */ +#endif /* SRC_REDECOMP_UTILS_BISECTREE_HPP_ */ diff --git a/src/shared/mesh/MeshBuilder.hpp b/src/shared/mesh/MeshBuilder.hpp index b2a3c2f0..c95aa252 100644 --- a/src/shared/mesh/MeshBuilder.hpp +++ b/src/shared/mesh/MeshBuilder.hpp @@ -1,5 +1,5 @@ -#ifndef SRC_SHARED_MESHBUILDER_HPP_ -#define SRC_SHARED_MESHBUILDER_HPP_ +#ifndef SRC_SHARED_MESH_MESHBUILDER_HPP_ +#define SRC_SHARED_MESH_MESHBUILDER_HPP_ #include @@ -220,4 +220,4 @@ class ParMeshBuilder { } // namespace shared -#endif // SRC_SHARED_MESHBUILDER_HPP_ +#endif // SRC_SHARED_MESH_MESHBUILDER_HPP_ diff --git a/src/tests/tribol_common_plane_gpu.cpp b/src/tests/tribol_common_plane_gpu.cpp index c8226822..848dd4fe 100644 --- a/src/tests/tribol_common_plane_gpu.cpp +++ b/src/tests/tribol_common_plane_gpu.cpp @@ -3,6 +3,8 @@ // // SPDX-License-Identifier: (MIT) +#include "tribol/config.hpp" + #include "gtest/gtest.h" #include "mfem.hpp" diff --git a/src/tests/tribol_comp_geom_3d.cpp b/src/tests/tribol_comp_geom_3d.cpp index 151a76ab..5570442d 100644 --- a/src/tests/tribol_comp_geom_3d.cpp +++ b/src/tests/tribol_comp_geom_3d.cpp @@ -55,7 +55,7 @@ class CompGeomTest : public ::testing::Test { { // if element thickness is not registered by test, then register dummy // element thickness in order to use auto contact - RealT element_thickness[numCells]; + tribol::Array1D element_thickness( numCells ); if ( !m_isElementThicknessRegistered ) { for ( int i = 0; i < numCells; ++i ) { element_thickness[i] = 1.0; diff --git a/src/tests/tribol_coupling_scheme.cpp b/src/tests/tribol_coupling_scheme.cpp index 03b21c73..82f5985d 100644 --- a/src/tests/tribol_coupling_scheme.cpp +++ b/src/tests/tribol_coupling_scheme.cpp @@ -8,30 +8,16 @@ #include "tribol/utils/TestUtils.hpp" #include "tribol/utils/Math.hpp" #include "tribol/common/Parameters.hpp" -#include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/mesh/InterfacePairs.hpp" -#include "tribol/mesh/MeshData.hpp" -#include "tribol/geom/GeomUtilities.hpp" #ifdef TRIBOL_USE_UMPIRE // Umpire includes #include "umpire/ResourceManager.hpp" #endif -// Axom includes -#include "axom/slic.hpp" - // gtest includes #include "gtest/gtest.h" -// c++ includes -#include // std::abs, std::cos, std::sin -#include -#include -#include -#include - using RealT = tribol::RealT; /*! @@ -211,8 +197,8 @@ TEST_F( CouplingSchemeTest, single_mortar_2D ) // register dummy nodal fields so error doesn't return from field // registration - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); @@ -240,8 +226,8 @@ TEST_F( CouplingSchemeTest, aligned_mortar_2D ) // register dummy nodal fields so error doesn't return from field // registration - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); @@ -290,8 +276,8 @@ TEST_F( CouplingSchemeTest, single_mortar_3D_penalty ) // register dummy nodal fields so error doesn't return from field // registration - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); @@ -358,8 +344,8 @@ TEST_F( CouplingSchemeTest, mortar_tied ) // register dummy nodal fields so error doesn't return from field // registration - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); @@ -387,8 +373,8 @@ TEST_F( CouplingSchemeTest, mortar_coulomb ) // register dummy nodal fields so error doesn't return from field // registration - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); @@ -818,8 +804,8 @@ TEST_F( CouplingSchemeTest, single_mortar_null_response_pointers ) registerDummy3DMesh( 0, numCells, setResponse ); registerDummy3DMesh( 1, numCells, setResponse ); - RealT gaps[this->m_lengthNodalData]; - RealT pressures[this->m_lengthNodalData]; + tribol::Array1D gaps( this->m_lengthNodalData ); + tribol::Array1D pressures( this->m_lengthNodalData ); tribol::registerMortarGaps( 1, &gaps[0] ); tribol::registerMortarPressures( 1, &pressures[0] ); diff --git a/src/tests/tribol_enzyme_nodal_normal.cpp b/src/tests/tribol_enzyme_nodal_normal.cpp index d8ffda54..b399fcbd 100644 --- a/src/tests/tribol_enzyme_nodal_normal.cpp +++ b/src/tests/tribol_enzyme_nodal_normal.cpp @@ -67,7 +67,7 @@ class EnzymeNodalNormalTest : public testing::Test { mfem::DenseMatrix dndx_fd( num_dofs ); auto mesh_view = mesh.getView(); - Array2D n_base = mesh_view.getNodalNormals(); + tribol::Array2D n_base = mesh_view.getNodalNormals(); for ( int dx{ 0 }; dx < mesh.spatialDimension(); ++dx ) { for ( int nx{ 0 }; nx < mesh.numberOfNodes(); ++nx ) { auto x_idx = dx * mesh.numberOfNodes() + nx; diff --git a/src/tests/tribol_iso_integ.cpp b/src/tests/tribol_iso_integ.cpp index fd750520..d34fdb23 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -3,22 +3,19 @@ // // SPDX-License-Identifier: (MIT) +// c++ includes +#include // std::abs + +// gtest includes +#include "gtest/gtest.h" + // Tribol includes -#include "tribol/mesh/MeshData.hpp" +#include "tribol/common/ArrayTypes.hpp" #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/integ/Integration.hpp" #include "tribol/geom/GeomUtilities.hpp" #include "tribol/integ/FE.hpp" -// Axom includes -#include "axom/slic.hpp" - -// gtest includes -#include "gtest/gtest.h" - -// c++ includes -#include // std::abs - using RealT = tribol::RealT; /*! @@ -34,7 +31,7 @@ using RealT = tribol::RealT; class IsoIntegTest : public ::testing::Test { public: int numNodes; - int dim; + static constexpr int dim = 3; RealT* getXCoords() { return x; } @@ -44,34 +41,20 @@ class IsoIntegTest : public ::testing::Test { bool integrate( RealT const tol ) { - RealT xyz[this->dim * this->numNodes]; - RealT* xy = xyz; - - RealT* x = this->x; - RealT* y = this->y; - RealT* z = this->z; + tribol::Array2D xyz( this->numNodes, dim ); // generate stacked coordinate array for ( int j = 0; j < this->numNodes; ++j ) { - for ( int k = 0; k < this->dim; ++k ) { - switch ( k ) { - case 0: - xy[this->dim * j + k] = x[j]; - break; - case 1: - xy[this->dim * j + k] = y[j]; - break; - case 2: - xy[this->dim * j + k] = z[j]; - break; - } // end switch - } // end loop over dimension + xyz( j, 0 ) = x[j]; + xyz( j, 1 ) = y[j]; + xyz( j, 2 ) = z[j]; } // end loop over nodes // instantiate SurfaceContactElem struct. Note, this object is instantiated // using face 1 as face 2, but these faces are not used in this test so this // is ok. - tribol::SurfaceContactElem elem( this->dim, xy, xy, xy, this->numNodes, this->numNodes, nullptr, nullptr, 0, 0 ); + tribol::SurfaceContactElem elem( this->dim, xyz.data(), xyz.data(), xyz.data(), this->numNodes, this->numNodes, + nullptr, nullptr, 0, 0 ); // instantiate integration object tribol::IntegPts integ; @@ -108,7 +91,6 @@ class IsoIntegTest : public ::testing::Test { void SetUp() override { this->numNodes = 4; - this->dim = 3; if ( this->x == nullptr ) { this->x = new RealT[this->numNodes]; diff --git a/src/tests/tribol_mortar_force.cpp b/src/tests/tribol_mortar_force.cpp index 6834c482..0d3b5e43 100644 --- a/src/tests/tribol_mortar_force.cpp +++ b/src/tests/tribol_mortar_force.cpp @@ -7,11 +7,6 @@ #include "tribol/interface/tribol.hpp" #include "tribol/common/Parameters.hpp" #include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/mesh/MeshData.hpp" -#include "tribol/physics/Mortar.hpp" -#include "tribol/physics/AlignedMortar.hpp" -#include "tribol/geom/GeomUtilities.hpp" -#include "tribol/utils/TestUtils.hpp" #ifdef TRIBOL_USE_UMPIRE // Umpire includes @@ -80,25 +75,13 @@ class MortarForceTest : public ::testing::Test { tribol::registerMesh( mortarMeshId, 1, this->numNodes, conn1, cellType, x, y, z, tribol::MemorySpace::Host ); tribol::registerMesh( nonmortarMeshId, 1, this->numNodes, conn2, cellType, x, y, z, tribol::MemorySpace::Host ); - // register nodal forces - RealT *fx1, *fy1, *fz1; - RealT *fx2, *fy2, *fz2; + tribol::Array1D fx1( this->numNodes ); + tribol::Array1D fy1( this->numNodes ); + tribol::Array1D fz1( this->numNodes ); - RealT forceX1[this->numNodes]; - RealT forceY1[this->numNodes]; - RealT forceZ1[this->numNodes]; - - RealT forceX2[this->numNodes]; - RealT forceY2[this->numNodes]; - RealT forceZ2[this->numNodes]; - - fx1 = forceX1; - fy1 = forceY1; - fz1 = forceZ1; - - fx2 = forceX2; - fy2 = forceY2; - fz2 = forceZ2; + tribol::Array1D fx2( this->numNodes ); + tribol::Array1D fy2( this->numNodes ); + tribol::Array1D fz2( this->numNodes ); // initialize force arrays for ( int i = 0; i < this->numNodes; ++i ) { @@ -110,8 +93,8 @@ class MortarForceTest : public ::testing::Test { fz2[i] = 0.; } - tribol::registerNodalResponse( mortarMeshId, fx1, fy1, fz1 ); - tribol::registerNodalResponse( nonmortarMeshId, fx2, fy2, fz2 ); + tribol::registerNodalResponse( mortarMeshId, fx1.data(), fy1.data(), fz1.data() ); + tribol::registerNodalResponse( nonmortarMeshId, fx2.data(), fy2.data(), fz2.data() ); // register nodal pressure and nodal gap array for the nonmortar mesh RealT *gaps, *pressures; @@ -348,15 +331,15 @@ TEST_F( MortarForceTest, parallel_misaligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->checkMortarForces( &conn1[0], &conn2[0], tribol::SINGLE_MORTAR ); + this->checkMortarForces( conn1.data(), conn2.data(), tribol::SINGLE_MORTAR ); } TEST_F( MortarForceTest, parallel_aligned ) @@ -416,15 +399,15 @@ TEST_F( MortarForceTest, parallel_aligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->checkMortarForces( &conn1[0], &conn2[0], tribol::SINGLE_MORTAR ); + this->checkMortarForces( conn1.data(), conn2.data(), tribol::SINGLE_MORTAR ); } TEST_F( MortarForceTest, non_parallel_misaligned ) @@ -484,15 +467,15 @@ TEST_F( MortarForceTest, non_parallel_misaligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->checkMortarForces( &conn1[0], &conn2[0], tribol::SINGLE_MORTAR ); + this->checkMortarForces( conn1.data(), conn2.data(), tribol::SINGLE_MORTAR ); } TEST_F( MortarForceTest, non_parallel_aligned ) @@ -552,15 +535,15 @@ TEST_F( MortarForceTest, non_parallel_aligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->checkMortarForces( &conn1[0], &conn2[0], tribol::SINGLE_MORTAR ); + this->checkMortarForces( conn1.data(), conn2.data(), tribol::SINGLE_MORTAR ); } TEST_F( MortarForceTest, parallel_simple_aligned ) @@ -620,15 +603,15 @@ TEST_F( MortarForceTest, parallel_simple_aligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->checkMortarForces( &conn1[0], &conn2[0], tribol::ALIGNED_MORTAR ); + this->checkMortarForces( conn1.data(), conn2.data(), tribol::ALIGNED_MORTAR ); } int main( int argc, char* argv[] ) diff --git a/src/tests/tribol_mortar_gap.cpp b/src/tests/tribol_mortar_gap.cpp index 244f714f..4b86f55d 100644 --- a/src/tests/tribol_mortar_gap.cpp +++ b/src/tests/tribol_mortar_gap.cpp @@ -13,8 +13,6 @@ #include "tribol/physics/AlignedMortar.hpp" #include "tribol/geom/NodalNormal.hpp" #include "tribol/geom/ElementNormal.hpp" -#include "tribol/geom/GeomUtilities.hpp" -#include "tribol/utils/TestUtils.hpp" #ifdef TRIBOL_USE_UMPIRE // Umpire includes @@ -83,52 +81,28 @@ class MortarGapTest : public ::testing::Test { { // declare arrays to hold stacked coordinates for each // face used in initializing a SurfaceContactElem struct - RealT xyz1[this->dim * this->numNodesPerFace]; - RealT xyz2[this->dim * this->numNodesPerFace]; + tribol::Array2D xyz1( this->numNodesPerFace, this->dim ); + tribol::Array2D xyz2( this->numNodesPerFace, this->dim ); // declare array to hold overlap vertices used for // initializing a SurfaceContactElem struct - RealT xyzOverlap[this->dim * this->numOverlapNodes]; - - // assign pointers to arrays - RealT* xy1 = xyz1; - RealT* xy2 = xyz2; - RealT* xyOverlap = xyzOverlap; - - // grab coordinate data - RealT* x1 = this->x1; - RealT* y1 = this->y1; - RealT* z1 = this->z1; - RealT* x2 = this->x2; - RealT* y2 = this->y2; - RealT* z2 = this->z2; - RealT* xo = this->xOverlap; - RealT* yo = this->yOverlap; - RealT* zo = this->zOverlap; + tribol::Array2D xyzOverlap( this->numOverlapNodes, this->dim ); // generate stacked coordinate array for ( int j = 0; j < this->numNodesPerFace; ++j ) { - for ( int k = 0; k < this->dim; ++k ) { - int id = this->dim * j + k; - // int id2 = this->dim * this->numNodesPerFace + id; - switch ( k ) { - case 0: - xy1[id] = x1[j]; - xy2[id] = x2[j]; - xyOverlap[id] = xo[j]; - break; - case 1: - xy1[id] = y1[j]; - xy2[id] = y2[j]; - xyOverlap[id] = yo[j]; - break; - case 2: - xy1[id] = z1[j]; - xy2[id] = z2[j]; - xyOverlap[id] = zo[j]; - break; - } // end switch - } // end loop over dimension + xyz1( j, 0 ) = x1[j]; + xyz1( j, 1 ) = y1[j]; + xyz1( j, 2 ) = z1[j]; + + xyz2( j, 0 ) = x2[j]; + xyz2( j, 1 ) = y2[j]; + xyz2( j, 2 ) = z2[j]; + } + + for ( int j = 0; j < this->numOverlapNodes; ++j ) { + xyzOverlap( j, 0 ) = xOverlap[j]; + xyzOverlap( j, 1 ) = yOverlap[j]; + xyzOverlap( j, 2 ) = zOverlap[j]; } // end loop over nodes // register the mesh with tribol @@ -180,8 +154,8 @@ class MortarGapTest : public ::testing::Test { // instantiate SurfaceContactElem struct. Note, this object is instantiated // using face 1, face 2, and the set overlap polygon. Note, the mesh ids are set // equal to 0, and the face ids are 0 and 1, respectively. - tribol::SurfaceContactElem elem( this->dim, xy1, xy2, xyOverlap, this->numNodesPerFace, this->numOverlapNodes, - &mortarView, &nonmortarView, 0, 0 ); + tribol::SurfaceContactElem elem( this->dim, xyz1.data(), xyz2.data(), xyzOverlap.data(), this->numNodesPerFace, + this->numOverlapNodes, &mortarView, &nonmortarView, 0, 0 ); // compute the mortar weights to be stored on // the surface contact element struct. @@ -399,8 +373,8 @@ TEST_F( MortarGapTest, parallel_misaligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; @@ -487,8 +461,8 @@ TEST_F( MortarGapTest, parallel_aligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; @@ -575,8 +549,8 @@ TEST_F( MortarGapTest, parallel_simple_aligned ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; diff --git a/src/tests/tribol_mortar_jacobian.cpp b/src/tests/tribol_mortar_jacobian.cpp index a9d4ac0b..e9b32fa9 100644 --- a/src/tests/tribol_mortar_jacobian.cpp +++ b/src/tests/tribol_mortar_jacobian.cpp @@ -6,13 +6,6 @@ // Tribol includes #include "tribol/interface/tribol.hpp" #include "tribol/common/Parameters.hpp" -#include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/mesh/MethodCouplingData.hpp" -#include "tribol/mesh/MeshData.hpp" -#include "tribol/physics/Mortar.hpp" -#include "tribol/physics/AlignedMortar.hpp" -#include "tribol/geom/GeomUtilities.hpp" -#include "tribol/utils/TestUtils.hpp" // Axom includes #include "axom/slic.hpp" @@ -24,10 +17,8 @@ #include "gtest/gtest.h" // c++ includes -#include // std::abs #include #include -#include #include using RealT = tribol::RealT; @@ -82,30 +73,13 @@ class MortarJacTest : public ::testing::Test { tribol::registerMesh( mortarMeshId, 1, this->numNodes, conn1, cellType, x, y, z, tribol::MemorySpace::Host ); tribol::registerMesh( nonmortarMeshId, 1, this->numNodes, conn2, cellType, x, y, z, tribol::MemorySpace::Host ); - // register nodal forces. Note, I was getting a seg fault when - // registering the same pointer to a single set of force arrays - // for both calls to tribol::registerNodalResponse(). As a result, - // I created two sets of nodal force arrays with their own pointers - // to the data that are registered with tribol and there is no longer - // a seg fault. - RealT *fx1, *fy1, *fz1; - RealT *fx2, *fy2, *fz2; + tribol::Array1D fx1( this->numNodes ); + tribol::Array1D fy1( this->numNodes ); + tribol::Array1D fz1( this->numNodes ); - RealT forceX1[this->numNodes]; - RealT forceY1[this->numNodes]; - RealT forceZ1[this->numNodes]; - - RealT forceX2[this->numNodes]; - RealT forceY2[this->numNodes]; - RealT forceZ2[this->numNodes]; - - fx1 = forceX1; - fy1 = forceY1; - fz1 = forceZ1; - - fx2 = forceX2; - fy2 = forceY2; - fz2 = forceZ2; + tribol::Array1D fx2( this->numNodes ); + tribol::Array1D fy2( this->numNodes ); + tribol::Array1D fz2( this->numNodes ); // initialize force arrays for ( int i = 0; i < this->numNodes; ++i ) { @@ -118,8 +92,8 @@ class MortarJacTest : public ::testing::Test { fz2[i] = 0.; } - tribol::registerNodalResponse( mortarMeshId, fx1, fy1, fz1 ); - tribol::registerNodalResponse( nonmortarMeshId, fx2, fy2, fz2 ); + tribol::registerNodalResponse( mortarMeshId, fx1.data(), fy1.data(), fz1.data() ); + tribol::registerNodalResponse( nonmortarMeshId, fx2.data(), fy2.data(), fz2.data() ); gaps = tribol::ArrayT( this->numNodes, this->numNodes ); // length of total mesh to use global connectivity to index @@ -245,15 +219,15 @@ TEST_F( MortarJacTest, jac_input_test ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->setupTribol( &conn1[0], &conn2[0], tribol::ALIGNED_MORTAR ); + this->setupTribol( conn1.data(), conn2.data(), tribol::ALIGNED_MORTAR ); // check Jacobian sparse matrix mfem::SparseMatrix* jac{ nullptr }; @@ -320,15 +294,15 @@ TEST_F( MortarJacTest, update_jac_test ) // register a tribol mesh for computing mortar gaps int numNodesPerFace = 4; - tribol::IndexT conn1[numNodesPerFace]; - tribol::IndexT conn2[numNodesPerFace]; + tribol::Array1D conn1( numNodesPerFace ); + tribol::Array1D conn2( numNodesPerFace ); for ( int i = 0; i < numNodesPerFace; ++i ) { conn1[i] = i; conn2[i] = numNodesPerFace + i; } - this->setupTribol( &conn1[0], &conn2[0], tribol::ALIGNED_MORTAR ); + this->setupTribol( conn1.data(), conn2.data(), tribol::ALIGNED_MORTAR ); // check Jacobian sparse matrix mfem::SparseMatrix* jac{ nullptr }; diff --git a/src/tests/tribol_mortar_lm_patch_test.cpp b/src/tests/tribol_mortar_lm_patch_test.cpp index 0ba9fbf5..7700e56e 100644 --- a/src/tests/tribol_mortar_lm_patch_test.cpp +++ b/src/tests/tribol_mortar_lm_patch_test.cpp @@ -7,12 +7,6 @@ #include "tribol/interface/tribol.hpp" #include "tribol/utils/TestUtils.hpp" #include "tribol/common/Parameters.hpp" -#include "tribol/mesh/MethodCouplingData.hpp" -#include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/mesh/MeshData.hpp" -#include "tribol/physics/Mortar.hpp" -#include "tribol/physics/AlignedMortar.hpp" -#include "tribol/geom/GeomUtilities.hpp" #ifdef TRIBOL_USE_UMPIRE // Umpire includes @@ -126,16 +120,16 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl // setup a temporary stacked array of nodal coordinates for use in the instantiation // of a reference configuration grid function. The nodal coordinates are stacked x, then // y, then z. Also setup a local interleaved array of incremental nodal displacements. - RealT xyz[this->m_mesh.dim * this->m_mesh.numTotalNodes]; - RealT xyz_inc[this->m_mesh.dim * this->m_mesh.numTotalNodes]; + tribol::Array2D xyz( m_mesh.dim, m_mesh.numTotalNodes ); + tribol::Array2D xyz_inc( m_mesh.dim, m_mesh.numTotalNodes ); for ( int i = 0; i < this->m_mesh.numTotalNodes; ++i ) { - xyz[i] = this->m_mesh.x[i]; - xyz[this->m_mesh.numTotalNodes + i] = this->m_mesh.y[i]; - xyz[2 * this->m_mesh.numTotalNodes + i] = this->m_mesh.z[i]; + xyz( 0, i ) = this->m_mesh.x[i]; + xyz( 1, i ) = this->m_mesh.y[i]; + xyz( 2, i ) = this->m_mesh.z[i]; - xyz_inc[this->m_mesh.dim * i] = 0.; - xyz_inc[this->m_mesh.dim * i + 1] = 0.; - xyz_inc[this->m_mesh.dim * i + 2] = 0.; + xyz_inc( 0, i ) = 0.0; + xyz_inc( 1, i ) = 0.0; + xyz_inc( 2, i ) = 0.0; } // define the FE collection and finite element space @@ -177,15 +171,9 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl // (space dimension) x (total number of mesh nodes) + (number of nonmortar nodes in contact), // where the last addition is for the pressure lagrange multiplier field int rhs_size = this->m_mesh.dim * this->m_mesh.numTotalNodes + this->m_mesh.numNonmortarSurfaceNodes; - RealT b[rhs_size]; - - // initialize b vector - for ( int i = 0; i < rhs_size; ++i ) { - b[i] = 0.; - } - // instantiate mfem vector for right hand side - mfem::Vector rhs( &b[0], rhs_size ); + mfem::Vector rhs( rhs_size ); + rhs = 0.0; //////////////////////////////////////////////////// // // @@ -210,7 +198,7 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl // called again with the updated contact solution, so doing so would actually // mess with the RHS in a negative way if ( contact ) { - this->m_mesh.getGapEvals( &b[0] ); + this->m_mesh.getGapEvals( rhs.GetData() ); } // zero out all Dirichlet BC components for each block @@ -293,13 +281,13 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl // mesh array as that of the reference configuration for ( int i = 0; i < this->m_mesh.numTotalNodes; ++i ) { for ( int j = 0; j < this->m_mesh.dim; ++j ) { - xyz[this->m_mesh.numTotalNodes * j + i] += sol_data[this->m_mesh.dim * i + j]; - xyz_inc[this->m_mesh.numTotalNodes * j + i] = sol_data[this->m_mesh.dim * i + j]; + xyz( j, i ) += sol_data[this->m_mesh.dim * i + j]; + xyz_inc( j, i ) = sol_data[this->m_mesh.dim * i + j]; } } // compute stress update - mfem::GridFunction u( fe_space, &xyz_inc[0] ); + mfem::GridFunction u( fe_space, xyz_inc.data() ); const int tdim = this->m_mesh.dim * ( this->m_mesh.dim + 1 ) / 2; mfem::FiniteElementSpace flux_fespace( this->m_mesh.mfem_mesh, &fe_coll, tdim ); mfem::GridFunction stress( &flux_fespace ); @@ -343,7 +331,7 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl this->m_mesh.mfem_mesh->PrintVTK( mesh_ref_ofs, 0, 0 ); // set the current configuration vector and mesh node grid function for output - mfem::Vector x_cur( &xyz[0], this->m_mesh.dim * this->m_mesh.numTotalNodes ); + mfem::Vector x_cur( xyz.data(), this->m_mesh.dim * this->m_mesh.numTotalNodes ); this->m_mesh.mfem_mesh->SetNodes( x_cur ); // print current mesh @@ -381,7 +369,7 @@ void MortarLMPatchTest::computeContactSolution( int nMortarElemsX, int nMortarEl rhs_vec.open( "rhs_" + suffix_rhs.str() ); for ( int i = 0; i < jac.NumRows(); ++i ) { - rhs_vec << b[i] << "\n"; + rhs_vec << rhs[i] << "\n"; sol_vec << sol_data[i] << "\n"; for ( int j = 0; j < jac.NumCols(); ++j ) { RealT val = dJac( i, j ); diff --git a/src/tests/tribol_mortar_sparse_weights.cpp b/src/tests/tribol_mortar_sparse_weights.cpp index ffef07cc..039e77e8 100644 --- a/src/tests/tribol_mortar_sparse_weights.cpp +++ b/src/tests/tribol_mortar_sparse_weights.cpp @@ -13,9 +13,6 @@ #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/mesh/CouplingScheme.hpp" #include "tribol/mesh/MeshData.hpp" -#include "tribol/physics/Mortar.hpp" -#include "tribol/physics/AlignedMortar.hpp" -#include "tribol/geom/GeomUtilities.hpp" // Axom includes #include "axom/slic.hpp" @@ -31,13 +28,6 @@ // gtest includes #include "gtest/gtest.h" -// c++ includes -#include // std::abs, std::cos, std::sin -#include -#include -#include -#include - using RealT = tribol::RealT; void computeGapsFromSparseWts( tribol::CouplingScheme* cs, RealT* gaps ) @@ -80,7 +70,7 @@ void computeGapsFromSparseWts( tribol::CouplingScheme* cs, RealT* gaps ) int numTotalNodes = static_cast( cs->getMethodData() )->m_numTotalNodes; for ( int a = 0; a < numTotalNodes; ++a ) { // get nonmortar nodal normal - RealT nrml_a[dim]; + tribol::Array1D nrml_a( dim ); nrml_a[0] = nonmortarMesh.getNodalNormals()[0][a]; // array is global length; no index out (?) nrml_a[1] = nonmortarMesh.getNodalNormals()[1][a]; if ( dim == 3 ) { @@ -90,8 +80,8 @@ void computeGapsFromSparseWts( tribol::CouplingScheme* cs, RealT* gaps ) // loop over range of nonzero column entries for ( int b = I[a]; b < I[a + 1]; ++b ) { // get face coordinates for node J[b], i.e. column node id - RealT mortar_xyz[dim]; - RealT nonmortar_xyz[dim]; + tribol::Array1D mortar_xyz( dim ); + tribol::Array1D nonmortar_xyz( dim ); for ( int i = 0; i < dim; ++i ) { mortar_xyz[i] = 0.; diff --git a/src/tests/tribol_mortar_wts.cpp b/src/tests/tribol_mortar_wts.cpp index 474456cd..d56ac766 100644 --- a/src/tests/tribol_mortar_wts.cpp +++ b/src/tests/tribol_mortar_wts.cpp @@ -4,10 +4,8 @@ // SPDX-License-Identifier: (MIT) // Tribol includes -#include "tribol/mesh/MeshData.hpp" #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/physics/Mortar.hpp" -#include "tribol/geom/GeomUtilities.hpp" // Axom includes #include "axom/slic.hpp" @@ -72,51 +70,32 @@ class MortarWeightTest : public ::testing::Test { SLIC_ERROR( "checkMortarWts: number of nodes per face not equal to 4." ); } - RealT xyz1[this->dim * this->numNodesPerFace]; - RealT xyz2[this->dim * this->numNodesPerFace]; - RealT xyzOverlap[this->dim * this->numOverlapNodes]; - RealT* xy1 = xyz1; - RealT* xy2 = xyz2; - RealT* xyOverlap = xyzOverlap; - - RealT* x1 = this->x1; - RealT* y1 = this->y1; - RealT* z1 = this->z1; - RealT* x2 = this->x2; - RealT* y2 = this->y2; - RealT* z2 = this->z2; - RealT* xo = this->xOverlap; - RealT* yo = this->yOverlap; - RealT* zo = this->zOverlap; + tribol::Array2D xyz1( this->numNodesPerFace, this->dim ); + tribol::Array2D xyz2( this->numNodesPerFace, this->dim ); + tribol::Array2D xyzOverlap( this->numOverlapNodes, this->dim ); // generate stacked coordinate array for ( int j = 0; j < this->numNodesPerFace; ++j ) { - for ( int k = 0; k < this->dim; ++k ) { - switch ( k ) { - case 0: - xy1[this->dim * j + k] = x1[j]; - xy2[this->dim * j + k] = x2[j]; - xyOverlap[this->dim * j + k] = xo[j]; - break; - case 1: - xy1[this->dim * j + k] = y1[j]; - xy2[this->dim * j + k] = y2[j]; - xyOverlap[this->dim * j + k] = yo[j]; - break; - case 2: - xy1[this->dim * j + k] = z1[j]; - xy2[this->dim * j + k] = z2[j]; - xyOverlap[this->dim * j + k] = zo[j]; - break; - } // end switch - } // end loop over dimension + xyz1( j, 0 ) = x1[j]; + xyz1( j, 1 ) = y1[j]; + xyz1( j, 2 ) = z1[j]; + + xyz2( j, 0 ) = x2[j]; + xyz2( j, 1 ) = y2[j]; + xyz2( j, 2 ) = z2[j]; + } + + for ( int j = 0; j < this->numOverlapNodes; ++j ) { + xyzOverlap( j, 0 ) = xOverlap[j]; + xyzOverlap( j, 1 ) = yOverlap[j]; + xyzOverlap( j, 2 ) = zOverlap[j]; } // end loop over nodes // instantiate SurfaceContactElem struct. Note, this object is instantiated // using face 1, face 2, and the set overlap polygon. Note, the mesh ids are set // equal to 0, and the face ids are 0 and 1, respectively. - tribol::SurfaceContactElem elem( this->dim, xy1, xy2, xyOverlap, this->numNodesPerFace, this->numOverlapNodes, 0, 0, - 0, 1 ); + tribol::SurfaceContactElem elem( this->dim, xyz1.data(), xyz2.data(), xyzOverlap.data(), this->numNodesPerFace, + this->numOverlapNodes, 0, 0, 0, 1 ); // compute the mortar weights to be stored on // the surface contact element struct. diff --git a/src/tests/tribol_nodal_nrmls.cpp b/src/tests/tribol_nodal_nrmls.cpp index 4b848273..570f236e 100644 --- a/src/tests/tribol_nodal_nrmls.cpp +++ b/src/tests/tribol_nodal_nrmls.cpp @@ -7,7 +7,6 @@ #include "tribol/interface/tribol.hpp" #include "tribol/mesh/MeshData.hpp" #include "tribol/geom/ElementNormal.hpp" -#include "tribol/geom/GeomUtilities.hpp" #include "tribol/geom/NodalNormal.hpp" #include "tribol/utils/Math.hpp" @@ -16,9 +15,6 @@ #include "umpire/ResourceManager.hpp" #endif -// Axom includes -#include "axom/slic.hpp" - // gtest includes #include "gtest/gtest.h" @@ -37,7 +33,7 @@ class NodalNormalTest : public ::testing::Test { int dim; void computeNodalNormals( int cell_type, RealT const* const x, RealT const* const y, RealT const* const z, - const tribol::IndexT* conn, int const numCells, int const numNodes, int const dim ) + const tribol::IndexT* conn, int const numCells, int const numNodes ) { // register the mesh with tribol const tribol::IndexT mesh_id = 0; @@ -76,23 +72,22 @@ TEST_F( NodalNormalTest, two_quad_inverted_v ) int numFaces = 2; int numNodesPerFace = 4; int cellType = (int)( tribol::LINEAR_QUAD ); - tribol::IndexT conn[numFaces * numNodesPerFace]; + tribol::Array2D conn( numFaces, numNodesPerFace ); // setup connectivity for the two faces ensuring nodes // are ordered consistent with an outward unit normal for ( int i = 0; i < numNodesPerFace; ++i ) { - conn[i] = i; + conn( 0, i ) = i; } - - conn[4] = 0; - conn[5] = 5; - conn[6] = 4; - conn[7] = 1; + conn( 1, 0 ) = 0; + conn( 1, 1 ) = 5; + conn( 1, 2 ) = 4; + conn( 1, 3 ) = 1; // setup the nodal coordinates of the mesh - RealT x[numNodesPerFace + 2]; - RealT y[numNodesPerFace + 2]; - RealT z[numNodesPerFace + 2]; + tribol::ArrayT x( numNodesPerFace + 2 ); + tribol::ArrayT y( numNodesPerFace + 2 ); + tribol::ArrayT z( numNodesPerFace + 2 ); x[0] = 0.; x[1] = 0.; @@ -116,7 +111,7 @@ TEST_F( NodalNormalTest, two_quad_inverted_v ) z[5] = -1.; // compute the nodal normals - computeNodalNormals( cellType, &x[0], &y[0], &z[0], &conn[0], numFaces, 6, 3 ); + computeNodalNormals( cellType, x.data(), y.data(), z.data(), conn.data(), numFaces, 6 ); tribol::MeshManager& meshManager = tribol::MeshManager::getInstance(); auto mesh = meshManager.at( 0 ).getView(); diff --git a/src/tests/tribol_quad_integ.cpp b/src/tests/tribol_quad_integ.cpp index ff080611..b63b79ff 100644 --- a/src/tests/tribol_quad_integ.cpp +++ b/src/tests/tribol_quad_integ.cpp @@ -4,15 +4,11 @@ // SPDX-License-Identifier: (MIT) // Tribol includes -#include "tribol/mesh/MeshData.hpp" #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/integ/Integration.hpp" #include "tribol/geom/GeomUtilities.hpp" #include "tribol/integ/FE.hpp" -// Axom includes -#include "axom/slic.hpp" - // gtest includes #include "gtest/gtest.h" @@ -39,34 +35,20 @@ class IsoIntegTest : public ::testing::Test { bool integrate( RealT const tol ) { - RealT xyz[this->dim * this->numNodes]; - RealT* xy = xyz; - - RealT* x = this->x; - RealT* y = this->y; - RealT* z = this->z; + tribol::Array2D xyz( this->numNodes, this->dim ); // generate stacked coordinate array for ( int j = 0; j < this->numNodes; ++j ) { - for ( int k = 0; k < this->dim; ++k ) { - switch ( k ) { - case 0: - xy[this->dim * j + k] = x[j]; - break; - case 1: - xy[this->dim * j + k] = y[j]; - break; - case 2: - xy[this->dim * j + k] = z[j]; - break; - } // end switch - } // end loop over dimension + xyz( j, 0 ) = x[j]; + xyz( j, 1 ) = y[j]; + xyz( j, 2 ) = z[j]; } // end loop over nodes // instantiate SurfaceContactElem struct. Note, this object is instantiated // using face 1 as face 2, but these faces are not used in this test so this // is ok. - tribol::SurfaceContactElem elem( this->dim, xy, xy, xy, this->numNodes, this->numNodes, nullptr, nullptr, 0, 0 ); + tribol::SurfaceContactElem elem( this->dim, xyz.data(), xyz.data(), xyz.data(), this->numNodes, this->numNodes, + nullptr, nullptr, 0, 0 ); // instantiate integration object tribol::IntegPts integ; diff --git a/src/tests/tribol_twb_integ.cpp b/src/tests/tribol_twb_integ.cpp index 853c20ca..5b180021 100644 --- a/src/tests/tribol_twb_integ.cpp +++ b/src/tests/tribol_twb_integ.cpp @@ -4,15 +4,11 @@ // SPDX-License-Identifier: (MIT) // Tribol includes -#include "tribol/mesh/MeshData.hpp" #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/integ/Integration.hpp" #include "tribol/geom/GeomUtilities.hpp" #include "tribol/integ/FE.hpp" -// Axom includes -#include "axom/slic.hpp" - // gtest includes #include "gtest/gtest.h" @@ -40,34 +36,20 @@ class TWBIntegTest : public ::testing::Test { bool integrate( RealT const tol ) { - RealT xyz[this->dim * this->numNodes]; - RealT* xy = xyz; - - RealT* x = this->x; - RealT* y = this->y; - RealT* z = this->z; + tribol::Array2D xyz( this->numNodes, this->dim ); // generate stacked coordinate array for ( int j = 0; j < this->numNodes; ++j ) { - for ( int k = 0; k < this->dim; ++k ) { - switch ( k ) { - case 0: - xy[this->dim * j + k] = x[j]; - break; - case 1: - xy[this->dim * j + k] = y[j]; - break; - case 2: - xy[this->dim * j + k] = z[j]; - break; - } // end switch - } // end loop over dimension + xyz( j, 0 ) = x[j]; + xyz( j, 1 ) = y[j]; + xyz( j, 2 ) = z[j]; } // end loop over nodes // instantiate SurfaceContactElem struct. Note, this object is instantiated // using face 1 as face 2, but these faces are not used in this test so this // is ok. - tribol::SurfaceContactElem elem( this->dim, xy, xy, xy, this->numNodes, this->numNodes, nullptr, nullptr, 0, 0 ); + tribol::SurfaceContactElem elem( this->dim, xyz.data(), xyz.data(), xyz.data(), this->numNodes, this->numNodes, + nullptr, nullptr, 0, 0 ); // instantiate integration object tribol::IntegPts integ; @@ -82,8 +64,8 @@ class TWBIntegTest : public ::testing::Test { for ( int a = 0; a < this->numNodes; ++a ) { for ( int ip = 0; ip < integ.numIPs; ++ip ) { - tribol::WachspressBasis( xy, integ.xy[dim * ip], integ.xy[dim * ip + 1], integ.xy[dim * ip + 2], this->numNodes, - a, phi ); + tribol::WachspressBasis( xyz.data(), integ.xy[dim * ip], integ.xy[dim * ip + 1], integ.xy[dim * ip + 2], + this->numNodes, a, phi ); areaTest += integ.wts[ip] * phi; } diff --git a/src/tribol/CMakeLists.txt b/src/tribol/CMakeLists.txt index 2d26f6f6..01caae75 100644 --- a/src/tribol/CMakeLists.txt +++ b/src/tribol/CMakeLists.txt @@ -14,6 +14,7 @@ set(tribol_headers common/ArrayTypes.hpp common/BasicTypes.hpp common/Containers.hpp + common/Enzyme.hpp common/ExecModel.hpp common/LoopExec.hpp common/Parameters.hpp @@ -21,6 +22,7 @@ set(tribol_headers geom/CompGeom.hpp geom/ElementNormal.hpp geom/GeomUtilities.hpp + geom/Hyperplane.hpp geom/NodalNormal.hpp integ/FE.hpp @@ -48,7 +50,7 @@ set(tribol_headers utils/DataManager.hpp utils/Math.hpp utils/TestUtils.hpp - ) +) ## list of sources set(tribol_sources diff --git a/src/tribol/common/ArrayTypes.hpp b/src/tribol/common/ArrayTypes.hpp index 5c50092f..b36c9549 100644 --- a/src/tribol/common/ArrayTypes.hpp +++ b/src/tribol/common/ArrayTypes.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_COMMON_ARRAYTYPES_HPP_ -#define TRIBOL_COMMON_ARRAYTYPES_HPP_ +#ifndef SRC_TRIBOL_COMMON_ARRAYTYPES_HPP_ +#define SRC_TRIBOL_COMMON_ARRAYTYPES_HPP_ // Tribol includes #include "tribol/common/ExecModel.hpp" @@ -70,4 +70,4 @@ using MultiViewArrayView = ArrayViewT, 1, SPACE>; } // namespace tribol -#endif /* TRIBOL_COMMON_ARRAYTYPES_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_ARRAYTYPES_HPP_ */ diff --git a/src/tribol/common/BasicTypes.hpp b/src/tribol/common/BasicTypes.hpp index 41cc1a02..49617b31 100644 --- a/src/tribol/common/BasicTypes.hpp +++ b/src/tribol/common/BasicTypes.hpp @@ -3,20 +3,23 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_COMMON_BASICTYPES_HPP_ -#define TRIBOL_COMMON_BASICTYPES_HPP_ +#ifndef SRC_TRIBOL_COMMON_BASICTYPES_HPP_ +#define SRC_TRIBOL_COMMON_BASICTYPES_HPP_ -// Tribol includes +// Tribol config include #include "tribol/config.hpp" -// Axom includes -#include "axom/core/Types.hpp" +// C includes +#include // MPI includes #ifdef TRIBOL_USE_MPI #include #endif +// Axom includes +#include "axom/core/Types.hpp" + namespace tribol { #ifdef TRIBOL_USE_MPI @@ -36,6 +39,9 @@ using CommT = int; // match index type used in axom (since data is held in axom data structures) using IndexT = axom::IndexType; +// size type matching size of addressable memory +using SizeT = size_t; + #ifdef TRIBOL_USE_SINGLE_PRECISION #error "Tribol does not support single precision." @@ -71,11 +77,23 @@ using RealT = double; #define TRIBOL_DEFAULT_HOST_DEVICE #endif -// Define variable when loops are computed on host +// Defined when Tribol doesn't have a device available #if !( defined( TRIBOL_USE_CUDA ) || defined( TRIBOL_USE_HIP ) ) #define TRIBOL_USE_HOST #endif +// Define variable when in device code +#if defined( __CUDA_ARCH__ ) || defined( __HIP_DEVICE_COMPILE__ ) +#define TRIBOL_DEVICE_CODE +#endif + +// Ignore host code in __host__ __device__ code warning on NVCC +#ifdef TRIBOL_USE_CUDA +#define TRIBOL_NVCC_EXEC_CHECK_DISABLE #pragma nv_exec_check_disable +#else +#define TRIBOL_NVCC_EXEC_CHECK_DISABLE +#endif + } // namespace tribol -#endif /* TRIBOL_COMMON_BASICTYPES_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_BASICTYPES_HPP_ */ diff --git a/src/tribol/common/Containers.hpp b/src/tribol/common/Containers.hpp index 9ef8743a..e0af5f98 100644 --- a/src/tribol/common/Containers.hpp +++ b/src/tribol/common/Containers.hpp @@ -3,8 +3,11 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_COMMON_CONTAINERS_HPP_ -#define TRIBOL_COMMON_CONTAINERS_HPP_ +#ifndef SRC_TRIBOL_COMMON_CONTAINERS_HPP_ +#define SRC_TRIBOL_COMMON_CONTAINERS_HPP_ + +// C++ includes +#include // Tribol includes #include "tribol/common/BasicTypes.hpp" @@ -218,4 +221,4 @@ class StackArray { }; } // namespace tribol -#endif /* TRIBOL_COMMON_CONTAINERS_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_CONTAINERS_HPP_ */ diff --git a/src/tribol/common/Enzyme.hpp b/src/tribol/common/Enzyme.hpp new file mode 100644 index 00000000..a384d9f5 --- /dev/null +++ b/src/tribol/common/Enzyme.hpp @@ -0,0 +1,55 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Tribol Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (MIT) + +#ifndef SRC_TRIBOL_COMMON_ENZYME_HPP_ +#define SRC_TRIBOL_COMMON_ENZYME_HPP_ + +#include "tribol/config.hpp" + +#include "tribol/common/BasicTypes.hpp" + +#ifdef TRIBOL_USE_ENZYME +#include "mfem/general/enzyme.hpp" + +#if !defined( TRIBOL_USE_HOST ) && !defined( TRIBOL_DEVICE_CODE ) +// When compiling with NVCC or HIPCC, the compiler performs multiple passes. +// In the HOST pass, reading the __device__ or __constant__ globals (like +// enzyme_const) triggers warnings. +// +// To fix this, we declare host-side aliases that map to the same +// underlying assembly symbols required by the Enzyme pass, but without +// the device attributes. +extern "C" { +extern int tribol_host_enzyme_const asm( "enzyme_const" ); +extern int tribol_host_enzyme_dup asm( "enzyme_dup" ); +extern int tribol_host_enzyme_out asm( "enzyme_out" ); +extern int tribol_host_enzyme_dupnoneed asm( "enzyme_dupnoneed" ); +} + +#define TRIBOL_ENZYME_CONST tribol_host_enzyme_const +#define TRIBOL_ENZYME_DUP tribol_host_enzyme_dup +#define TRIBOL_ENZYME_OUT tribol_host_enzyme_out +#define TRIBOL_ENZYME_DUPNONEED tribol_host_enzyme_dupnoneed + +#else +// We are either: +// 1. In a device compilation pass (__CUDA_ARCH__ or __HIP_DEVICE_COMPILE__ defined). +// 2. Using a standard host compiler (GCC, Clang, etc.). +// In these cases, we use the symbols provided by the Enzyme/MFEM headers. +#define TRIBOL_ENZYME_CONST enzyme_const +#define TRIBOL_ENZYME_DUP enzyme_dup +#define TRIBOL_ENZYME_OUT enzyme_out +#define TRIBOL_ENZYME_DUPNONEED enzyme_dupnoneed +#endif + +#else +// Fallback definitions if Enzyme is disabled +#define TRIBOL_ENZYME_CONST 0 +#define TRIBOL_ENZYME_DUP 0 +#define TRIBOL_ENZYME_OUT 0 +#define TRIBOL_ENZYME_DUPNONEED 0 +#endif + +#endif /* SRC_TRIBOL_COMMON_ENZYME_HPP_ */ \ No newline at end of file diff --git a/src/tribol/common/ExecModel.hpp b/src/tribol/common/ExecModel.hpp index ae644333..ef4ebed2 100644 --- a/src/tribol/common/ExecModel.hpp +++ b/src/tribol/common/ExecModel.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_COMMON_EXECMODEL_HPP_ -#define SRC_COMMON_EXECMODEL_HPP_ +#ifndef SRC_TRIBOL_COMMON_EXECMODEL_HPP_ +#define SRC_TRIBOL_COMMON_EXECMODEL_HPP_ // Tribol includes #include "tribol/common/BasicTypes.hpp" @@ -102,6 +102,8 @@ inline umpire::resource::MemoryResourceType toUmpireMemoryType( MemorySpace mem_ #endif +inline int getDefaultAllocatorID() { return axom::getDefaultAllocatorID(); } + inline int getResourceAllocatorID( MemorySpace mem_space ) { int allocator_id = axom::getDefaultAllocatorID(); @@ -142,4 +144,4 @@ inline bool isOnDevice( ExecutionMode exec ) } // namespace tribol -#endif /* SRC_COMMON_EXECMODEL_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_EXECMODEL_HPP_ */ diff --git a/src/tribol/common/LoopExec.hpp b/src/tribol/common/LoopExec.hpp index ad692a33..5edf3bdf 100644 --- a/src/tribol/common/LoopExec.hpp +++ b/src/tribol/common/LoopExec.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_COMMON_LOOPEXEC_HPP_ -#define SRC_COMMON_LOOPEXEC_HPP_ +#ifndef SRC_TRIBOL_COMMON_LOOPEXEC_HPP_ +#define SRC_TRIBOL_COMMON_LOOPEXEC_HPP_ // C++ includes #include @@ -185,4 +185,4 @@ void forAllExec( ExecutionMode exec_mode, IndexT N, BODY&& body ) } // namespace tribol -#endif /* SRC_COMMON_LOOPEXEC_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_LOOPEXEC_HPP_ */ diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 364858ff..c944a135 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_PARAMETERS_HPP_ -#define TRIBOL_PARAMETERS_HPP_ +#ifndef SRC_TRIBOL_COMMON_PARAMETERS_HPP_ +#define SRC_TRIBOL_COMMON_PARAMETERS_HPP_ // Tribol includes #include "tribol/common/BasicTypes.hpp" @@ -500,4 +500,4 @@ struct Parameters { } // namespace tribol -#endif /* TRIBOL_PARAMETERS_HPP_ */ +#endif /* SRC_TRIBOL_COMMON_PARAMETERS_HPP_ */ diff --git a/src/tribol/config.hpp.in b/src/tribol/config.hpp.in index 52df83a1..8b5c36e8 100644 --- a/src/tribol/config.hpp.in +++ b/src/tribol/config.hpp.in @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_CONFIG_HPP_ -#define TRIBOL_CONFIG_HPP_ +#ifndef SRC_TRIBOL_CONFIG_HPP_ +#define SRC_TRIBOL_CONFIG_HPP_ /* * Note: Use only C-style comments in this file since it might be included from a C file @@ -40,4 +40,4 @@ #cmakedefine TRIBOL_USE_OPENMP #cmakedefine BUILD_REDECOMP -#endif /* TRIBOL_CONFIG_HPP_ */ +#endif /* SRC_TRIBOL_CONFIG_HPP_ */ diff --git a/src/tribol/geom/CompGeom.hpp b/src/tribol/geom/CompGeom.hpp index e90a485a..9ee22e5f 100644 --- a/src/tribol/geom/CompGeom.hpp +++ b/src/tribol/geom/CompGeom.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_GEOM_COMPGEOM_HPP_ -#define SRC_GEOM_COMPGEOM_HPP_ +#ifndef SRC_TRIBOL_GEOM_COMPGEOM_HPP_ +#define SRC_TRIBOL_GEOM_COMPGEOM_HPP_ #include "tribol/mesh/MeshData.hpp" #include "tribol/mesh/InterfacePairs.hpp" @@ -950,4 +950,4 @@ TRIBOL_HOST_DEVICE FaceGeomException CheckInterfacePair( InterfacePair& pair, co } // namespace tribol -#endif /* SRC_GEOM_COMPGEOM_HPP_ */ +#endif /* SRC_TRIBOL_GEOM_COMPGEOM_HPP_ */ diff --git a/src/tribol/geom/ElementNormal.hpp b/src/tribol/geom/ElementNormal.hpp index 3cabfc43..0e9c340d 100644 --- a/src/tribol/geom/ElementNormal.hpp +++ b/src/tribol/geom/ElementNormal.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_GEOM_ELEMENTNORMAL_HPP_ -#define SRC_GEOM_ELEMENTNORMAL_HPP_ +#ifndef SRC_TRIBOL_GEOM_ELEMENTNORMAL_HPP_ +#define SRC_TRIBOL_GEOM_ELEMENTNORMAL_HPP_ #include "tribol/common/BasicTypes.hpp" @@ -76,4 +76,4 @@ class ElementCentroidNormal : public ElementNormal { } // namespace tribol -#endif /* SRC_GEOM_ELEMENTNORMAL_HPP_ */ +#endif /* SRC_TRIBOL_GEOM_ELEMENTNORMAL_HPP_ */ diff --git a/src/tribol/geom/GeomUtilities.cpp b/src/tribol/geom/GeomUtilities.cpp index 5a80add9..4dea0e7e 100644 --- a/src/tribol/geom/GeomUtilities.cpp +++ b/src/tribol/geom/GeomUtilities.cpp @@ -17,7 +17,6 @@ #include #include -#include namespace tribol { diff --git a/src/tribol/geom/GeomUtilities.hpp b/src/tribol/geom/GeomUtilities.hpp index 7f8d08d8..b03af41b 100644 --- a/src/tribol/geom/GeomUtilities.hpp +++ b/src/tribol/geom/GeomUtilities.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_GEOM_GEOMUTILITIES_HPP_ -#define SRC_GEOM_GEOMUTILITIES_HPP_ +#ifndef SRC_TRIBOL_GEOM_GEOMUTILITIES_HPP_ +#define SRC_TRIBOL_GEOM_GEOMUTILITIES_HPP_ #include "tribol/common/Parameters.hpp" #include "tribol/mesh/MeshData.hpp" @@ -1824,4 +1824,4 @@ TRIBOL_HOST_DEVICE bool IsConvex( const RealT* const x, const RealT* const y, co } // namespace tribol -#endif /* SRC_GEOM_GEOMUTILITIES_HPP_ */ +#endif /* SRC_TRIBOL_GEOM_GEOMUTILITIES_HPP_ */ diff --git a/src/tribol/geom/Hyperplane.hpp b/src/tribol/geom/Hyperplane.hpp new file mode 100644 index 00000000..1f3e016d --- /dev/null +++ b/src/tribol/geom/Hyperplane.hpp @@ -0,0 +1,101 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Tribol Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (MIT) + +#ifndef SRC_TRIBOL_GEOM_HYPERPLANE_HPP_ +#define SRC_TRIBOL_GEOM_HYPERPLANE_HPP_ + +// Tribol config include +#include "tribol/config.hpp" + +// Axom includes +#include "axom/primal.hpp" + +// Tribol includes +#include "tribol/common/BasicTypes.hpp" +#include "tribol/geom/Vector.hpp" + +namespace tribol { + +/** + * @brief Represents a hyperplane (affine subspace) in N-dimensional space. + * + * The `Hyperplane` type models a geometric hyperplane defined by an origin point and a normal vector. It provides + * basic operations such as projecting a point onto the hyperplane. The template parameter `_VectorT` allows using + * different vector storage/backing types. + * + * @tparam _T Element numeric type for vector components. + * @tparam _VectorT Vector type used for points and normals (defaults to `Vector<_T>`). + */ +template > +class Hyperplane { + public: + /// @brief Alias for the vector type used by this hyperplane. + using VectorT_ = _VectorT; + + /// @brief Scalar value type for components (derived from the vector type). + using ValueT_ = typename VectorT_::ValueT_; + + static_assert( std::is_same<_T, ValueT_>::value, "Hyperplane must be used with the same type as the vector type" ); + + /** + * @brief Construct a hyperplane from an origin point and a normal vector. + * + * @param origin A point on the hyperplane. + * @param normal A normal vector to the hyperplane (must be non-zero and same dimension as origin). + */ + TRIBOL_HOST_DEVICE Hyperplane( const VectorT_& origin, const VectorT_& normal ) : origin_( origin ), normal_( normal ) + { + assert( normal_.norm() > 0.0 ); + assert( origin_.dim() == normal_.dim() ); + } + /** + * @brief Move-construct a hyperplane from rvalue origin and normal vectors. + * + * @param origin Rvalue reference to a point on the hyperplane (moved into the object). + * @param normal Rvalue reference to the normal vector (moved into the object). + */ + TRIBOL_HOST_DEVICE Hyperplane( VectorT_&& origin, VectorT_&& normal ) + : origin_( std::move( origin ) ), normal_( std::move( normal ) ) + { + assert( normal_.norm() > 0.0 ); + assert( origin_.dim() == normal_.dim() ); + } + + /** + * @brief Return the dimension of the hyperplane (same as underlying vectors). + * @return The dimension (number of components) of the origin/normal vectors. + */ + TRIBOL_HOST_DEVICE constexpr SizeT dim() const { return origin_.dim(); } + + /** + * @brief Project a point orthogonally onto the hyperplane. + * + * The function computes the orthogonal projection of `point` onto the hyperplane defined by `origin_` and + * `normal_`. + * + * @param point The point to project (must have same dimension as the hyperplane). + * @return The projected point on the hyperplane. + */ + TRIBOL_HOST_DEVICE VectorT_ projectPoint( const VectorT_& point ) const + { + assert( point.dim() == dim() ); + + // Calculate the projection of the point onto the hyperplane + VectorT_ diff = point - origin_; + ValueT_ dist = diff.dot( normal_ ) / normal_.norm(); + return point - dist * normal_; + } + + private: + /// @brief A point that lies on the hyperplane. + VectorT_ origin_; + + /// @brief The normal vector to the hyperplane (should be non-zero). + VectorT_ normal_; +}; + +} // namespace tribol + +#endif /* SRC_TRIBOL_GEOM_HYPERPLANE_HPP_ */ diff --git a/src/tribol/geom/NodalNormal.cpp b/src/tribol/geom/NodalNormal.cpp index ded86d88..1a96ebf2 100644 --- a/src/tribol/geom/NodalNormal.cpp +++ b/src/tribol/geom/NodalNormal.cpp @@ -5,13 +5,10 @@ #include "NodalNormal.hpp" +#include "tribol/common/Enzyme.hpp" #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/utils/Math.hpp" -#ifdef TRIBOL_USE_ENZYME -#include "mfem/general/enzyme.hpp" -#endif - namespace tribol { // forward declare free functions for enzyme. these shouldn't be used outside the class, so no need to put them in the @@ -103,7 +100,7 @@ void EdgeAvgNodalNormal::Compute( MeshData& mesh, MethodData* jacobian_data ) mesh.allocateNodalNormals(); - auto n0 = Array2D( { 3, mesh.numberOfNodes() }, mesh.getAllocatorId() ); + auto n0 = ArrayT( { 3, mesh.numberOfNodes() }, mesh.getAllocatorId() ); n0.fill( 0.0 ); if ( jacobian_data != nullptr ) { @@ -222,8 +219,9 @@ void ElementEdgeAvgNodalNormalJacobian( [[maybe_unused]] const RealT* x, [[maybe RealT x_dot[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; for ( int i{ 0 }; i < num_nodes_per_elem * 3; ++i ) { x_dot[i] = 1.0; - __enzyme_fwddiff( (void*)ElementEdgeAvgNodalNormal, enzyme_dup, x, x_dot, enzyme_const, xref, enzyme_dup, n, - &dndx[num_nodes_per_elem * 3 * i], enzyme_const, num_nodes_per_elem ); + __enzyme_fwddiff( (void*)ElementEdgeAvgNodalNormal, TRIBOL_ENZYME_DUP, x, x_dot, TRIBOL_ENZYME_CONST, xref, + TRIBOL_ENZYME_DUP, n, &dndx[num_nodes_per_elem * 3 * i], TRIBOL_ENZYME_CONST, + num_nodes_per_elem ); x_dot[i] = 0.0; } #else diff --git a/src/tribol/geom/NodalNormal.hpp b/src/tribol/geom/NodalNormal.hpp index 6fb31435..00214f8c 100644 --- a/src/tribol/geom/NodalNormal.hpp +++ b/src/tribol/geom/NodalNormal.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_GEOM_NODALNORMAL_HPP_ -#define SRC_GEOM_NODALNORMAL_HPP_ +#ifndef SRC_TRIBOL_GEOM_NODALNORMAL_HPP_ +#define SRC_TRIBOL_GEOM_NODALNORMAL_HPP_ #include "tribol/config.hpp" @@ -65,4 +65,4 @@ class EdgeAvgNodalNormal : public NodalNormal { } // namespace tribol -#endif /* SRC_GEOM_NODALNORMAL_HPP_ */ +#endif /* SRC_TRIBOL_GEOM_NODALNORMAL_HPP_ */ diff --git a/src/tribol/integ/FE.hpp b/src/tribol/integ/FE.hpp index 44ef19f5..54e8f76d 100644 --- a/src/tribol/integ/FE.hpp +++ b/src/tribol/integ/FE.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_INTEG_FE_HPP_ -#define SRC_INTEG_FE_HPP_ +#ifndef SRC_TRIBOL_INTEG_FE_HPP_ +#define SRC_TRIBOL_INTEG_FE_HPP_ #include @@ -210,13 +210,17 @@ inline void InvIso( const RealT x[3], const RealT* xA, const RealT* yA, const Re xi[0] = x_sol[0]; xi[1] = x_sol[1]; - // check to make sure point is inside isoparametric quad + // check to make sure point is inside isoparametric quad_wt +#if !defined( TRIBOL_USE_ENZYME ) bool in_quad = true; +#endif if ( std::abs( xi[0] ) > 1. || std::abs( xi[1] ) > 1. ) { if ( std::abs( xi[0] ) > 1. + 100. * xtol || std::abs( xi[1] ) > 1. + 100. * xtol ) // should have some tolerance dependent conv tol? { +#if !defined( TRIBOL_USE_ENZYME ) in_quad = false; +#endif } else { xi[0] = std::min( xi[0], 1. ); xi[1] = std::min( xi[1], 1. ); @@ -412,4 +416,4 @@ void DetJQuad( const RealT xi, const RealT eta, const RealT* x, const int dim, R } // namespace tribol -#endif /* SRC_INTEG_FE_HPP_ */ +#endif /* SRC_TRIBOL_INTEG_FE_HPP_ */ diff --git a/src/tribol/integ/Integration.cpp b/src/tribol/integ/Integration.cpp index 00c837dd..beaf7d28 100644 --- a/src/tribol/integ/Integration.cpp +++ b/src/tribol/integ/Integration.cpp @@ -121,9 +121,8 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ) integ.initialize( 3, numTotalPoints ); - // declare local array to hold barycentric coordinates for each - // triangle - RealT bary[elem.dim * numTriPoints]; + // declare local array to hold barycentric coordinates for each triangle + Array2D bary( numTriPoints, elem.dim ); switch ( k ) { case 2: @@ -133,17 +132,17 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ) } // first barycentric point - bary[0] = 0.1666666666667; - bary[1] = 0.6666666666667; - bary[2] = 1. - bary[0] - bary[1]; + bary( 0, 0 ) = 0.1666666666667; + bary( 0, 1 ) = 0.6666666666667; + bary( 0, 2 ) = 1. - bary( 0, 0 ) - bary( 0, 1 ); // second barycentric point - bary[3] = 0.6666666666667; - bary[4] = 0.1666666666667; - bary[5] = 1. - bary[3] - bary[4]; + bary( 1, 0 ) = 0.6666666666667; + bary( 1, 1 ) = 0.1666666666667; + bary( 1, 2 ) = 1. - bary( 1, 0 ) - bary( 1, 1 ); // third barycentric point - bary[6] = 0.1666666666667; - bary[7] = 0.1666666666667; - bary[8] = 1. - bary[6] - bary[7]; + bary( 2, 0 ) = 0.1666666666667; + bary( 2, 1 ) = 0.1666666666667; + bary( 2, 2 ) = 1. - bary( 2, 0 ) - bary( 2, 1 ); break; case 3: @@ -166,36 +165,36 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ) } // first barycentric point - bary[0] = 0.0915762135098; - bary[1] = bary[0]; - bary[2] = 1. - bary[0] - bary[1]; + bary( 0, 0 ) = 0.0915762135098; + bary( 0, 1 ) = bary( 0, 0 ); + bary( 0, 2 ) = 1. - bary( 0, 0 ) - bary( 0, 1 ); // second barycentric point - bary[3] = 0.8168475729805; - bary[4] = bary[0]; - bary[5] = 1. - bary[3] - bary[4]; + bary( 1, 0 ) = 0.8168475729805; + bary( 1, 1 ) = bary( 0, 0 ); + bary( 1, 2 ) = 1. - bary( 1, 0 ) - bary( 1, 1 ); // third barycentric point - bary[6] = bary[0]; - bary[7] = bary[3]; - bary[8] = 1. - bary[6] - bary[7]; + bary( 2, 0 ) = bary( 0, 0 ); + bary( 2, 1 ) = bary( 1, 0 ); + bary( 2, 2 ) = 1. - bary( 2, 0 ) - bary( 2, 1 ); // fourth barycentric point - bary[9] = 0.1081030181681; - bary[10] = 0.4459484909160; - bary[11] = 1. - bary[9] - bary[10]; + bary( 3, 0 ) = 0.1081030181681; + bary( 3, 1 ) = 0.4459484909160; + bary( 3, 2 ) = 1. - bary( 3, 0 ) - bary( 3, 1 ); // fifth barycentric point - bary[12] = bary[10]; - bary[13] = bary[9]; - bary[14] = 1. - bary[12] - bary[13]; + bary( 4, 0 ) = bary( 3, 1 ); + bary( 4, 1 ) = bary( 3, 0 ); + bary( 4, 2 ) = 1. - bary( 4, 0 ) - bary( 4, 1 ); // sixth barycentric point - bary[15] = bary[10]; - bary[16] = bary[10]; - bary[17] = 1. - bary[15] - bary[16]; + bary( 5, 0 ) = bary( 3, 1 ); + bary( 5, 1 ) = bary( 3, 1 ); + bary( 5, 2 ) = 1. - bary( 5, 0 ) - bary( 5, 1 ); break; } // compute the vertex averaged centroid of the overlap polygon. Note // that the coordinates of the overlap polygon are always assumed to be // 3D - RealT xc[elem.dim]; + Array1D xc( elem.dim ); for ( int i = 0; i < elem.dim; ++i ) { xc[i] = 0.; } @@ -235,14 +234,14 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ) for ( int m = 0; m < numTriPoints; ++m ) { integ.wts[numTriPoints * k + m] *= 0.5 * area; integ.xy[( elem.dim * numTriPoints ) * k + ( elem.dim * m )] = - bary[elem.dim * m] * elem.overlapCoords[elem.dim * k] + - bary[elem.dim * m + 1] * elem.overlapCoords[elem.dim * ( k + 1 )] + bary[elem.dim * m + 2] * xc[0]; + bary( m, 0 ) * elem.overlapCoords[elem.dim * k] + bary( m, 1 ) * elem.overlapCoords[elem.dim * ( k + 1 )] + + bary( m, 2 ) * xc[0]; integ.xy[( elem.dim * numTriPoints ) * k + ( elem.dim * m ) + 1] = - bary[elem.dim * m] * elem.overlapCoords[elem.dim * k + 1] + - bary[elem.dim * m + 1] * elem.overlapCoords[elem.dim * ( k + 1 ) + 1] + bary[elem.dim * m + 2] * xc[1]; + bary( m, 0 ) * elem.overlapCoords[elem.dim * k + 1] + + bary( m, 1 ) * elem.overlapCoords[elem.dim * ( k + 1 ) + 1] + bary( m, 2 ) * xc[1]; integ.xy[( elem.dim * numTriPoints ) * k + ( elem.dim * m ) + 2] = - bary[elem.dim * m] * elem.overlapCoords[elem.dim * k + 2] + - bary[elem.dim * m + 1] * elem.overlapCoords[elem.dim * ( k + 1 ) + 2] + bary[elem.dim * m + 2] * xc[2]; + bary( m, 0 ) * elem.overlapCoords[elem.dim * k + 2] + + bary( m, 1 ) * elem.overlapCoords[elem.dim * ( k + 1 ) + 2] + bary( m, 2 ) * xc[2]; } // end loop over number of points per triangle } // end loop over (n-1) number of triangles @@ -265,14 +264,14 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ) for ( int i = 0; i < numTriPoints; ++i ) { integ.wts[numTriPoints * ( elem.numPolyVert - 1 ) + i] *= 0.5 * area; integ.xy[elem.dim * numTriPoints * ( elem.numPolyVert - 1 ) + ( elem.dim * i )] = - bary[elem.dim * i] * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 )] + - bary[elem.dim * i + 1] * elem.overlapCoords[0] + bary[elem.dim * i + 2] * xc[0]; + bary( i, 0 ) * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 )] + bary( i, 1 ) * elem.overlapCoords[0] + + bary( i, 2 ) * xc[0]; integ.xy[elem.dim * numTriPoints * ( elem.numPolyVert - 1 ) + ( elem.dim * i ) + 1] = - bary[elem.dim * i] * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 ) + 1] + - bary[elem.dim * i + 1] * elem.overlapCoords[1] + bary[elem.dim * i + 2] * xc[1]; + bary( i, 0 ) * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 ) + 1] + + bary( i, 1 ) * elem.overlapCoords[1] + bary( i, 2 ) * xc[1]; integ.xy[elem.dim * numTriPoints * ( elem.numPolyVert - 1 ) + ( elem.dim * i ) + 2] = - bary[elem.dim * i] * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 ) + 2] + - bary[elem.dim * i + 1] * elem.overlapCoords[2] + bary[elem.dim * i + 2] * xc[2]; + bary( i, 0 ) * elem.overlapCoords[elem.dim * ( elem.numPolyVert - 1 ) + 2] + + bary( i, 1 ) * elem.overlapCoords[2] + bary( i, 2 ) * xc[2]; } // end loop over numTriPoints for last triangle return; } diff --git a/src/tribol/integ/Integration.hpp b/src/tribol/integ/Integration.hpp index 98804dc6..61f90987 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_INTEG_INTEGRATION_HPP_ -#define SRC_INTEG_INTEGRATION_HPP_ +#ifndef SRC_TRIBOL_INTEG_INTEGRATION_HPP_ +#define SRC_TRIBOL_INTEG_INTEGRATION_HPP_ #include "tribol/common/Parameters.hpp" @@ -168,4 +168,4 @@ int NumTWBPointsPoly( SurfaceContactElem const& elem, int k ); int NumTWBPointsPerTri( int order ); } // end namespace tribol -#endif /* SRC_INTEG_INTEGRATION_HPP_ */ +#endif /* SRC_TRIBOL_INTEG_INTEGRATION_HPP_ */ diff --git a/src/tribol/interface/mfem_tribol.hpp b/src/tribol/interface/mfem_tribol.hpp index 1e0c1285..35388270 100644 --- a/src/tribol/interface/mfem_tribol.hpp +++ b/src/tribol/interface/mfem_tribol.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef MFEM_TRIBOL_HPP_ -#define MFEM_TRIBOL_HPP_ +#ifndef SRC_TRIBOL_INTERFACE_MFEM_TRIBOL_HPP_ +#define SRC_TRIBOL_INTERFACE_MFEM_TRIBOL_HPP_ #include "tribol/config.hpp" @@ -308,4 +308,4 @@ void saveRedecompMesh( int output_id ); #endif /* BUILD_REDECOMP */ -#endif /* MFEM_TRIBOL_HPP_ */ +#endif /* SRC_TRIBOL_INTERFACE_MFEM_TRIBOL_HPP_ */ diff --git a/src/tribol/interface/simple_tribol.hpp b/src/tribol/interface/simple_tribol.hpp index 8dee30b7..692aba9a 100644 --- a/src/tribol/interface/simple_tribol.hpp +++ b/src/tribol/interface/simple_tribol.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SIMPLE_TRIBOL_HPP_ -#define SIMPLE_TRIBOL_HPP_ +#ifndef SRC_TRIBOL_INTERFACE_SIMPLE_TRIBOL_HPP_ +#define SRC_TRIBOL_INTERFACE_SIMPLE_TRIBOL_HPP_ #include "tribol/common/Parameters.hpp" @@ -90,4 +90,4 @@ int GetSimpleCouplingCSR( int** I, int** J, double** vals, int* n_offsets, int* /// @} -#endif /* SIMPLE_TRIBOL_HPP_ */ +#endif /* SRC_TRIBOL_INTERFACE_SIMPLE_TRIBOL_HPP_ */ diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 14f0d766..4ac9e3be 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_HPP_ -#define TRIBOL_HPP_ +#ifndef SRC_TRIBOL_INTERFACE_TRIBOL_HPP_ +#define SRC_TRIBOL_INTERFACE_TRIBOL_HPP_ #include "tribol/common/ExecModel.hpp" #include "tribol/common/ArrayTypes.hpp" @@ -516,4 +516,4 @@ void finalize(); } /* namespace tribol */ -#endif /* TRIBOL_HPP_ */ +#endif /* SRC_TRIBOL_INTERFACE_TRIBOL_HPP_ */ diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 3f5d2d2e..bcce23a5 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -367,7 +367,9 @@ const ContactPlanePair& CouplingScheme::getContactPlanePair( IndexT id ) const break; } default: { - // no-op + // won't get here, but just return something to suppress warnings + return m_cg_pairs.getMortarPlane( id ); + break; } } // end switch } @@ -1815,7 +1817,9 @@ TRIBOL_HOST_DEVICE ContactPlanePair& CouplingScheme::Viewer::getContactPlanePair break; } default: { - // no-op + // won't get here, but just return something to suppress warnings + return m_cg_pairs.getMortarPlane( id ); + break; } } // end switch } diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index 72e69d27..35b2be63 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -2,8 +2,8 @@ // other Tribol Project Developers. See the top-level LICENSE file for details. // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_COUPLINGSCHEME_HPP_ -#define TRIBOL_COUPLINGSCHEME_HPP_ +#ifndef SRC_TRIBOL_MESH_COUPLINGSCHEME_HPP_ +#define SRC_TRIBOL_MESH_COUPLINGSCHEME_HPP_ // Tribol includes #include "tribol/common/BasicTypes.hpp" @@ -952,4 +952,4 @@ using CouplingSchemeManager = DataManager; } /* namespace tribol */ -#endif /* SRC_MESH_COUPLINGSCHEME_HPP_ */ +#endif /* SRC_TRIBOL_MESH_COUPLINGSCHEME_HPP_ */ diff --git a/src/tribol/mesh/InterfacePairs.hpp b/src/tribol/mesh/InterfacePairs.hpp index fdd38f3f..6c7aa4e2 100644 --- a/src/tribol/mesh/InterfacePairs.hpp +++ b/src/tribol/mesh/InterfacePairs.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_MESH_INTERFACE_PAIRS_HPP_ -#define SRC_MESH_INTERFACE_PAIRS_HPP_ +#ifndef SRC_TRIBOL_MESH_INTERFACE_PAIRS_HPP_ +#define SRC_TRIBOL_MESH_INTERFACE_PAIRS_HPP_ #include "tribol/common/BasicTypes.hpp" @@ -35,4 +35,4 @@ struct InterfacePair { } /* namespace tribol */ -#endif /* SRC_MESH_INTERFACE_PAIRS_HPP_ */ +#endif /* SRC_TRIBOL_MESH_INTERFACE_PAIRS_HPP_ */ diff --git a/src/tribol/mesh/MeshData.cpp b/src/tribol/mesh/MeshData.cpp index bc2b5be0..a93bbc9f 100644 --- a/src/tribol/mesh/MeshData.cpp +++ b/src/tribol/mesh/MeshData.cpp @@ -324,12 +324,12 @@ bool MeshData::computeFaceData( ExecutionMode exec_mode, ElemNormalMethod elem_n constexpr RealT nrml_mag_tol = 1.0e-15; // allocate m_c - m_c = Array2D( { m_dim, numberOfElements() }, m_allocator_id ); + m_c = ArrayT( { m_dim, numberOfElements() }, m_allocator_id ); // initialize m_c m_c.fill( 0.0 ); // allocate m_n - m_n = Array2D( { m_dim, numberOfElements() }, m_allocator_id ); + m_n = ArrayT( { m_dim, numberOfElements() }, m_allocator_id ); // initialize m_n m_n.fill( 0.0 ); @@ -476,7 +476,7 @@ RealT MeshData::computeEdgeLength( int faceId ) //------------------------------------------------------------------------------ void MeshData::allocateNodalNormals() { - m_node_n = Array2D( { m_dim, m_num_nodes }, m_allocator_id ); + m_node_n = ArrayT( { m_dim, m_num_nodes }, m_allocator_id ); m_node_n.fill( 0.0 ); } // end MeshData::allocateNodalNormals() diff --git a/src/tribol/mesh/MeshData.hpp b/src/tribol/mesh/MeshData.hpp index 4d764f31..6992d9d6 100644 --- a/src/tribol/mesh/MeshData.hpp +++ b/src/tribol/mesh/MeshData.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_MESH_MESHDATA_HPP_ -#define SRC_MESH_MESHDATA_HPP_ +#ifndef SRC_TRIBOL_MESH_MESHDATA_HPP_ +#define SRC_TRIBOL_MESH_MESHDATA_HPP_ #include @@ -677,13 +677,13 @@ class MeshData { MultiArrayView m_vel; ///< Nodal velocity MultiArrayView m_response; ///< Nodal responses (forces) - Array2D m_node_n; ///< Outward unit node normals + ArrayT m_node_n; ///< Outward unit node normals // Element field data Array2DView m_connectivity; ///< Element connectivity arrays - Array2D m_c; ///< Vertex averaged element centroids - Array2D m_n; ///< Outward unit element normals + ArrayT m_c; ///< Vertex averaged element centroids + ArrayT m_n; ///< Outward unit element normals Array1D m_face_radius; ///< Face radius used in low level proximity check Array1D m_area; ///< Element areas @@ -753,4 +753,4 @@ using MeshManager = DataManager; /// \a ostream operator to print a \a MeshData instance to \a os std::ostream& operator<<( std::ostream& os, const tribol::MeshData& md ); -#endif /* SRC_MESH_MESHDATA_HPP_ */ +#endif /* SRC_TRIBOL_MESH_MESHDATA_HPP_ */ diff --git a/src/tribol/mesh/MethodCouplingData.cpp b/src/tribol/mesh/MethodCouplingData.cpp index 2f1dd4ad..f52cb6e7 100644 --- a/src/tribol/mesh/MethodCouplingData.cpp +++ b/src/tribol/mesh/MethodCouplingData.cpp @@ -6,16 +6,10 @@ #include "tribol/mesh/MethodCouplingData.hpp" #include "tribol/common/Parameters.hpp" #include "tribol/mesh/MeshData.hpp" -#include "tribol/mesh/InterfacePairs.hpp" // Axom includes #include "axom/slic.hpp" -#include -#include -#include -#include - namespace tribol { //////////////////////////////////////////////// @@ -189,7 +183,7 @@ void MethodData::storeElemBlockJ( ArrayT&& blockJElemIds, const StackArray< // convert to mfem::DenseMatrix auto& block = blockJ( i, j ); // this DenseMatrix is a "view" of block - mfem::DenseMatrix block_densemat( block.data(), block.height(), block.width() ); + mfem::DenseMatrix block_densemat( const_cast( block.data() ), block.height(), block.width() ); // deep copy should happen here m_blockJ( blockIdxI, blockIdxJ ).push_back( block_densemat ); } @@ -282,29 +276,31 @@ void MortarData::assembleJacobian( SurfaceContactElem& elem, SparseMode s_mode ) // Mortar-Lagrange multiplier block (0, 2) this->m_smat->Add( i, j, elem.blockJ( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId] ); + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId] ); this->m_smat->Add( i + 1, j, elem.blockJ( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId + dimOffset] ); - this->m_smat->Add( - i + 2, j, - elem.blockJ( - static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId + 2 * dimOffset] ); // assume 3D for now + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId + dimOffset] ); + this->m_smat->Add( i + 2, j, + elem.blockJ( static_cast( BlockSpace::MORTAR ), + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId + 2 * dimOffset] ); // assume 3D for now // Nonmortar-Lagrange Multiplier block (1, 2) i = elem.dim * nonmortarNodeIdA; // nonmortar row contributions this->m_smat->Add( i, j, elem.blockJ( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId] ); + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId] ); this->m_smat->Add( i + 1, j, elem.blockJ( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId + dimOffset] ); - this->m_smat->Add( - i + 2, j, - elem.blockJ( - static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[localId + 2 * dimOffset] ); // assume 3D for now + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId + dimOffset] ); + this->m_smat->Add( i + 2, j, + elem.blockJ( static_cast( BlockSpace::NONMORTAR ), + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[localId + 2 * dimOffset] ); // assume 3D for now //////////////////////////////////////////////////////////////// // Assemble Jgu contributions (lower-left off-diagonal block) // @@ -321,27 +317,31 @@ void MortarData::assembleJacobian( SurfaceContactElem& elem, SparseMode s_mode ) // Lagrange-multiplier-mortar block (2, 0) this->m_smat->Add( i, j, elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[localId] ); + static_cast( BlockSpace::MORTAR ) ) + .data()[localId] ); this->m_smat->Add( i, j + 1, elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[localId + dimOffset] ); - this->m_smat->Add( - i, j + 2, - elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[localId + 2 * dimOffset] ); // assume 3D for now + static_cast( BlockSpace::MORTAR ) ) + .data()[localId + dimOffset] ); + this->m_smat->Add( i, j + 2, + elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), + static_cast( BlockSpace::MORTAR ) ) + .data()[localId + 2 * dimOffset] ); // assume 3D for now // Lagrange multiplier-nonmortar block (2, 1) j = elem.dim * nonmortarNodeIdA; // nonmortar column contributions this->m_smat->Add( i, j, elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[localId] ); + static_cast( BlockSpace::NONMORTAR ) ) + .data()[localId] ); this->m_smat->Add( i, j + 1, elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[localId + dimOffset] ); - this->m_smat->Add( - i, j + 2, - elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[localId + 2 * dimOffset] ); // assume 3D for now + static_cast( BlockSpace::NONMORTAR ) ) + .data()[localId + dimOffset] ); + this->m_smat->Add( i, j + 2, + elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), + static_cast( BlockSpace::NONMORTAR ) ) + .data()[localId + 2 * dimOffset] ); // assume 3D for now ////////////////////////////////////////////////// // Assemble Jru contributions (matrix 11-block) // diff --git a/src/tribol/mesh/MethodCouplingData.hpp b/src/tribol/mesh/MethodCouplingData.hpp index 0e937d84..1ca558e0 100644 --- a/src/tribol/mesh/MethodCouplingData.hpp +++ b/src/tribol/mesh/MethodCouplingData.hpp @@ -3,12 +3,12 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_MESH_METHODCOUPLINGDATA_HPP_ -#define SRC_MESH_METHODCOUPLINGDATA_HPP_ +#ifndef SRC_TRIBOL_MESH_METHODCOUPLINGDATA_HPP_ +#define SRC_TRIBOL_MESH_METHODCOUPLINGDATA_HPP_ #include "tribol/common/ArrayTypes.hpp" -#include "tribol/common/Parameters.hpp" #include "tribol/common/Containers.hpp" +#include "tribol/common/Parameters.hpp" #include "tribol/mesh/MeshData.hpp" // Axom includes @@ -350,4 +350,4 @@ class MortarData : public MethodData { }; } // end namespace tribol -#endif /* SRC_MESH_METHODCOUPLINGDATA_HPP_ */ +#endif /* SRC_TRIBOL_MESH_METHODCOUPLINGDATA_HPP_ */ diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index de809434..6ac22dfc 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_MESH_MFEMDATA_HPP_ -#define SRC_MESH_MFEMDATA_HPP_ +#ifndef SRC_TRIBOL_MESH_MFEMDATA_HPP_ +#define SRC_TRIBOL_MESH_MFEMDATA_HPP_ // Tribol config include #include "tribol/config.hpp" @@ -1749,4 +1749,4 @@ class MfemJacobianData { #endif /* BUILD_REDECOMP */ -#endif /* SRC_MESH_MFEMDATA_HPP_ */ +#endif /* SRC_TRIBOL_MESH_MFEMDATA_HPP_ */ diff --git a/src/tribol/physics/AlignedMortar.cpp b/src/tribol/physics/AlignedMortar.cpp index efa83c7c..a9aac9ec 100644 --- a/src/tribol/physics/AlignedMortar.cpp +++ b/src/tribol/physics/AlignedMortar.cpp @@ -14,14 +14,9 @@ #include "tribol/common/Parameters.hpp" #include "tribol/integ/Integration.hpp" #include "tribol/integ/FE.hpp" -#include "tribol/utils/ContactPlaneOutput.hpp" #include "tribol/utils/Algorithm.hpp" #include "tribol/utils/Math.hpp" -#include -#include -#include - namespace tribol { void ComputeAlignedMortarWeights( SurfaceContactElem& elem ) @@ -106,7 +101,7 @@ void ComputeNodalGap( SurfaceContactElem& elem ) // loop over nodes on nonmortar side for ( int a = 0; a < elem.numFaceVert; ++a ) { // get global nonmortar node number from connectivity - RealT nrml_a[elem.dim]; + Array1D nrml_a( elem.dim ); int glbId = nonmortarConn[elem.numFaceVert * elem.faceId2 + a]; nrml_a[0] = nonmortarMesh.getNodalNormals()[0][glbId]; nrml_a[1] = nonmortarMesh.getNodalNormals()[1][glbId]; @@ -170,8 +165,8 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) const IndexT numNodesPerFace = mortarMesh.numberOfNodesPerElement(); // arrays to store face coords - RealT mortarX[dim * numNodesPerFace]; - RealT nonmortarX[dim * numNodesPerFace]; + Array2D mortarX( numNodesPerFace, dim ); + Array2D nonmortarX( numNodesPerFace, dim ); //////////////////////////// // compute nonmortar gaps // @@ -187,7 +182,7 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) auto& cg_pairs = cs->getCompGeom(); auto& plane = cg_pairs.getAlignedMortarPlane( cpID ); - RealT overlapX[dim * plane.m_numPolyVert]; + Array2D overlapX( plane.m_numPolyVert, dim ); // get pair indices IndexT index1 = pair.m_element_id1; @@ -198,18 +193,18 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) // onto the common plane, since the aligned mortar gap // calculation uses the current configuration nodal coordinates // themselves - plane.getFace1Coords( &mortarX[0], numNodesPerFace ); - plane.getFace2Coords( &nonmortarX[0], numNodesPerFace ); + plane.getFace1Coords( mortarX.data(), numNodesPerFace ); + plane.getFace2Coords( nonmortarX.data(), numNodesPerFace ); - plane.getOverlapVertices( &overlapX[0] ); + plane.getOverlapVertices( overlapX.data() ); // instantiate SurfaceContactElem struct. Note, this is done with // the projected area of overlap, but with the actual current // configuration face coordinates. We need the current // configuration face coordinates here in order to correctly // compute the mortar gaps. - SurfaceContactElem elem_for_gap( dim, mortarX, nonmortarX, overlapX, numNodesPerFace, plane.m_numPolyVert, - &mortarMesh, &nonmortarMesh, index1, index2 ); + SurfaceContactElem elem_for_gap( dim, mortarX.data(), nonmortarX.data(), overlapX.data(), numNodesPerFace, + plane.m_numPolyVert, &mortarMesh, &nonmortarMesh, index1, index2 ); ///////////////////////// // compute mortar gaps // @@ -343,16 +338,16 @@ int ApplyNormal( CouplingScheme* cs ) auto& plane = cg_pairs.getAlignedMortarPlane( cpID ); // get projected face coords and overlap coords - RealT mortarX_bar[dim * numNodesPerFace]; - RealT nonmortarX_bar[dim * numNodesPerFace]; - RealT overlapX[dim * plane.m_numPolyVert]; - plane.getFace1ProjectedCoords( &mortarX_bar[0], numNodesPerFace ); - plane.getFace2ProjectedCoords( &nonmortarX_bar[0], numNodesPerFace ); - plane.getOverlapVertices( &overlapX[0] ); + Array2D mortarX_bar( numNodesPerFace, dim ); + Array2D nonmortarX_bar( numNodesPerFace, dim ); + Array2D overlapX( plane.m_numPolyVert, dim ); + plane.getFace1ProjectedCoords( mortarX_bar.data(), numNodesPerFace ); + plane.getFace2ProjectedCoords( nonmortarX_bar.data(), numNodesPerFace ); + plane.getOverlapVertices( overlapX.data() ); // instantiate a new surface contact element with projected face // coordinates - SurfaceContactElem elem_for_jac( dim, &mortarX_bar[0], &nonmortarX_bar[0], &overlapX[0], numNodesPerFace, + SurfaceContactElem elem_for_jac( dim, mortarX_bar.data(), nonmortarX_bar.data(), overlapX.data(), numNodesPerFace, plane.m_numPolyVert, &mortarMesh, &nonmortarMesh, index1, index2 ); // HAVE TO set the number of active constraints. For now set to diff --git a/src/tribol/physics/AlignedMortar.hpp b/src/tribol/physics/AlignedMortar.hpp index eee60f7d..c9998149 100644 --- a/src/tribol/physics/AlignedMortar.hpp +++ b/src/tribol/physics/AlignedMortar.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_PHYSICS_ALIGNEDMORTAR_HPP_ -#define SRC_PHYSICS_ALIGNEDMORTAR_HPP_ +#ifndef SRC_TRIBOL_PHYSICS_ALIGNEDMORTAR_HPP_ +#define SRC_TRIBOL_PHYSICS_ALIGNEDMORTAR_HPP_ #include "Mortar.hpp" #include "Physics.hpp" @@ -132,4 +132,4 @@ void ComputeAlignedMortarJacobian( SurfaceContactElem& elem ); } // namespace tribol -#endif /* SRC_PHYSICS_ALIGNEDMORTAR_HPP_ */ +#endif /* SRC_TRIBOL_PHYSICS_ALIGNEDMORTAR_HPP_ */ diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 6f45d15b..1aaf5b59 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -427,11 +427,9 @@ int ApplyTangential( CouplingScheme* /////////////////////////////// // loop over interface pairs // /////////////////////////////// - ArrayT err_data( { 0 }, cs->getAllocatorId() ); - ArrayViewT err = err_data; auto cs_view = cs->getView(); const auto num_pairs = cs->getNumActivePairs(); - forAllExec( cs->getExecutionMode(), num_pairs, [cs_view, err] TRIBOL_HOST_DEVICE( IndexT i ) { + forAllExec( cs->getExecutionMode(), num_pairs, [cs_view] TRIBOL_HOST_DEVICE( IndexT i ) { auto& cg_view = cs_view.getCompGeomView(); auto& plane = cg_view.getCommonPlane( i ); @@ -615,8 +613,7 @@ int ApplyTangential( CouplingScheme* } // end for loop over face nodes } ); - ArrayT err_host( err_data ); - return err_host[0]; + return 0; } //------------------------------------------------------------------------------ diff --git a/src/tribol/physics/CommonPlane.hpp b/src/tribol/physics/CommonPlane.hpp index 560c04c1..2b1d82ac 100644 --- a/src/tribol/physics/CommonPlane.hpp +++ b/src/tribol/physics/CommonPlane.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_PHYSICS_COMMONPLANE_HPP_ -#define SRC_PHYSICS_COMMONPLANE_HPP_ +#ifndef SRC_TRIBOL_PHYSICS_COMMONPLANE_HPP_ +#define SRC_TRIBOL_PHYSICS_COMMONPLANE_HPP_ #include "Physics.hpp" @@ -50,4 +50,4 @@ int ApplyTangential( CouplingScheme* } // end namespace tribol -#endif /* SRC_PHYSICS_COMMONPLANE_HPP_ */ +#endif /* SRC_TRIBOL_PHYSICS_COMMONPLANE_HPP_ */ diff --git a/src/tribol/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 94433da7..91e85f1d 100644 --- a/src/tribol/physics/Mortar.cpp +++ b/src/tribol/physics/Mortar.cpp @@ -11,23 +11,16 @@ #include "tribol/geom/CompGeom.hpp" #include "tribol/geom/GeomUtilities.hpp" #include "tribol/geom/NodalNormal.hpp" +#include "tribol/common/ArrayTypes.hpp" +#include "tribol/common/Enzyme.hpp" #include "tribol/common/Parameters.hpp" #include "tribol/integ/Integration.hpp" #include "tribol/integ/FE.hpp" -#include "tribol/utils/ContactPlaneOutput.hpp" #include "tribol/utils/Math.hpp" -#include "tribol/utils/Algorithm.hpp" // Axom includes #include "axom/slic.hpp" -#include -#include - -#ifdef TRIBOL_USE_ENZYME -#include "mfem/general/enzyme.hpp" -#endif - namespace tribol { void ComputeMortarWeights( SurfaceContactElem& elem ) @@ -43,12 +36,12 @@ void ComputeMortarWeights( SurfaceContactElem& elem ) // TWBPolyInt( elem, integ, 3 ); // get individual arrays of coordinates for each face - RealT x1[elem.numFaceVert]; - RealT y1[elem.numFaceVert]; - RealT z1[elem.numFaceVert]; - RealT x2[elem.numFaceVert]; - RealT y2[elem.numFaceVert]; - RealT z2[elem.numFaceVert]; + Array1D x1( elem.numFaceVert ); + Array1D y1( elem.numFaceVert ); + Array1D z1( elem.numFaceVert ); + Array1D x2( elem.numFaceVert ); + Array1D y2( elem.numFaceVert ); + Array1D z2( elem.numFaceVert ); for ( int i = 0; i < elem.numFaceVert; ++i ) { x1[i] = elem.faceCoords1[elem.dim * i]; @@ -85,14 +78,14 @@ void ComputeMortarWeights( SurfaceContactElem& elem ) RealT xp[3] = { integ.xy[elem.dim * ip], integ.xy[elem.dim * ip + 1], integ.xy[elem.dim * ip + 2] }; RealT xi[2] = { 0., 0. }; - InvIso( xp, x1, y1, z1, elem.numFaceVert, xi ); + InvIso( xp, x1.data(), y1.data(), z1.data(), elem.numFaceVert, xi ); if ( elem.numFaceVert == 4 ) { LinIsoQuadShapeFunc( xi[0], xi[1], a, phiMortarA ); } else if ( elem.numFaceVert == 3 ) { LinIsoTriShapeFunc( xi[0], xi[1], a, phiMortarA ); } - InvIso( xp, x2, y2, z2, elem.numFaceVert, xi ); + InvIso( xp, x2.data(), y2.data(), z2.data(), elem.numFaceVert, xi ); if ( elem.numFaceVert == 4 ) { LinIsoQuadShapeFunc( xi[0], xi[1], a, phiNonmortarA ); LinIsoQuadShapeFunc( xi[0], xi[1], b, phiNonmortarB ); @@ -148,7 +141,7 @@ void ComputeNodalGap( SurfaceContactElem& elem ) RealT g2 = 0.; // get global nonmortar node number from connectivity - RealT nrml_a[elem.dim]; + Array1D nrml_a( elem.dim ); int glbId = nonmortarConn[elem.numFaceVert * elem.faceId2 + a]; nrml_a[0] = nonmortarMesh.getNodalNormals()[0][glbId]; nrml_a[1] = nonmortarMesh.getNodalNormals()[1][glbId]; @@ -166,8 +159,8 @@ void ComputeNodalGap( SurfaceContactElem& elem ) RealT nab_1 = elem.getNonmortarMortarWt( a, b ); // nonmortar-mortar weight RealT nab_2 = elem.getNonmortarNonmortarWt( a, b ); // nonmortar-nonmortar weight - g1 += dotProd( &nrml_a[0], &elem.faceCoords1[elem.dim * b], elem.dim ) * nab_1; - g2 += dotProd( &nrml_a[0], &elem.faceCoords2[elem.dim * b], elem.dim ) * nab_2; + g1 += dotProd( nrml_a.data(), &elem.faceCoords1[elem.dim * b], elem.dim ) * nab_1; + g2 += dotProd( nrml_a.data(), &elem.faceCoords2[elem.dim * b], elem.dim ) * nab_2; } // store local gap @@ -202,11 +195,10 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) // declare local variables to hold face nodal coordinates // and overlap vertex coordinates int const dim = cs->spatialDimension(); - IndexT size = dim * numNodesPerFace; - RealT mortarX[size]; - RealT nonmortarX[size]; - RealT mortarX_bar[size]; - RealT nonmortarX_bar[size]; + Array2D mortarX( numNodesPerFace, dim ); + Array2D nonmortarX( numNodesPerFace, dim ); + Array2D mortarX_bar( numNodesPerFace, dim ); + Array2D nonmortarX_bar( numNodesPerFace, dim ); //////////////////////////////////////////////////////////////////// // compute nonmortar gaps to determine active set of contact dofs // @@ -222,7 +214,7 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) auto& cg_pairs = cs->getCompGeom(); auto& plane = cg_pairs.getMortarPlane( cpID ); - RealT overlapX[dim * plane.m_numPolyVert]; + Array2D overlapX( plane.m_numPolyVert, dim ); // get pair indices IndexT index1 = pair.m_element_id1; @@ -230,19 +222,19 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) // populate the current configuration nodal coordinates for the // two faces; stored on the contact plane object - plane.getFace1Coords( &mortarX[0], numNodesPerFace ); - plane.getFace2Coords( &nonmortarX[0], numNodesPerFace ); + plane.getFace1Coords( mortarX.data(), numNodesPerFace ); + plane.getFace2Coords( nonmortarX.data(), numNodesPerFace ); // get face coordinates projected onto contact plane - plane.getFace1ProjectedCoords( &mortarX_bar[0], numNodesPerFace ); - plane.getFace2ProjectedCoords( &nonmortarX_bar[0], numNodesPerFace ); + plane.getFace1ProjectedCoords( mortarX_bar.data(), numNodesPerFace ); + plane.getFace2ProjectedCoords( nonmortarX_bar.data(), numNodesPerFace ); // get overlap vertices - plane.getOverlapVertices( &overlapX[0] ); + plane.getOverlapVertices( overlapX.data() ); // instantiate contact surface element for purposes of computing // mortar weights. Note, this uses projected face coords - SurfaceContactElem elem( dim, &mortarX_bar[0], &nonmortarX_bar[0], &overlapX[0], numNodesPerFace, + SurfaceContactElem elem( dim, mortarX_bar.data(), nonmortarX_bar.data(), overlapX.data(), numNodesPerFace, plane.m_numPolyVert, &mortarMesh, &nonmortarMesh, index1, index2 ); // compute the mortar weights to be stored on the surface @@ -252,8 +244,8 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) // compute mortar gaps. Note, we have to now use current configuration // nodal coordinates on the contact element - elem.faceCoords1 = &mortarX[0]; - elem.faceCoords2 = &nonmortarX[0]; + elem.faceCoords1 = mortarX.data(); + elem.faceCoords2 = nonmortarX.data(); ComputeNodalGap( elem ); @@ -330,9 +322,8 @@ int ApplyNormal( CouplingScheme* cs ) } // declare local variables to hold projected face nodal coordinates - IndexT size = dim * numNodesPerFace; - RealT mortarX_bar[size]; - RealT nonmortarX_bar[size]; + Array2D mortarX_bar( numNodesPerFace, dim ); + Array2D nonmortarX_bar( numNodesPerFace, dim ); //////////////////////////////////////////////////////////////// // // @@ -350,22 +341,22 @@ int ApplyNormal( CouplingScheme* cs ) auto& cg_pairs = cs->getCompGeom(); auto& plane = cg_pairs.getMortarPlane( cpID ); - RealT overlapX[dim * plane.m_numPolyVert]; + Array2D overlapX( plane.m_numPolyVert, dim ); // get pair indices IndexT index1 = pair.m_element_id1; IndexT index2 = pair.m_element_id2; // get face coordinates projected onto contact plane - plane.getFace1ProjectedCoords( &mortarX_bar[0], numNodesPerFace ); - plane.getFace2ProjectedCoords( &nonmortarX_bar[0], numNodesPerFace ); + plane.getFace1ProjectedCoords( mortarX_bar.data(), numNodesPerFace ); + plane.getFace2ProjectedCoords( nonmortarX_bar.data(), numNodesPerFace ); // get overlap coords - plane.getOverlapVertices( &overlapX[0] ); + plane.getOverlapVertices( overlapX.data() ); // instantiate contact surface element for purposes of computing // mortar weights. Note, this uses projected face coords - SurfaceContactElem elem( dim, &mortarX_bar[0], &nonmortarX_bar[0], &overlapX[0], numNodesPerFace, + SurfaceContactElem elem( dim, mortarX_bar.data(), nonmortarX_bar.data(), overlapX.data(), numNodesPerFace, plane.m_numPolyVert, &mortarMesh, &nonmortarMesh, index1, index2 ); ////////////////////////////////// @@ -467,7 +458,7 @@ void ComputeResidualJacobian( SurfaceContactElem& elem ) for ( int b = 0; b < elem.numFaceVert; ++b ) { // get global nonmortar node id to index into nodal normals on // nonmortar mesh - RealT nrml_b[elem.dim]; + Array1D nrml_b( elem.dim ); int glbId = nonmortarConn[elem.numFaceVert * elem.faceId2 + b]; // We assemble ALL nonmortar node contributions, even if gap is in separation. @@ -491,24 +482,23 @@ void ComputeResidualJacobian( SurfaceContactElem& elem ) // Fill block (0, 2) int elem_xdof = elem.getJacobianIndex( SurfaceContactElem::JrpBlock, a, b ); int dim_offset = elem.getJacobianDimOffset( SurfaceContactElem::JrpBlock ); - elem.blockJ( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof] += nrml_b[0] * n_mortar_b; - elem.blockJ( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof + dim_offset] += - nrml_b[1] * n_mortar_b; - elem.blockJ( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof + 2 * dim_offset] += - nrml_b[2] * n_mortar_b; + elem.blockJ( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof] += nrml_b[0] * n_mortar_b; + elem.blockJ( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof + dim_offset] += nrml_b[1] * n_mortar_b; + elem.blockJ( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof + 2 * dim_offset] += nrml_b[2] * n_mortar_b; // Fill block (1, 2) elem.blockJ( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof] -= nrml_b[0] * n_nonmortar_b; + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof] -= nrml_b[0] * n_nonmortar_b; elem.blockJ( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof + dim_offset] -= - nrml_b[1] * n_nonmortar_b; + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof + dim_offset] -= nrml_b[1] * n_nonmortar_b; elem.blockJ( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) )[elem_xdof + 2 * dim_offset] -= - nrml_b[2] * n_nonmortar_b; + static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) + .data()[elem_xdof + 2 * dim_offset] -= nrml_b[2] * n_nonmortar_b; } // end loop over b nodes @@ -529,7 +519,7 @@ void ComputeConstraintJacobian( SurfaceContactElem& elem for ( int a = 0; a < elem.numFaceVert; ++a ) { // get global nonmortar node id to index into nodal normals on // nonmortar mesh - RealT nrml_a[elem.dim]; + Array1D nrml_a( elem.dim ); int glbId = nonmortarConn[elem.numFaceVert * elem.faceId2 + a]; // We assemble ALL nonmortar node contributions even if gap is in separation. @@ -556,21 +546,23 @@ void ComputeConstraintJacobian( SurfaceContactElem& elem // Fill block (2, 0) int dim_offset = elem.getJacobianDimOffset( SurfaceContactElem::JguBlock ); int elem_xdof = elem.getJacobianIndex( SurfaceContactElem::JguBlock, a, b ); - elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[elem_xdof] += nrml_a[0] * n_mortar_a; - elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[elem_xdof + dim_offset] += nrml_a[1] * n_mortar_a; - elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) )[elem_xdof + 2 * dim_offset] += nrml_a[2] * n_mortar_a; + elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), static_cast( BlockSpace::MORTAR ) ) + .data()[elem_xdof] += nrml_a[0] * n_mortar_a; + elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), static_cast( BlockSpace::MORTAR ) ) + .data()[elem_xdof + dim_offset] += nrml_a[1] * n_mortar_a; + elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), static_cast( BlockSpace::MORTAR ) ) + .data()[elem_xdof + 2 * dim_offset] += nrml_a[2] * n_mortar_a; // Fill block (2, 1) elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[elem_xdof] -= nrml_a[0] * n_nonmortar_a; + static_cast( BlockSpace::NONMORTAR ) ) + .data()[elem_xdof] -= nrml_a[0] * n_nonmortar_a; elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[elem_xdof + dim_offset] -= nrml_a[1] * n_nonmortar_a; + static_cast( BlockSpace::NONMORTAR ) ) + .data()[elem_xdof + dim_offset] -= nrml_a[1] * n_nonmortar_a; elem.blockJ( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) )[elem_xdof + 2 * dim_offset] -= - nrml_a[2] * n_nonmortar_a; + static_cast( BlockSpace::NONMORTAR ) ) + .data()[elem_xdof + 2 * dim_offset] -= nrml_a[2] * n_nonmortar_a; } // end loop over b nodes @@ -1019,15 +1011,15 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* x1_dot[i] = 1.0; // clang-format off __enzyme_fwddiff((void*)ComputeMortarForceEnzyme, - enzyme_dup, x1, x1_dot, - enzyme_const, n1, - enzyme_const, p1, - enzyme_dup, f1, &df1dx1[size1*3*i], - enzyme_dup, g1, &dg1dx1[size1*i], - enzyme_const, size1, - enzyme_const, x2, - enzyme_dup, f2, &df2dx1[size1*3*i], - enzyme_const, size2); + TRIBOL_ENZYME_DUP, x1, x1_dot, + TRIBOL_ENZYME_CONST, n1, + TRIBOL_ENZYME_CONST, p1, + TRIBOL_ENZYME_DUP, f1, &df1dx1[size1*3*i], + TRIBOL_ENZYME_DUP, g1, &dg1dx1[size1*i], + TRIBOL_ENZYME_CONST, size1, + TRIBOL_ENZYME_CONST, x2, + TRIBOL_ENZYME_DUP, f2, &df2dx1[size1*3*i], + TRIBOL_ENZYME_CONST, size2); // clang-format on x1_dot[i] = 0.0; } @@ -1036,15 +1028,15 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* n1_dot[i] = 1.0; // clang-format off __enzyme_fwddiff((void*)ComputeMortarForceEnzyme, - enzyme_const, x1, - enzyme_dup, n1, n1_dot, - enzyme_const, p1, - enzyme_dup, f1, &df1dn1[size1*3*i], - enzyme_dup, g1, &dg1dn1[size1*i], - enzyme_const, size1, - enzyme_const, x2, - enzyme_dup, f2, &df2dn1[size1*3*i], - enzyme_const, size2); + TRIBOL_ENZYME_CONST, x1, + TRIBOL_ENZYME_DUP, n1, n1_dot, + TRIBOL_ENZYME_CONST, p1, + TRIBOL_ENZYME_DUP, f1, &df1dn1[size1*3*i], + TRIBOL_ENZYME_DUP, g1, &dg1dn1[size1*i], + TRIBOL_ENZYME_CONST, size1, + TRIBOL_ENZYME_CONST, x2, + TRIBOL_ENZYME_DUP, f2, &df2dn1[size1*3*i], + TRIBOL_ENZYME_CONST, size2); // clang-format on n1_dot[i] = 0.0; } @@ -1053,15 +1045,15 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* p1_dot[i] = 1.0; // clang-format off __enzyme_fwddiff((void*)ComputeMortarForceEnzyme, - enzyme_const, x1, - enzyme_const, n1, - enzyme_dup, p1, p1_dot, - enzyme_dup, f1, &df1dp1[size1*3*i], - enzyme_const, g1, - enzyme_const, size1, - enzyme_const, x2, - enzyme_dup, f2, &df2dp1[size1*3*i], - enzyme_const, size2); + TRIBOL_ENZYME_CONST, x1, + TRIBOL_ENZYME_CONST, n1, + TRIBOL_ENZYME_DUP, p1, p1_dot, + TRIBOL_ENZYME_DUP, f1, &df1dp1[size1*3*i], + TRIBOL_ENZYME_CONST, g1, + TRIBOL_ENZYME_CONST, size1, + TRIBOL_ENZYME_CONST, x2, + TRIBOL_ENZYME_DUP, f2, &df2dp1[size1*3*i], + TRIBOL_ENZYME_CONST, size2); // clang-format on p1_dot[i] = 0.0; } @@ -1070,15 +1062,15 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* x2_dot[i] = 1.0; // clang-format off __enzyme_fwddiff((void*)ComputeMortarForceEnzyme, - enzyme_const, x1, - enzyme_const, n1, - enzyme_const, p1, - enzyme_dup, f1, &df1dx2[size2*3*i], - enzyme_dup, g1, &dg1dx2[size2*i], - enzyme_const, size1, - enzyme_dup, x2, x2_dot, - enzyme_dup, f2, &df2dx2[size2*3*i], - enzyme_const, size2); + TRIBOL_ENZYME_CONST, x1, + TRIBOL_ENZYME_CONST, n1, + TRIBOL_ENZYME_CONST, p1, + TRIBOL_ENZYME_DUP, f1, &df1dx2[size2*3*i], + TRIBOL_ENZYME_DUP, g1, &dg1dx2[size2*i], + TRIBOL_ENZYME_CONST, size1, + TRIBOL_ENZYME_DUP, x2, x2_dot, + TRIBOL_ENZYME_DUP, f2, &df2dx2[size2*3*i], + TRIBOL_ENZYME_CONST, size2); // clang-format on x2_dot[i] = 0.0; } @@ -1105,9 +1097,8 @@ int GetMethodData( CouplingScheme* cs ) auto nonmortarMesh = cs->getMesh2().getView(); IndexT const numNodesPerFace = mortarMesh.numberOfNodesPerElement(); - const int size = dim * numNodesPerFace; - RealT mortarX_bar[size]; - RealT nonmortarX_bar[size]; + Array2D mortarX_bar( numNodesPerFace, dim ); + Array2D nonmortarX_bar( numNodesPerFace, dim ); int numRows = cs->getNumTotalNodes(); static_cast( cs->getMethodData() )->allocateMfemSparseMatrix( numRows ); @@ -1129,22 +1120,22 @@ int GetMethodData( CouplingScheme* cs ) auto& cg_pairs = cs->getCompGeom(); auto& plane = cg_pairs.getMortarPlane( cpID ); - RealT overlapX[dim * plane.m_numPolyVert]; + Array2D overlapX( plane.m_numPolyVert, dim ); // get pair indices IndexT index1 = pair.m_element_id1; IndexT index2 = pair.m_element_id2; // get face coordinates projected onto contact plane - plane.getFace1ProjectedCoords( &mortarX_bar[0], numNodesPerFace ); - plane.getFace2ProjectedCoords( &nonmortarX_bar[0], numNodesPerFace ); + plane.getFace1ProjectedCoords( mortarX_bar.data(), numNodesPerFace ); + plane.getFace2ProjectedCoords( nonmortarX_bar.data(), numNodesPerFace ); // construct array of polygon overlap vertex coordinates - plane.getOverlapVertices( &overlapX[0] ); + plane.getOverlapVertices( overlapX.data() ); // instantiate contact surface element for purposes of computing // mortar weights. Note, this uses projected face coords - SurfaceContactElem elem( dim, &mortarX_bar[0], &nonmortarX_bar[0], &overlapX[0], numNodesPerFace, + SurfaceContactElem elem( dim, mortarX_bar.data(), nonmortarX_bar.data(), overlapX.data(), numNodesPerFace, plane.m_numPolyVert, &mortarMesh, &nonmortarMesh, index1, index2 ); // compute the mortar weights to be stored on the surface diff --git a/src/tribol/physics/Mortar.hpp b/src/tribol/physics/Mortar.hpp index df98efa9..efec3a9f 100644 --- a/src/tribol/physics/Mortar.hpp +++ b/src/tribol/physics/Mortar.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_PHYSICS_MORTAR_HPP_ -#define SRC_PHYSICS_MORTAR_HPP_ +#ifndef SRC_TRIBOL_PHYSICS_MORTAR_HPP_ +#define SRC_TRIBOL_PHYSICS_MORTAR_HPP_ #include "tribol/common/Parameters.hpp" #include "Physics.hpp" @@ -247,4 +247,4 @@ int GetMethodData( CouplingScheme* cs ); } // namespace tribol -#endif /* SRC_PHYSICS_MORTAR_HPP_ */ +#endif /* SRC_TRIBOL_PHYSICS_MORTAR_HPP_ */ diff --git a/src/tribol/physics/Physics.hpp b/src/tribol/physics/Physics.hpp index 56cf2254..753237a9 100644 --- a/src/tribol/physics/Physics.hpp +++ b/src/tribol/physics/Physics.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_PHYSICS_PHYSICS_HPP_ -#define SRC_PHYSICS_PHYSICS_HPP_ +#ifndef SRC_TRIBOL_PHYSICS_PHYSICS_HPP_ +#define SRC_TRIBOL_PHYSICS_PHYSICS_HPP_ #include "tribol/common/Parameters.hpp" @@ -63,4 +63,4 @@ int GetMethodData( CouplingScheme* cs ); } // namespace tribol -#endif /* SRC_PHYSICS_PHYSICS_HPP_ */ +#endif /* SRC_TRIBOL_PHYSICS_PHYSICS_HPP_ */ diff --git a/src/tribol/search/InterfacePairFinder.hpp b/src/tribol/search/InterfacePairFinder.hpp index 46231c74..e9886b94 100644 --- a/src/tribol/search/InterfacePairFinder.hpp +++ b/src/tribol/search/InterfacePairFinder.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ -#define TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ +#ifndef SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ +#define SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ #include "tribol/common/Parameters.hpp" #include "tribol/mesh/MeshData.hpp" @@ -57,4 +57,4 @@ class InterfacePairFinder { } // end namespace tribol -#endif /* TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ */ +#endif /* SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ */ diff --git a/src/tribol/utils/Algorithm.hpp b/src/tribol/utils/Algorithm.hpp index 13379115..b6c33b6d 100644 --- a/src/tribol/utils/Algorithm.hpp +++ b/src/tribol/utils/Algorithm.hpp @@ -3,11 +3,14 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_UTILS_ALGORITHM_HPP_ -#define SRC_UTILS_ALGORITHM_HPP_ +#ifndef SRC_TRIBOL_UTILS_ALGORITHM_HPP_ +#define SRC_TRIBOL_UTILS_ALGORITHM_HPP_ #include "tribol/common/ArrayTypes.hpp" +#include +#include + namespace tribol { namespace algorithm { @@ -127,8 +130,71 @@ TRIBOL_HOST_DEVICE void transpose( const ArrayT& in, ArrayT +TRIBOL_HOST_DEVICE void bubbleSort( Container& c, Compare comp ) +{ + using std::swap; + + const auto size = c.size(); + if ( size < 2 ) { + return; + } + + for ( decltype( size ) i = 0; i < size - 1; ++i ) { + for ( decltype( size ) j = 0; j < size - i - 1; ++j ) { + if ( comp( c[j + 1], c[j] ) ) { + swap( c[j], c[j + 1] ); + } + } + } +} + +/** + * @brief Computes the product of all elements in a container. + * + * @tparam Container Type of container + * @param c Container to compute product of + * @return The product of all elements in the container + */ +template +TRIBOL_HOST_DEVICE typename Container::value_type product( const Container& c ) +{ + typename Container::value_type result = 1; + for ( auto val : c ) { + result *= val; + } + return result; +} + +/** + * @brief Finds the minimum element in a container. + * + * @tparam Container Type of container + * @param c Container to find the minimum element of + * @return The minimum element in the container + */ +template +TRIBOL_HOST_DEVICE typename Container::value_type min( const Container& c ) +{ + typename Container::value_type result = std::numeric_limits::max(); + for ( auto val : c ) { + if ( val < result ) { + result = val; + } + } + return result; +} + } // namespace algorithm } // namespace tribol -#endif /* SRC_UTILS_ALGORITHM_HPP_ */ \ No newline at end of file +#endif /* SRC_TRIBOL_UTILS_ALGORITHM_HPP_ */ diff --git a/src/tribol/utils/ContactPlaneOutput.hpp b/src/tribol/utils/ContactPlaneOutput.hpp index 4e5ba1a1..18982d33 100644 --- a/src/tribol/utils/ContactPlaneOutput.hpp +++ b/src/tribol/utils/ContactPlaneOutput.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_UTILS_CONTACTPLANEOUTPUT_HPP_ -#define SRC_UTILS_CONTACTPLANEOUTPUT_HPP_ +#ifndef SRC_TRIBOL_UTILS_CONTACTPLANEOUTPUT_HPP_ +#define SRC_TRIBOL_UTILS_CONTACTPLANEOUTPUT_HPP_ #include "tribol/common/BasicTypes.hpp" #include "tribol/common/Parameters.hpp" @@ -37,4 +37,4 @@ void WriteContactPlaneMeshToVtk( const std::string& dir, const VisType v_type, c } // namespace tribol -#endif /* SRC_UTILS_CONTACTPLANEOUTPUT_HPP_ */ +#endif /* SRC_TRIBOL_UTILS_CONTACTPLANEOUTPUT_HPP_ */ diff --git a/src/tribol/utils/DataManager.hpp b/src/tribol/utils/DataManager.hpp index b6e03397..095c5270 100644 --- a/src/tribol/utils/DataManager.hpp +++ b/src/tribol/utils/DataManager.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_UTILS_DATAMANAGER_HPP_ -#define SRC_UTILS_DATAMANAGER_HPP_ +#ifndef SRC_TRIBOL_UTILS_DATAMANAGER_HPP_ +#define SRC_TRIBOL_UTILS_DATAMANAGER_HPP_ // C/C++ includes #include @@ -138,4 +138,4 @@ class DataManager { } // end namespace tribol -#endif /* SRC_UTILS_DATAMANAGER_HPP_ */ +#endif /* SRC_TRIBOL_UTILS_DATAMANAGER_HPP_ */ diff --git a/src/tribol/utils/Math.hpp b/src/tribol/utils/Math.hpp index a00606a1..ffb0d799 100644 --- a/src/tribol/utils/Math.hpp +++ b/src/tribol/utils/Math.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_UTILS_MATH_HPP_ -#define SRC_UTILS_MATH_HPP_ +#ifndef SRC_TRIBOL_UTILS_MATH_HPP_ +#define SRC_TRIBOL_UTILS_MATH_HPP_ // C++ includes #include @@ -143,4 +143,4 @@ TRIBOL_HOST_DEVICE inline void initBoolArray( bool* arr, int length, bool init_v } // namespace tribol -#endif /* SRC_UTILS_MATH_HPP_ */ +#endif /* SRC_TRIBOL_UTILS_MATH_HPP_ */ diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 9b6eff19..110fa74d 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -792,20 +792,28 @@ void TestMesh::allocateAndSetVelocities( IndexT mesh_id, RealT valX, RealT valY, "TestMesh::allocateAndSetVelocities(): please set unique " << "mortarMeshId and nonmortarMeshId prior to calling this routine." ); - // check to see if pointers have been set +// check to see if pointers have been set +#ifdef TRIBOL_DEBUG bool deleteVels = false; +#endif if ( mesh_id == this->mortarMeshId ) { if ( this->vx1 != nullptr ) { delete[] this->vx1; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } if ( this->vy1 != nullptr ) { delete[] this->vy1; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } if ( this->vz1 != nullptr ) { delete[] this->vz1; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } allocRealArray( &this->vx1, this->numTotalNodes, valX ); @@ -816,15 +824,21 @@ void TestMesh::allocateAndSetVelocities( IndexT mesh_id, RealT valX, RealT valY, } else if ( mesh_id == this->nonmortarMeshId ) { if ( this->vx2 != nullptr ) { delete[] this->vx2; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } if ( this->vy2 != nullptr ) { delete[] this->vy2; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } if ( this->vz2 != nullptr ) { delete[] this->vz2; +#ifdef TRIBOL_DEBUG deleteVels = true; +#endif } allocRealArray( &this->vx2, this->numTotalNodes, valX ); @@ -852,18 +866,24 @@ void TestMesh::allocateAndSetBulkModulus( IndexT mesh_id, RealT val ) << "mortarMeshId and nonmortarMeshId prior to calling this routine." ); // check to see if pointers have been set +#ifdef TRIBOL_DEBUG bool deleteData = false; +#endif if ( mesh_id == this->mortarMeshId ) { if ( this->mortar_bulk_mod != nullptr ) { delete[] this->mortar_bulk_mod; +#ifdef TRIBOL_DEBUG deleteData = true; +#endif } allocRealArray( &this->mortar_bulk_mod, this->numMortarFaces, val ); } else if ( mesh_id == this->nonmortarMeshId ) { if ( this->nonmortar_bulk_mod != nullptr ) { delete[] this->nonmortar_bulk_mod; +#ifdef TRIBOL_DEBUG deleteData = true; +#endif } allocRealArray( &this->nonmortar_bulk_mod, this->numNonmortarFaces, val ); @@ -880,18 +900,24 @@ void TestMesh::allocateAndSetBulkModulus( IndexT mesh_id, RealT val ) void TestMesh::allocateAndSetElementThickness( IndexT mesh_id, RealT t ) { // check to see if pointers have been set +#ifdef TRIBOL_DEBUG bool deleteData = false; +#endif if ( mesh_id == this->mortarMeshId ) { if ( this->mortar_element_thickness != nullptr ) { delete[] this->mortar_element_thickness; +#ifdef TRIBOL_DEBUG deleteData = true; +#endif } allocRealArray( &this->mortar_element_thickness, this->numMortarFaces, t ); } else if ( mesh_id == this->nonmortarMeshId ) { if ( this->nonmortar_element_thickness != nullptr ) { delete[] this->nonmortar_element_thickness; +#ifdef TRIBOL_DEBUG deleteData = true; +#endif } allocRealArray( &this->nonmortar_element_thickness, this->numNonmortarFaces, t ); @@ -1359,7 +1385,7 @@ void TestMesh::setupMfemMesh( bool fix_orientation ) // add mortar elements and vertices. Not sure if order of adding // elements matters, but adding vertices should probably correspond // to the global contiguous id system - int mConn[this->numNodesPerElement]; + Array1D mConn( this->numNodesPerElement ); for ( int iel = 0; iel < this->numMortarElements; ++iel ) { for ( int idx = 0; idx < this->numNodesPerElement; ++idx ) { int index = iel * this->numNodesPerElement + idx; @@ -1367,11 +1393,11 @@ void TestMesh::setupMfemMesh( bool fix_orientation ) } switch ( this->cellType ) { case LINEAR_TRIANGLE: { - this->mfem_mesh->AddTet( &mConn[0] ); + this->mfem_mesh->AddTet( mConn.data() ); break; } case LINEAR_QUAD: { - this->mfem_mesh->AddHex( &mConn[0] ); + this->mfem_mesh->AddHex( mConn.data() ); break; } default: { @@ -1386,13 +1412,13 @@ void TestMesh::setupMfemMesh( bool fix_orientation ) vert[1] = this->y[i]; vert[2] = this->z[i]; - this->mfem_mesh->AddVertex( &vert[0] ); + this->mfem_mesh->AddVertex( vert ); } // add nonmortar elements and vertices. Not sure if order of adding // elements matters, but adding vertices should probably correspond // to the global contiguous id system - int sConn[this->numNodesPerElement]; + Array1D sConn( this->numNodesPerElement ); for ( int iel = 0; iel < this->numNonmortarElements; ++iel ) { for ( int idx = 0; idx < this->numNodesPerElement; ++idx ) { int index = iel * this->numNodesPerElement + idx; diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 3255ad54..52272275 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -3,8 +3,8 @@ // // SPDX-License-Identifier: (MIT) -#ifndef SRC_UTILS_TESTUTILS_HPP_ -#define SRC_UTILS_TESTUTILS_HPP_ +#ifndef SRC_TRIBOL_UTILS_TESTUTILS_HPP_ +#define SRC_TRIBOL_UTILS_TESTUTILS_HPP_ #include "tribol/common/Parameters.hpp" @@ -562,4 +562,4 @@ class ExplicitMechanics : public mfem::SecondOrderTimeDependentOperator { } // namespace mfem_ext -#endif /* SRC_UTILS_TESTUTILS_HPP_ */ +#endif /* SRC_TRIBOL_UTILS_TESTUTILS_HPP_ */