diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index b380472a..5cad930c 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -67,7 +67,7 @@ void crackBranchingExample( const std::string filename ) // ==================================================== // Force model // ==================================================== - using model_type = CabanaPD::ForceModel; + using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); // ==================================================== diff --git a/examples/mechanics/elastic_wave.cpp b/examples/mechanics/elastic_wave.cpp index 79e7921f..d2dc19cf 100644 --- a/examples/mechanics/elastic_wave.cpp +++ b/examples/mechanics/elastic_wave.cpp @@ -55,7 +55,8 @@ void elasticWaveExample( const std::string filename ) // Force model // ==================================================== using model_type = - CabanaPD::ForceModel; + CabanaPD::ForceModel; model_type force_model( delta, K, G ); // ==================================================== @@ -100,7 +101,7 @@ void elasticWaveExample( const std::string filename ) // ==================================================== // Create solver // ==================================================== - auto cabana_pd = CabanaPD::createSolverElastic( + auto cabana_pd = CabanaPD::createSolverNoFracture( inputs, particles, force_model ); // ==================================================== diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index dfd270e1..d806d906 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -60,7 +60,7 @@ void fragmentingCylinderExample( const std::string filename ) // ==================================================== // Force model // ==================================================== - using model_type = CabanaPD::ForceModel; + using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); // using model_type = // CabanaPD::ForceModel; diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index b6a3372b..aed4f4e6 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -78,10 +78,10 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== // Force model // ==================================================== - using model_type = CabanaPD::ForceModel; + using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); // using model_type = - // CabanaPD::ForceModel; + // CabanaPD::ForceModel; // model_type force_model( delta, K, G, G0 ); // ==================================================== diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 74491ed3..f5a4e68a 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -81,13 +81,14 @@ void thermalDeformationExample( const std::string filename ) // ==================================================== // Force model // ==================================================== - auto force_model = CabanaPD::createForceModel( - model_type{}, CabanaPD::Elastic{}, *particles, delta, K, alpha, temp0 ); + auto force_model = + CabanaPD::createForceModel( model_type{}, CabanaPD::NoFracture{}, + *particles, delta, K, alpha, temp0 ); // ==================================================== // Create solver // ==================================================== - auto cabana_pd = CabanaPD::createSolverElastic( + auto cabana_pd = CabanaPD::createSolverNoFracture( inputs, particles, force_model ); // ==================================================== diff --git a/examples/thermomechanics/thermal_deformation_heat_transfer.cpp b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp index 1b0fdedc..4b2ddbbf 100644 --- a/examples/thermomechanics/thermal_deformation_heat_transfer.cpp +++ b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp @@ -91,13 +91,13 @@ void thermalDeformationHeatTransferExample( const std::string filename ) // Force model // ==================================================== auto force_model = CabanaPD::createForceModel( - model_type{}, CabanaPD::Elastic{}, *particles, delta, K, kappa, cp, + model_type{}, CabanaPD::NoFracture{}, *particles, delta, K, kappa, cp, alpha, temp0 ); // ==================================================== // Create solver // ==================================================== - auto cabana_pd = CabanaPD::createSolverElastic( + auto cabana_pd = CabanaPD::createSolverNoFracture( inputs, particles, force_model ); // ==================================================== diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index ea63c8fe..26f600a8 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -111,7 +111,8 @@ struct BaseDynamicTemperatureModel } }; -template struct ForceModel; diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp index 9b1e7ad9..f34b0bf9 100644 --- a/src/CabanaPD_HeatTransfer.hpp +++ b/src/CabanaPD_HeatTransfer.hpp @@ -36,7 +36,7 @@ class HeatTransfer : public Force using base_type::_neigh_list; using model_type = ModelType; static_assert( - std::is_same_v ); + std::is_same_v ); // Running with mechanics as well; no reason to rebuild neighbors. template diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index af642f53..de64e1d8 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -93,7 +93,7 @@ class SolverBase template -class SolverElastic +class SolverNoFracture { public: using memory_space = MemorySpace; @@ -114,9 +114,9 @@ class SolverElastic using contact_type = Force; using contact_model_type = ContactModelType; - SolverElastic( input_type _inputs, - std::shared_ptr _particles, - force_model_type force_model ) + SolverNoFracture( input_type _inputs, + std::shared_ptr _particles, + force_model_type force_model ) : inputs( _inputs ) , particles( _particles ) , _init_time( 0.0 ) @@ -124,10 +124,10 @@ class SolverElastic setup( force_model ); } - SolverElastic( input_type _inputs, - std::shared_ptr _particles, - force_model_type force_model, - contact_model_type contact_model ) + SolverNoFracture( input_type _inputs, + std::shared_ptr _particles, + force_model_type force_model, + contact_model_type contact_model ) : inputs( _inputs ) , particles( _particles ) , _init_time( 0.0 ) @@ -488,12 +488,12 @@ class SolverElastic template class SolverFracture - : public SolverElastic + : public SolverNoFracture { public: - using base_type = SolverElastic; + using base_type = SolverNoFracture; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; @@ -735,24 +735,26 @@ class SolverFracture template -auto createSolverElastic( InputsType inputs, - std::shared_ptr particles, - ForceModelType model ) +auto createSolverNoFracture( InputsType inputs, + std::shared_ptr particles, + ForceModelType model ) { - return std::make_shared< - SolverElastic>( + return std::make_shared>( inputs, particles, model ); } template -auto createSolverElastic( InputsType inputs, - std::shared_ptr particles, - ForceModelType model, ContactModelType contact_model ) +auto createSolverNoFracture( InputsType inputs, + std::shared_ptr particles, + ForceModelType model, + ContactModelType contact_model ) { - return std::make_shared>( - inputs, particles, model, contact_model ); + return std::make_shared< + SolverNoFracture>( inputs, particles, model, + contact_model ); } template -struct ForceModel : public BaseForceModel +struct ForceModel : public BaseForceModel { using base_type = BaseForceModel; using base_model = LPS; - using fracture_type = Elastic; + using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; using base_type::delta; @@ -34,7 +34,6 @@ struct ForceModel : public BaseForceModel double theta_coeff; double s_coeff; - ForceModel(){}; ForceModel( const double _delta, const double _K, const double _G, const int _influence = 0 ) : base_type( _delta ) @@ -63,9 +62,10 @@ struct ForceModel : public BaseForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; using thermal_type = base_type::thermal_type; @@ -80,7 +80,6 @@ struct ForceModel : public ForceModel double s0; double bond_break_coeff; - ForceModel() {} ForceModel( const double _delta, const double _K, const double _G, const double _G0, const int _influence = 0 ) : base_type( _delta, _K, _G, _influence ) @@ -106,9 +105,10 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; @@ -124,9 +124,10 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index f0200ec3..2bdc2774 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -19,11 +19,12 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel + : public BaseForceModel { using base_type = BaseForceModel; using base_model = PMB; - using fracture_type = Elastic; + using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; using base_type::delta; @@ -31,7 +32,6 @@ struct ForceModel : public BaseForceModel double c; double K; - ForceModel(){}; ForceModel( const double delta, const double K ) : base_type( delta ) { @@ -54,10 +54,10 @@ struct ForceModel : public BaseForceModel }; template <> -struct ForceModel - : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; using thermal_type = base_type::thermal_type; @@ -69,7 +69,6 @@ struct ForceModel double s0; double bond_break_coeff; - ForceModel() {} ForceModel( const double delta, const double K, const double G0 ) : base_type( delta, K ) { @@ -101,10 +100,10 @@ struct ForceModel }; template <> -struct ForceModel - : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; @@ -117,10 +116,10 @@ struct ForceModel }; template <> -struct ForceModel - : public ForceModel +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; @@ -137,14 +136,16 @@ struct ForceModel }; template -struct ForceModel - : public ForceModel, +struct ForceModel + : public ForceModel, BaseTemperatureModel { - using base_type = ForceModel; + using base_type = + ForceModel; using base_temperature_type = BaseTemperatureModel; using base_model = PMB; - using fracture_type = Elastic; + using fracture_type = NoFracture; using thermal_type = TemperatureDependent; using base_type::c; @@ -158,7 +159,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_temperature_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const TemperatureType _temp, const double _alpha, const double _temp0 = 0.0 ) @@ -176,22 +176,33 @@ struct ForceModel } }; +// Default to Fracture. template -auto createForceModel( PMB, Elastic, ParticleType particles, const double delta, +auto createForceModel( PMB model, ParticleType particles, const double delta, const double K, const double alpha, const double temp0 ) +{ + return createForceModel( model, Fracture{}, particles, delta, K, alpha, + temp0 ); +} + +template +auto createForceModel( PMB, NoFracture, ParticleType particles, + const double delta, const double K, const double alpha, + const double temp0 ) { auto temp = particles.sliceTemperature(); using temp_type = decltype( temp ); - return ForceModel( - delta, K, temp, alpha, temp0 ); + return ForceModel( delta, K, temp, alpha, temp0 ); } template -struct ForceModel - : public ForceModel, +struct ForceModel + : public ForceModel, BaseTemperatureModel { - using base_type = ForceModel; + using base_type = + ForceModel; using base_temperature_type = BaseTemperatureModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; @@ -213,7 +224,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_temperature_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const double _G0, const TemperatureType _temp, const double _alpha, const double _temp0 = 0.0 ) @@ -248,20 +258,21 @@ auto createForceModel( PMB, Fracture, ParticleType particles, { auto temp = particles.sliceTemperature(); using temp_type = decltype( temp ); - return ForceModel( + return ForceModel( delta, K, G0, temp, alpha, temp0 ); } template -struct ForceModel - : public ForceModel, +struct ForceModel + : public ForceModel, BaseDynamicTemperatureModel { - using base_type = - ForceModel; + using base_type = ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; using base_model = PMB; - using fracture_type = Elastic; + using fracture_type = NoFracture; using thermal_type = DynamicTemperature; using base_type::c; @@ -279,7 +290,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const TemperatureType _temp, const double _kappa, const double _cp, const double _alpha, @@ -304,14 +314,14 @@ struct ForceModel }; template -auto createForceModel( PMB, Elastic, ParticleType particles, const double delta, - const double K, const double kappa, const double cp, - const double alpha, const double temp0, +auto createForceModel( PMB, NoFracture, ParticleType particles, + const double delta, const double K, const double kappa, + const double cp, const double alpha, const double temp0, const bool constant_microconductivity = true ) { auto temp = particles.sliceTemperature(); using temp_type = decltype( temp ); - return ForceModel( + return ForceModel( delta, K, temp, kappa, cp, alpha, temp0, constant_microconductivity ); } @@ -346,7 +356,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const double _G0, const TemperatureType _temp, const double _kappa, const double _cp, const double _alpha, @@ -379,7 +388,7 @@ auto createForceModel( PMB, Fracture, ParticleType particles, { auto temp = particles.sliceTemperature(); using temp_type = decltype( temp ); - return ForceModel( + return ForceModel( delta, K, G0, temp, kappa, cp, alpha, temp0, constant_microconductivity ); } diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp index 90916153..8262e625 100644 --- a/src/force/CabanaPD_Force_LPS.hpp +++ b/src/force/CabanaPD_Force_LPS.hpp @@ -70,13 +70,14 @@ namespace CabanaPD { template -class Force> +class Force> : public Force { protected: using base_type = Force; using base_type::_half_neigh; - ForceModel _model; + using model_type = ForceModel; + model_type _model; using base_type::_energy_timer; using base_type::_timer; @@ -89,7 +90,7 @@ class Force> template Force( const bool half_neigh, ParticleType& particles, - const ForceModel model ) + const model_type model ) : base_type( half_neigh, model.delta, particles ) , _model( model ) { @@ -259,13 +260,14 @@ class Force> }; template -class Force> - : public Force> +class Force> + : public Force> { protected: - using base_type = Force>; + using base_type = Force>; using base_type::_half_neigh; - ForceModel _model; + using model_type = ForceModel; + model_type _model; using base_type::_energy_timer; using base_type::_timer; @@ -278,7 +280,7 @@ class Force> template Force( const bool half_neigh, const ParticleType& particles, - const ForceModel model ) + const model_type model ) : base_type( half_neigh, particles, model ) , _model( model ) { @@ -511,12 +513,13 @@ class Force> }; template -class Force> - : public Force> +class Force> + : public Force> { protected: - using base_type = Force>; - ForceModel _model; + using base_type = Force>; + using model_type = ForceModel; + model_type _model; using base_type::_energy_timer; using base_type::_timer; @@ -529,7 +532,7 @@ class Force> template Force( const bool half_neigh, ParticleType& particles, - const ForceModel model ) + const model_type model ) : base_type( half_neigh, particles, model ) , _model( model ) { diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 5d6cd2a4..8a099639 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -71,13 +71,13 @@ namespace CabanaPD { template -class Force> +class Force> : public Force { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; + using model_type = ForceModel; using base_type = Force; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; @@ -181,13 +181,13 @@ class Force> }; template -class Force> +class Force> : public Force { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; + using model_type = ForceModel; using base_type = Force; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; @@ -328,13 +328,15 @@ class Force> }; template -class Force> +class Force> : public Force { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; + using model_type = + ForceModel; using base_type = Force; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; diff --git a/src/force/CabanaPD_HertzianContact.hpp b/src/force/CabanaPD_HertzianContact.hpp index 56435da7..0f6f1f9f 100644 --- a/src/force/CabanaPD_HertzianContact.hpp +++ b/src/force/CabanaPD_HertzianContact.hpp @@ -24,7 +24,7 @@ struct HertzianModel : public ContactModel { // FIXME: This is for use as the primary force model. using base_model = PMB; - using fracture_type = Elastic; + using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; using ContactModel::Rc; // Contact horizon (should be > 2*radius) diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 0d57e5ea..036329e2 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -36,6 +36,8 @@ #include #include +#include + namespace Test { struct LinearTag @@ -54,10 +56,12 @@ struct QuadraticTag //---------------------------------------------------------------------------// // Get the PMB strain energy density (at the center point). // Simplified here because the stretch is constant. -template +template double computeReferenceStrainEnergyDensity( - LinearTag, CabanaPD::ForceModel model, - const int m, const double s0, const double ) + LinearTag, ModelType model, const int m, const double s0, const double, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double W = 0.0; double dx = model.delta / m; @@ -79,10 +83,13 @@ double computeReferenceStrainEnergyDensity( return W; } -template +template double computeReferenceStrainEnergyDensity( - QuadraticTag, CabanaPD::ForceModel model, - const int m, const double u11, const double x ) + QuadraticTag, ModelType model, const int m, const double u11, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double W = 0.0; double dx = model.delta / m; @@ -119,11 +126,13 @@ double computeReferenceForceX( LinearTag, ModelType, const int, const double, // Get the PMB force (at one point). // Assumes zero y/z displacement components. -template -double -computeReferenceForceX( QuadraticTag, - CabanaPD::ForceModel model, - const int m, const double u11, const double x ) +template +double computeReferenceForceX( + QuadraticTag, ModelType model, const int m, const double u11, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double fx = 0.0; double dx = model.delta / m; @@ -247,10 +256,12 @@ double computeReferenceNeighbors( const double delta, const int m ) } // Get the LPS strain energy density (at one point). -template +template double computeReferenceStrainEnergyDensity( - LinearTag, CabanaPD::ForceModel model, - const int m, const double s0, const double ) + LinearTag, ModelType model, const int m, const double s0, const double, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double W = 0.0; double dx = model.delta / m; @@ -283,10 +294,13 @@ double computeReferenceStrainEnergyDensity( return W; } -template +template double computeReferenceStrainEnergyDensity( - QuadraticTag, CabanaPD::ForceModel model, - const int m, const double u11, const double x ) + QuadraticTag, ModelType model, const int m, const double u11, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double W = 0.0; double dx = model.delta / m; @@ -330,11 +344,13 @@ double computeReferenceStrainEnergyDensity( // Get the LPS strain energy density (at one point). // Assumes zero y/z displacement components. -template -double -computeReferenceForceX( QuadraticTag, - CabanaPD::ForceModel model, - const int m, const double u11, const double x ) +template +double computeReferenceForceX( + QuadraticTag, ModelType model, const int m, const double u11, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double fx = 0.0; double dx = model.delta / m; @@ -532,10 +548,12 @@ void checkParticle( QuadraticTag tag, ModelType model, const double s0, checkAnalyticalStrainEnergy( tag, model, s0, W, x ); } -template +template void checkAnalyticalStrainEnergy( - LinearTag, CabanaPD::ForceModel model, - const double s0, const double W, const double ) + LinearTag, ModelType model, const double s0, const double W, const double, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { // Relatively large error for small m. double threshold = W * 0.15; @@ -543,20 +561,25 @@ void checkAnalyticalStrainEnergy( EXPECT_NEAR( W, analytical_W, threshold ); } -template +template void checkAnalyticalStrainEnergy( - LinearTag, CabanaPD::ForceModel model, - const double s0, const double W, const double ) + LinearTag, ModelType model, const double s0, const double W, const double, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { // LPS is exact. double analytical_W = 9.0 / 2.0 * model.K * s0 * s0; EXPECT_FLOAT_EQ( W, analytical_W ); } -template +template void checkAnalyticalStrainEnergy( - QuadraticTag, CabanaPD::ForceModel model, - const double u11, const double W, const double x ) + QuadraticTag, ModelType model, const double u11, const double W, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double threshold = W * 0.05; double analytical_W = @@ -565,10 +588,13 @@ void checkAnalyticalStrainEnergy( EXPECT_NEAR( W, analytical_W, threshold ); } -template +template void checkAnalyticalStrainEnergy( - QuadraticTag, CabanaPD::ForceModel model, - const double u11, const double W, const double x ) + QuadraticTag, ModelType model, const double u11, const double W, + const double x, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double threshold = W * 0.20; double analytical_W = @@ -578,36 +604,46 @@ void checkAnalyticalStrainEnergy( EXPECT_NEAR( W, analytical_W, threshold ); } -template +template void checkAnalyticalForce( - QuadraticTag, CabanaPD::ForceModel model, - const double s0, const double fx ) + QuadraticTag, ModelType model, const double s0, const double fx, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double threshold = fx * 0.10; double analytical_f = 18.0 / 5.0 * model.K * s0; EXPECT_NEAR( fx, analytical_f, threshold ); } -template +template void checkAnalyticalForce( - QuadraticTag, CabanaPD::ForceModel model, - const double s0, const double fx ) + QuadraticTag, ModelType model, const double s0, const double fx, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { double threshold = fx * 0.10; double analytical_f = 2.0 * ( model.K + 4.0 / 3.0 * model.G ) * s0; EXPECT_NEAR( fx, analytical_f, threshold ); } -template -void checkAnalyticalDilatation( CabanaPD::ForceModel, - LinearTag, const double, const double theta ) +template +void checkAnalyticalDilatation( + ModelType, LinearTag, const double, const double theta, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { EXPECT_FLOAT_EQ( 0.0, theta ); } -template -void checkAnalyticalDilatation( CabanaPD::ForceModel, - LinearTag, const double s0, const double theta ) +template +void checkAnalyticalDilatation( + ModelType, LinearTag, const double s0, const double theta, + typename std::enable_if< + ( std::is_same::value ), + int>::type* = 0 ) { EXPECT_FLOAT_EQ( 3 * s0, theta ); } @@ -618,15 +654,8 @@ void checkAnalyticalDilatation( ModelType, QuadraticTag, const double, { } -struct DamageTag -{ -}; -struct NoDamageTag -{ -}; - template -double computeEnergyAndForce( NoDamageTag, ForceType force, +double computeEnergyAndForce( CabanaPD::NoFracture, ForceType force, ParticleType& particles, const int ) { computeForce( force, particles, Cabana::SerialOpTag() ); @@ -634,7 +663,7 @@ double computeEnergyAndForce( NoDamageTag, ForceType force, return Phi; } template -double computeEnergyAndForce( DamageTag, ForceType force, +double computeEnergyAndForce( CabanaPD::Fracture, ForceType force, ParticleType& particles, const int max_neighbors ) { Kokkos::View mu( @@ -647,15 +676,22 @@ double computeEnergyAndForce( DamageTag, ForceType force, } template -void initializeForce( ModelType, ForceType& force, ParticleType& particles ) +void initializeForce( + ModelType, ForceType& force, ParticleType& particles, + typename std::enable_if<( std::is_same::value ), + int>::type* = 0 ) { force.computeWeightedVolume( particles, Cabana::SerialOpTag() ); force.computeDilatation( particles, Cabana::SerialOpTag() ); } -template -void initializeForce( CabanaPD::ForceModel, - ForceType& force, ParticleType& particles ) +template +void initializeForce( + ModelType, ForceType& force, ParticleType& particles, + typename std::enable_if<( std::is_same::value ), + int>::type* = 0 ) { auto max_neighbors = force.getMaxLocalNeighbors(); Kokkos::View mu( "broken_bonds", particles.numLocal(), @@ -683,10 +719,10 @@ void copyTheta( CabanaPD::LPS, ParticleType particles, AoSoAType aosoa_host ) //---------------------------------------------------------------------------// // Main test function. //---------------------------------------------------------------------------// -template -void testForce( ModelType model, const DamageType damage_tag, const double dx, - const double m, const double boundary_width, - const TestType test_tag, const double s0 ) +template +void testForce( ModelType model, const double dx, const double m, + const double boundary_width, const TestType test_tag, + const double s0 ) { auto particles = createParticles( model, test_tag, dx, s0 ); @@ -705,8 +741,9 @@ void testForce( ModelType model, const DamageType damage_tag, const double dx, unsigned int max_neighbors; unsigned long long total_neighbors; force.getNeighborStatistics( max_neighbors, total_neighbors ); - double Phi = - computeEnergyAndForce( damage_tag, force, particles, max_neighbors ); + using fracture_type = typename ModelType::fracture_type; + double Phi = computeEnergyAndForce( fracture_type{}, force, particles, + max_neighbors ); // Make a copy of final results on the host std::size_t num_particle = x.size(); @@ -745,9 +782,10 @@ TEST( TEST_CATEGORY, test_force_pmb ) double dx = 2.0 / 11.0; double delta = dx * m; double K = 1.0; - CabanaPD::ForceModel model( delta, K ); - testForce( model, NoDamageTag{}, dx, m, 1.1, LinearTag{}, 0.1 ); - testForce( model, NoDamageTag{}, dx, m, 1.1, QuadraticTag{}, 0.01 ); + CabanaPD::ForceModel + model( delta, K ); + testForce( model, dx, m, 1.1, LinearTag{}, 0.1 ); + testForce( model, dx, m, 1.1, QuadraticTag{}, 0.01 ); } TEST( TEST_CATEGORY, test_force_linear_pmb ) { @@ -755,9 +793,10 @@ TEST( TEST_CATEGORY, test_force_linear_pmb ) double dx = 2.0 / 11.0; double delta = dx * m; double K = 1.0; - CabanaPD::ForceModel model( delta, - K ); - testForce( model, NoDamageTag{}, dx, m, 1.1, LinearTag{}, 0.1 ); + CabanaPD::ForceModel + model( delta, K ); + testForce( model, dx, m, 1.1, LinearTag{}, 0.1 ); } TEST( TEST_CATEGORY, test_force_lps ) { @@ -767,10 +806,10 @@ TEST( TEST_CATEGORY, test_force_lps ) double delta = dx * m; double K = 1.0; double G = 0.5; - CabanaPD::ForceModel model( delta, K, G, - 1 ); - testForce( model, NoDamageTag{}, dx, m, 2.1, LinearTag{}, 0.1 ); - testForce( model, NoDamageTag{}, dx, m, 2.1, QuadraticTag{}, 0.01 ); + CabanaPD::ForceModel + model( delta, K, G, 1 ); + testForce( model, dx, m, 2.1, LinearTag{}, 0.1 ); + testForce( model, dx, m, 2.1, QuadraticTag{}, 0.01 ); } TEST( TEST_CATEGORY, test_force_linear_lps ) { @@ -779,9 +818,10 @@ TEST( TEST_CATEGORY, test_force_linear_lps ) double delta = dx * m; double K = 1.0; double G = 0.5; - CabanaPD::ForceModel model( - delta, K, G, 1 ); - testForce( model, NoDamageTag{}, dx, m, 2.1, LinearTag{}, 0.1 ); + CabanaPD::ForceModel + model( delta, K, G, 1 ); + testForce( model, dx, m, 2.1, LinearTag{}, 0.1 ); } // Tests without damage, but using damage models. @@ -793,10 +833,9 @@ TEST( TEST_CATEGORY, test_force_pmb_damage ) double K = 1.0; // Large value to make sure no bonds break. double G0 = 1000.0; - CabanaPD::ForceModel model( delta, K, - G0 ); - testForce( model, DamageTag{}, dx, m, 1.1, LinearTag{}, 0.1 ); - testForce( model, DamageTag{}, dx, m, 1.1, QuadraticTag{}, 0.01 ); + CabanaPD::ForceModel model( delta, K, G0 ); + testForce( model, dx, m, 1.1, LinearTag{}, 0.1 ); + testForce( model, dx, m, 1.1, QuadraticTag{}, 0.01 ); } TEST( TEST_CATEGORY, test_force_lps_damage ) { @@ -806,9 +845,8 @@ TEST( TEST_CATEGORY, test_force_lps_damage ) double K = 1.0; double G = 0.5; double G0 = 1000.0; - CabanaPD::ForceModel model( delta, K, G, - G0, 1 ); - testForce( model, DamageTag{}, dx, m, 2.1, LinearTag{}, 0.1 ); - testForce( model, DamageTag{}, dx, m, 2.1, QuadraticTag{}, 0.01 ); + CabanaPD::ForceModel model( delta, K, G, G0, 1 ); + testForce( model, dx, m, 2.1, LinearTag{}, 0.1 ); + testForce( model, dx, m, 2.1, QuadraticTag{}, 0.01 ); } } // end namespace Test diff --git a/unit_test/tstHertz.hpp b/unit_test/tstHertz.hpp index f9e08bb3..bee1f506 100644 --- a/unit_test/tstHertz.hpp +++ b/unit_test/tstHertz.hpp @@ -122,7 +122,7 @@ void testHertzianContact( const std::string filename ) // ==================================================== // Simulation run // ==================================================== - auto cabana_pd = CabanaPD::createSolverElastic( + auto cabana_pd = CabanaPD::createSolverNoFracture( inputs, particles, contact_model ); cabana_pd->init(); cabana_pd->run();