Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,8 @@ if ("${PROJECT_SOURCE_DIR}" STREQUAL "${CMAKE_SOURCE_DIR}")
endforeach()

set(BLT_REQUIRED_CLANGFORMAT_VERSION "19" CACHE STRING "")
set(ENABLE_ALL_WARNINGS ON CACHE BOOL "")
set(ENABLE_WARNINGS_AS_ERRORS ON CACHE BOOL "")
endif()

include(${BLT_SOURCE_DIR}/SetupBLT.cmake)
Expand Down
1 change: 0 additions & 1 deletion cmake/TribolCompilerFlags.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -19,4 +19,3 @@ string(APPEND CMAKE_EXE_LINKER_FLAGS " -rdynamic")
if(APPLE)
string(APPEND CMAKE_EXE_LINKER_FLAGS " -Wl,-no_warn_duplicate_libraries")
endif()

32 changes: 28 additions & 4 deletions cmake/thirdparty/SetupMFEM.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,21 @@ set(_MFEM_DIR ${MFEM_DIR}) # Save MFEM_DIR as a non-cache variable
find_package(MFEM CONFIG NO_DEFAULT_PATH PATHS "${MFEM_DIR}/lib/cmake/mfem")
# find_package will overwrite MFEM_DIR, so restore it here
set(MFEM_DIR ${_MFEM_DIR} CACHE PATH "" FORCE)
set(_mfem_uses_mpi FALSE)
set(_mfem_uses_openmp FALSE)

if(MFEM_FOUND)
# MFEM was built with CMake so use that config file
message(STATUS "Using MFEM's CMake config file")
set(MFEM_BUILT_WITH_CMAKE TRUE)
# It looks like include directories are not always built into the target
target_include_directories(mfem INTERFACE ${MFEM_INCLUDE_DIRS})
if(MFEM_USE_MPI)
set(_mfem_uses_mpi TRUE)
endif()
if(MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
set(_mfem_uses_openmp TRUE)
endif()
else()
set(MFEM_BUILT_WITH_CMAKE FALSE)
find_path(
Expand Down Expand Up @@ -68,6 +76,12 @@ else()

# read config.mk file
file(READ "${MFEM_CFG_DIR}/config.mk" mfem_cfg_file_txt)
if(mfem_cfg_file_txt MATCHES "MFEM_USE_MPI[ \\t]*\\+?=[ \\t]*YES")
set(_mfem_uses_mpi TRUE)
endif()
if(mfem_cfg_file_txt MATCHES "MFEM_USE_(LEGACY_)?OPENMP[ \\t]*\\+?=[ \\t]*YES")
set(_mfem_uses_openmp TRUE)
endif()

# parse include flags
string(REGEX MATCHALL "MFEM_TPLFLAGS [^\n]+\n" mfem_tpl_inc_flags ${mfem_cfg_file_txt})
Expand Down Expand Up @@ -128,16 +142,26 @@ else()
TREAT_INCLUDES_AS_SYSTEM ON
EXPORTABLE ON)

# Tribol edit
if(TRIBOL_USE_MPI)
endif()

if(_mfem_uses_mpi)
if(NOT TARGET blt::mpi)
message(FATAL_ERROR "MFEM was built with MPI support, but MPI is not enabled in BLT. Configure with ENABLE_MPI=ON.")
endif()
if(NOT MFEM_BUILT_WITH_CMAKE)
# Note: -lmpifort is being added to MFEM's link line w/o a -L<mpi lib dir>
list(GET MPI_C_LIBRARIES 0 _first_mpi_lib)
get_filename_component(_mpi_lib_dir ${_first_mpi_lib} DIRECTORY)
target_link_directories(mfem INTERFACE ${_mpi_lib_dir})
target_link_libraries(mfem INTERFACE mpi)
endif()
# End Tribol edit
target_link_libraries(mfem INTERFACE blt::mpi)
endif()

if(_mfem_uses_openmp)
if(NOT TARGET blt::openmp)
message(FATAL_ERROR "MFEM was built with OpenMP support, but OpenMP is not enabled in BLT. Configure with ENABLE_OPENMP=ON.")
endif()
target_link_libraries(mfem INTERFACE blt::openmp)
endif()

include(FindPackageHandleStandardArgs)
Expand Down
3 changes: 0 additions & 3 deletions src/tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -251,9 +251,6 @@ endif()
if( MFEM_FOUND )

set(mfem_smoke_depends mfem gtest)
# MFEM may require OpenMP even when Tribol does not, and FindMFEM.cmake does not propagate that dependency.
blt_list_append(TO mfem_smoke_depends ELEMENTS blt::openmp IF ENABLE_OPENMP )
blt_list_append(TO mfem_smoke_depends ELEMENTS blt::mpi IF TRIBOL_USE_MPI )
blt_add_executable(
NAME mfem_smoke_test
SOURCES mfem_smoke.cpp
Expand Down
4 changes: 3 additions & 1 deletion src/tests/tribol_comp_geom_3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -88,14 +88,16 @@ class CompGeomTest : public ::testing::Test {

tribol::setKinematicConstantPenalty( meshId, 1.0 );

return tribol::update( 1, 1., dt );
const int update_err = tribol::update( 1, 1., dt );

delete[] fx;
delete[] fy;
delete[] fz;
delete[] vx;
delete[] vy;
delete[] vz;

return update_err;
}

protected:
Expand Down
8 changes: 6 additions & 2 deletions src/tribol/geom/CompGeom.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -960,7 +960,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CheckInterfacePairByMethod(
InterfacePair& pair, const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, const Parameters& params,
ContactCase TRIBOL_UNUSED_PARAM( cCase ), bool& isInteracting, CompGeom::Viewer& cg, IndexT* plane_ct )
{
auto dim = static_cast<int>( mesh1.spatialDimension() );
auto dim = mesh1.spatialDimension();
T my_plane( &pair, params, dim );
FaceGeomException face_err = NO_FACE_GEOM_EXCEPTION;
if ( dim == 3 ) {
Expand Down Expand Up @@ -2151,7 +2151,8 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c
cx[2], 1, &xInter_local, &yInter_local );

// get the local coordinates of the other face's centroid
RealT cx_other_local, cy_other_local;
RealT cx_other_local = 0;
RealT cy_other_local = 0;
RealT cz = 0.; // dummy arg.
VertexAvgCentroid( &x_other_local[0], &y_other_local[0], nullptr, num_nodes_other, cx_other_local,
cy_other_local, cz );
Expand Down Expand Up @@ -2465,6 +2466,9 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::projectPointsAndCom
<< "input number of vertices to match number of nodes per element." );
}
}
#else
TRIBOL_UNUSED_VAR( m1 );

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

any reason you dont use [[maybe_unused]]?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This could be a follow up PR across tribol

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

No good reason. I agree -- worth cleaning up in a follow-on PR.

TRIBOL_UNUSED_VAR( m2 );
#endif

constexpr int max_nodes_per_clipped_face = 5;
Expand Down
4 changes: 2 additions & 2 deletions src/tribol/geom/GeomUtilities.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@

#include <float.h>
#include <cmath>
#include <limits>

#include "axom/core/NumericLimits.hpp"
#include "axom/slic.hpp"

#include "tribol/common/Parameters.hpp"
Expand Down Expand Up @@ -1895,7 +1895,7 @@ TRIBOL_HOST_DEVICE inline void PolyInterYCentroid( const int namax, const RealT*
RealT vol;

// calculate origin shift to avoid roundoff errors
RealT realt_max = std::numeric_limits<RealT>::max();
RealT realt_max = axom::numeric_limits<RealT>::max();

// clang-format off
RealT xorg = realt_max;
Expand Down
7 changes: 5 additions & 2 deletions src/tribol/integ/FE.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -472,6 +472,7 @@ TRIBOL_HOST_DEVICE inline void GalerkinEval( const RealT* const x, const RealT p
TRIBOL_HOST_DEVICE inline void EvalBasis( const RealT* const x, const RealT pX, const RealT pY, const RealT pZ,
const int numPoints, const int vertexId, RealT& phi )
{
phi = 0.;
if ( numPoints > 2 ) {
WachspressBasis( x, pX, pY, pZ, numPoints, vertexId, phi );
} else if ( numPoints == 2 ) {
Expand Down Expand Up @@ -524,14 +525,16 @@ TRIBOL_HOST_DEVICE inline void SegmentBasis( const RealT* const x, const RealT p
TRIBOL_HOST_DEVICE inline void WachspressBasis( const RealT* const x, const RealT pX, const RealT pY, const RealT pZ,
const int numPoints, const int vertexId, RealT& phi )
{
constexpr int max_nodes_per_elem = 4;
#ifdef TRIBOL_USE_HOST
SLIC_ERROR_IF( numPoints < 3, "WachspressBasis: numPoints < 3." );
SLIC_ERROR_IF( numPoints > max_nodes_per_elem, "WachspressBasis: numPoints > 4." );
SLIC_ERROR_IF( vertexId < 0 || vertexId >= numPoints, "WachspressBasis: vertexId is out of bounds." );
#endif

// first compute the areas of all the triangles formed by the i-1,i,i+1 vertices.
// These consist of all the numerators in the Wachspress formulation
// NOTE: this limits the routine to 4 noded quadrilaterals
constexpr int max_nodes_per_elem = 4;
RealT triVertArea[max_nodes_per_elem];
for ( int i = 0; i < numPoints; ++i ) {
// determine the i-1, i, i+1 vertices
Expand Down Expand Up @@ -585,7 +588,7 @@ TRIBOL_HOST_DEVICE inline void WachspressBasis( const RealT* const x, const Real

// third, compute all of the weights per Wachspress formulation
RealT weight[max_nodes_per_elem];
RealT myWeight;
RealT myWeight = 0.;
RealT weightSum = 0.;
for ( int i = 0; i < numPoints; ++i ) {
int vId = i;
Expand Down
7 changes: 4 additions & 3 deletions src/tribol/mesh/CouplingScheme.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ class CouplingScheme {
*
* @return spatial dimension
*/
TRIBOL_HOST_DEVICE int spatialDimension() const { return static_cast<int>( m_mesh1.spatialDimension() ); }
TRIBOL_HOST_DEVICE int spatialDimension() const { return m_mesh1.spatialDimension(); }

/**
* @brief Return a view of the first mesh in the coupling scheme
Expand Down Expand Up @@ -1072,7 +1072,7 @@ TRIBOL_HOST_DEVICE inline bool CouplingScheme::Viewer::pruneMethodFacePair( cons

auto& mesh1 = this->getMesh1View();
auto& mesh2 = this->getMesh2View();
int dim = static_cast<int>( mesh1.spatialDimension() );
int dim = mesh1.spatialDimension();
int num_nodes_face_1 = static_cast<int>( mesh1.numberOfNodesPerElement() );
int num_nodes_face_2 = static_cast<int>( mesh2.numberOfNodesPerElement() );

Expand Down Expand Up @@ -1111,7 +1111,8 @@ TRIBOL_HOST_DEVICE inline bool CouplingScheme::Viewer::pruneMethodFacePair( cons
}
}

RealT nrml[max_dim], cx[max_dim];
RealT nrml[max_dim]{};
RealT cx[max_dim]{};

switch ( m_contact_method ) {
case ALIGNED_MORTAR:
Expand Down
2 changes: 1 addition & 1 deletion src/tribol/mesh/MeshData.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ class MeshData {
*
* @return spatial dimension
*/
TRIBOL_HOST_DEVICE IndexT spatialDimension() const { return m_position.size(); }
TRIBOL_HOST_DEVICE int spatialDimension() const { return static_cast<int>( m_position.size() ); }

/**
* @brief Number of nodes in the mesh
Expand Down
17 changes: 11 additions & 6 deletions src/tribol/physics/AlignedMortar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,9 @@ void ComputeAlignedMortarWeights( SurfaceContactElem& elem )
// also initializes the array
elem.allocateMortarWts();

RealT phiNonmortarA, phiNonmortarB, phiMortarA;
RealT phiNonmortarA = 0.;
RealT phiNonmortarB = 0.;
RealT phiMortarA = 0.;

// loop over nodes "a", where node "a" can be a nonmortar node or a mortar node
for ( int a = 0; a < elem.numFaceVert; ++a ) {
Expand Down Expand Up @@ -84,8 +86,7 @@ void ComputeNodalGap<ALIGNED_MORTAR>( SurfaceContactElem& elem )

// set the distance magnitude tolerance as the longest edge of
// the mortar face
RealT magTol;
RealT magTest = 0.;
RealT magTol = 0.;
for ( int k = 0; k < elem.numFaceVert; ++k ) {
int idPlus = ( k == ( elem.numFaceVert - 1 ) ) ? 0 : k + 1;
RealT dx = elem.faceCoords1[elem.dim * idPlus] - elem.faceCoords1[elem.dim * k];
Expand All @@ -94,8 +95,7 @@ void ComputeNodalGap<ALIGNED_MORTAR>( SurfaceContactElem& elem )

RealT mag = magnitude( dx, dy, dz );

magTol = ( mag > magTest ) ? mag : magTest;
magTest = mag;
magTol = ( mag > magTol ) ? mag : magTol;
}

// loop over nodes on nonmortar side
Expand All @@ -113,7 +113,7 @@ void ComputeNodalGap<ALIGNED_MORTAR>( SurfaceContactElem& elem )
// determine which mortar node is aligned with
// nonmortar node "a"
//////////////////////////////////////////////
int mortarNodeId;
int mortarNodeId = -1;
RealT v[3] = { 0., 0., 0. };
RealT magTest = magTol;
// loop over nodes on the mortar side
Expand All @@ -130,6 +130,11 @@ void ComputeNodalGap<ALIGNED_MORTAR>( SurfaceContactElem& elem )
}
}

SLIC_ERROR_IF( mortarNodeId < 0, "ComputeNodalGap< ALIGNED_MORTAR >: unable to find an aligned mortar node." );
if ( mortarNodeId < 0 ) {
continue;
}

// store local gap
v[0] = elem.faceCoords1[elem.dim * mortarNodeId] - elem.faceCoords2[elem.dim * a];
v[1] = elem.faceCoords1[elem.dim * mortarNodeId + 1] - elem.faceCoords2[elem.dim * a + 1];
Expand Down
8 changes: 4 additions & 4 deletions src/tribol/physics/ContactFormulationFactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,18 +18,18 @@ std::unique_ptr<ContactFormulation> createContactFormulation( CouplingScheme* cs
}

if ( cs->getContactMethod() == ENERGY_MORTAR ) {
#if defined( TRIBOL_USE_ENZYME ) && defined( BUILD_REDECOMP )
// Default parameters for now, or extract from CouplingScheme if available
double k = 1000.0;
double delta = 0.1;
int N = 3;
bool enzyme_quadrature = true;

// ENERGY_MORTAR supports a penalty-style mode driven by the kinematic penalty parameters, even if the coupling
// scheme is registered with LM enforcement (which is often done to enable submesh/pressure infrastructure).
const auto& penalty_opts = cs->getEnforcementOptions().penalty_options;
// TODO: Figure out how contact formulations interact with coupling scheme duplication (SRW)
bool use_penalty_ = penalty_opts.kinematic_calc_set;
bool use_penalty = penalty_opts.kinematic_calc_set;

#if defined( TRIBOL_USE_ENZYME ) && defined( BUILD_REDECOMP )
if ( cs->hasMfemData() ) {
// Attempt to get penalty from MfemMeshData if available
auto* k1_ptr = cs->getMfemMeshData()->GetMesh1KinematicConstantPenalty();
Expand All @@ -45,7 +45,7 @@ std::unique_ptr<ContactFormulation> createContactFormulation( CouplingScheme* cs

return std::make_unique<EnergyMortarAdapter>( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(),
*cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature,
use_penalty_ );
use_penalty );
#else
SLIC_ERROR_ROOT( "ENERGY_MORTAR requires Enzyme and redecomp to be built." );
return nullptr;
Expand Down
4 changes: 3 additions & 1 deletion src/tribol/physics/Mortar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,9 @@ void ComputeMortarWeights( SurfaceContactElem& elem )
// also initializes the array
elem.allocateMortarWts();

RealT phiNonmortarA, phiNonmortarB, phiMortarA;
RealT phiNonmortarA = 0.;
RealT phiNonmortarB = 0.;
RealT phiMortarA = 0.;

// loop over number of nodes on the nonmortar or mortar depending on whether forming
// nonmortar/nonmortar or mortar/nonmortar weights
Expand Down
4 changes: 2 additions & 2 deletions src/tribol/search/InterfacePairFinder.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -490,7 +490,7 @@ class BvhSearch : public SearchBase {
{
auto boxes_view = boxes.view();
forAllExec( m_coupling_scheme->getExecutionMode(), mesh.numberOfElements(),
[this, mesh, boxes_view, binning_proximity] TRIBOL_HOST_DEVICE( IndexT i ) {
[mesh, boxes_view, binning_proximity] TRIBOL_HOST_DEVICE( IndexT i ) {
BoxT box;
auto num_nodes_per_elem = mesh.numberOfNodesPerElement();
for ( IndexT j{ 0 }; j < num_nodes_per_elem; ++j ) {
Expand All @@ -516,7 +516,7 @@ class BvhSearch : public SearchBase {
* Expands bounding box by projecting the face normal by a distance
* equal to the effective face radius
*/
TRIBOL_HOST_DEVICE void expandBBoxNormal( BoxT& bbox, const VectorT& faceNormal, const RealT faceRadius )
TRIBOL_HOST_DEVICE static void expandBBoxNormal( BoxT& bbox, const VectorT& faceNormal, const RealT faceRadius )
{
PointT p0 = bbox.getCentroid();
RayT outwardRay( p0, faceNormal );
Expand Down
2 changes: 1 addition & 1 deletion src/tribol/utils/TestUtils.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1767,7 +1767,7 @@ void TestMesh::testMeshToVtk( const std::string& dir, int cycle, RealT time )
// specify integer id for each cell type.
mesh << "CELL_TYPES " << this->numTotalElements << std::endl;

int element_id;
int element_id = -1;
switch ( this->numNodesPerElement ) {
case 8:
element_id = 12; // vtk 8-node hexahedron
Expand Down
Loading