From 0134b5226f63e3b10a8ce73125c1863790c8faf2 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Wed, 4 Sep 2024 14:15:25 -0600 Subject: [PATCH] FEAT: Added option for the constrained incrementation solver --- docs/sphinx/changelog.rst | 1 + src/cpp/tardigrade_hydra.cpp | 23 +++++++- src/cpp/tardigrade_hydra.h | 6 ++ src/cpp/tests/test_tardigrade_hydra.cpp | 77 +++++++++++++++++++++++++ 4 files changed, 105 insertions(+), 2 deletions(-) diff --git a/docs/sphinx/changelog.rst b/docs/sphinx/changelog.rst index aa330b6..68589a2 100644 --- a/docs/sphinx/changelog.rst +++ b/docs/sphinx/changelog.rst @@ -14,6 +14,7 @@ New Features - Added a data storage set object that will automatically set an object when the destructor is called (:pull:`143`). By `Nathan Miller`_. - Added a setDataStorage to the macros for data storage definitions (:pull:`144`). By `Nathan Miller`_. - Added a micromorphic Drucker-Prager plasticity model based on constrained optimization (:pull:`166`). By `Nathan Miller`_. +- Constrained active set solver added as option for incrementation (:pull:`167`). By `Nathan Miller`_. Internal Changes ================ diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index f022687..7f78575 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -1864,10 +1864,21 @@ namespace tardigradeHydra{ floatVector X0 = *getUnknownVector( ); - solveNewtonUpdate( deltaX ); - setBaseQuantities( ); + if ( *getUseSQPSolver( ) ){ + + std::fill( deltaX.begin( ), deltaX.end( ), 0 ); + + solveConstrainedQP( deltaX ); + + } + else{ + + solveNewtonUpdate( deltaX ); + + } + updateUnknownVector( X0 + *getLambda( ) * deltaX ); // Refine the estimate if the new point has a higher residual @@ -2346,6 +2357,14 @@ namespace tardigradeHydra{ active_constraints = std::vector< bool >( getNumConstraints( ), false ); + for ( auto c = getConstraints( )->begin( ); c != getConstraints( )->end( ); c++ ){ + + unsigned int index = ( unsigned int )( c - getConstraints( )->begin( ) ); + + active_constraints[ index ] = ( ( *c ) < 0. ); + + } + } void hydraBase::solveConstrainedQP( floatVector &dx, const unsigned int kmax ){ diff --git a/src/cpp/tardigrade_hydra.h b/src/cpp/tardigrade_hydra.h index 297881e..af95b78 100644 --- a/src/cpp/tardigrade_hydra.h +++ b/src/cpp/tardigrade_hydra.h @@ -1403,6 +1403,8 @@ namespace tardigradeHydra{ unsigned int getNumGrad( ){ /*! Get the number of gradient descent steps performed */ return _NUM_GRAD; } + const bool *getUseSQPSolver( ){ return &_useSQPSolver; } + protected: // Setters that the user may need to access but not override @@ -1633,6 +1635,8 @@ namespace tardigradeHydra{ virtual void initializeActiveConstraints( std::vector< bool > &active_constraints ); + void setUseSQPSolver( const unsigned int &value ){ _useSQPSolver = value; } + private: // Friend classes @@ -1756,6 +1760,8 @@ namespace tardigradeHydra{ floatType _lambda = 1; + bool _useSQPSolver = false; + void setFirstConfigurationJacobians( ); void setPreviousFirstConfigurationJacobians( ); diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index 9df2c37..e89524f 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -6841,3 +6841,80 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_solveConstrainedQP, * boost::unit_test::tol BOOST_TEST( ( result + hydra.initialUnknownVector ) == answer, CHECK_PER_ELEMENT ); } + +BOOST_AUTO_TEST_CASE( test_hydraBase_initializeActiveConstraints, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ + + class hydraBaseMock : public tardigradeHydra::hydraBase{ + + public: + + using tardigradeHydra::hydraBase::hydraBase; + + using tardigradeHydra::hydraBase::setResidualClasses; + + void public_initializeActiveConstraints( std::vector< bool > &active_constraints ){ + + initializeActiveConstraints( active_constraints ); + + } + + protected: + + virtual void setConstraints( ) override{ + + auto constraints = get_setDataStorage_constraints( ); + + *constraints.value = { 2, 6, -2, 0.1, 2 }; + + } + + virtual const unsigned int getNumConstraints( ) override{ return 5; } + + }; + + floatType time = 1.1; + + floatType deltaTime = 2.2; + + floatType temperature = 5.3; + + floatType previousTemperature = 23.4; + + floatVector deformationGradient = { 0.39293837, -0.42772133, -0.54629709, + 0.10262954, 0.43893794, -0.15378708, + 0.9615284 , 0.36965948, -0.0381362 }; + + floatVector previousDeformationGradient = { -0.21576496, -0.31364397, 0.45809941, + -0.12285551, -0.88064421, -0.20391149, + 0.47599081, -0.63501654, -0.64909649 }; + + floatVector previousStateVariables = { 0.53155137, 0.53182759, 0.63440096, 0.84943179, 0.72445532, + 0.61102351, 0.72244338, 0.32295891, 0.36178866, 0.22826323, + 0.29371405, 0.63097612, 0.09210494, 0.43370117, 0.43086276, + 0.4936851 , 0.42583029, 0.31226122, 0.42635131, 0.89338916, + 0.94416002, 0.50183668, 0.62395295, 0.1156184 , 0.31728548, + 0.41482621, 0.86630916, 0.25045537, 0.48303426, 0.98555979, + 0.51948512, 0.61289453, 0.12062867, 0.8263408 , 0.60306013, + 0.54506801, 0.34276383, 0.30412079, 0.0, 0.1 }; + + floatVector parameters = { 1, 2, 3, 4, 5 }; + + unsigned int numConfigurations = 4; + + unsigned int numNonLinearSolveStateVariables = 5; + + unsigned int dimension = 3; + + hydraBaseMock hydra( time, deltaTime, temperature, previousTemperature, deformationGradient, previousDeformationGradient, + { }, { }, + previousStateVariables, parameters, numConfigurations, numNonLinearSolveStateVariables, dimension ); + + std::vector< bool > answer = { false, false, true, false, false }; + + std::vector< bool > result; + + hydra.public_initializeActiveConstraints( result ); + + BOOST_TEST( result == answer, CHECK_PER_ELEMENT ); + +}