From 2b938a7bcf5a0ffe15e9b7a81d3f194b7d75ffc3 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Fri, 30 Aug 2024 11:08:39 -0600 Subject: [PATCH 01/10] MAINT: Added ability to identify the number of constraint equations --- src/cpp/tardigrade_hydra.h | 9 ++++++++- src/cpp/tests/test_tardigrade_hydra.cpp | 22 +++++++++++++++++++++- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/cpp/tardigrade_hydra.h b/src/cpp/tardigrade_hydra.h index c5f34aa..09958bc 100644 --- a/src/cpp/tardigrade_hydra.h +++ b/src/cpp/tardigrade_hydra.h @@ -718,7 +718,7 @@ namespace tardigradeHydra{ * * \param &r: The residual to be copied */ - residualBase( residualBase &r ) : hydra( r.hydra ), _numEquations( *r.getNumEquations( ) ){ } + residualBase( residualBase &r ) : hydra( r.hydra ), _numEquations( *r.getNumEquations( ) ), _numConstraints( *r.getNumConstraints( ) ){ } hydraBase* hydra; //!< The hydra class which owns the residualBase object @@ -884,6 +884,9 @@ namespace tardigradeHydra{ //! Get the number of equations the residual defined const unsigned int* getNumEquations( ){ return &_numEquations; } + //! Get the number of constraints the residual defined + const unsigned int* getNumConstraints( ){ return &_numConstraints; } + void addIterationData( dataBase *data ); template @@ -935,6 +938,8 @@ namespace tardigradeHydra{ protected: + void setNumConstraints( const unsigned int numConstraints ){ _numConstraints = numConstraints; } + void setPenaltyIndices( const std::vector< unsigned int > &indices ){ _penalty_indices = indices; } void unexpectedError( ){ @@ -989,6 +994,8 @@ namespace tardigradeHydra{ unsigned int _numEquations; //!< The number of residual equations + unsigned int _numConstraints = 0; //!< The number of constraint equations + bool _useProjection = false; //!< Flag for whether to use the projection or not std::vector< unsigned int > _penalty_indices; //!< The indices of the variables which should be penalized for negative values diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index 85a8332..4bd9b76 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -2441,12 +2441,32 @@ BOOST_AUTO_TEST_CASE( test_residualBase_residualBase, * boost::unit_test::tolera unsigned int numEquations = 3; - tardigradeHydra::residualBase residual( &hydra, numEquations ); + unsigned int numConstraints = 5; + + class residualBaseMock : public tardigradeHydra::residualBase{ + + public: + + using tardigradeHydra::residualBase::residualBase; + + void public_setNumConstraints( const unsigned int &val ){ + + setNumConstraints( val ); + + } + + }; + + residualBaseMock residual( &hydra, numEquations ); + + residual.public_setNumConstraints( numConstraints ); BOOST_CHECK( residual.hydra == &hydra ); BOOST_CHECK( *residual.getNumEquations( ) == numEquations ); + BOOST_CHECK( *residual.getNumConstraints( ) == numConstraints ); + } BOOST_AUTO_TEST_CASE( test_residualBase_checkDefaults, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ From b297aad111ff6074590ee4cc1afc501a974aff1f Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Fri, 30 Aug 2024 14:01:27 -0600 Subject: [PATCH 02/10] WIP: Initial commit of an active set Newton-Raphson solver --- src/cpp/tardigrade_hydra.cpp | 322 +++++++++++++++++++++++++++++++++++ src/cpp/tardigrade_hydra.h | 40 ++++- 2 files changed, 357 insertions(+), 5 deletions(-) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index d11c3ef..f8c4188 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2203,6 +2203,328 @@ namespace tardigradeHydra{ } + void hydraBase::assembleKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ){ + /*! + * Assemble the Karush-Kuhn-Tucker matrix for an inequality constrained Newton-Raphson solve + * + * \param &KKTMatrix: The Karush-Kuhn-Tucker matrix + * \param &active_constraints: The vector of currently active constraints. + */ + + const unsigned int numUnknowns = getNumUnknowns( ); + + const unsigned int numConstraints = active_constraints.size( ); + + KKTMatrix = floatVector( ( numUnknowns + numConstraints ) * ( numUnknowns + numConstraints ), 0 ); + + + Eigen::Map< Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > K( KKTMatrix.data( ), ( numUnknowns + numConstraints ), ( numUnknowns + numConstraints ) ); + + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > J( getFlatJacobian( )->data( ), numUnknowns, numUnknowns ); + + K.block( 0, 0, numUnknowns, numUnknowns ) = ( J.transpose( ) * J ).eval( ); + + for ( unsigned int I = 0; I < numUnknowns; I++ ){ + + KKTMatrix[ ( numUnknowns + numConstraints ) * I + I ] = ( *getMuk( ) ); + + } + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( active_constraints[ i ] ){ + + for ( unsigned int I = 0; I < numUnknowns; I++ ){ + + KKTMatrix[ ( numUnknowns + numConstraints ) * ( I ) + numUnknowns + i ] = ( *getConstraintJacobians( ) )[ numUnknowns * i + I ]; + KKTMatrix[ ( numUnknowns + numConstraints ) * ( numUnknowns + i ) + I ] = ( *getConstraintJacobians( ) )[ numUnknowns * i + I ]; + + } + + } + else{ + + KKTMatrix[ ( numUnknowns + numConstraints ) * ( numUnknowns + i ) + numUnknowns + i ] = 1; + + } + + } + + } + + void hydraBase::updateKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ){ + /*! + * Update the KKTMatrix if the active constraints have changed + * + * \param &KKTMatrix: The Karush-Kuhn-Tucker matrix + * \param &active_constraints: The vector of currently active constraints. + */ + + const unsigned int numUnknowns = getNumUnknowns( ); + + const unsigned int numConstraints = active_constraints.size( ); + + Eigen::Map< Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > K( KKTMatrix.data( ), ( numUnknowns + numConstraints ), ( numUnknowns + numConstraints ) ); + + K.block( 0, numUnknowns, numUnknowns, numConstraints ).setZero( ); + + K.block( numUnknowns, 0, numConstraints, numUnknowns ).setZero( ); + + K.block( numUnknowns, numUnknowns, numConstraints, numConstraints ).setZero( ); + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( active_constraints[ i ] ){ + + for ( unsigned int I = 0; I < numUnknowns; I++ ){ + + KKTMatrix[ ( numUnknowns + numConstraints ) * ( I ) + numUnknowns + i ] = ( *getConstraintJacobians( ) )[ numUnknowns * i + I ]; + KKTMatrix[ ( numUnknowns + numConstraints ) * ( numUnknowns + i ) + I ] = ( *getConstraintJacobians( ) )[ numUnknowns * i + I ]; + + } + + } + else{ + + KKTMatrix[ ( numUnknowns + numConstraints ) * ( numUnknowns + i ) + numUnknowns + i ] = 1; + + } + + } + + } + + void hydraBase::assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, std::vector< bool > &active_constraints ){ + /*! + * Assemble the right hand side vector for the KKT matrix + * + * \param &dx: The delta vector being solved for + * \param &KKTRHSVector: The right hand size vector for the KKT matrix + * \param &active_constraints: The active constraint vector + */ + + const unsigned int numUnknowns = getNumUnknowns( ); + + const unsigned int numConstraints = getNumConstraints( ); + + KKTRHSVector = floatVector( numUnknowns + numConstraints, 0 ); + + Eigen::Map< const Eigen::Vector< floatType, -1 > > _dx( dx.data( ), numUnknowns ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > RHS( KKTRHSVector.data( ), ( numUnknowns + numConstraints ), ( numUnknowns + numConstraints ) ); + + Eigen::Map< const Eigen::Vector< floatType, -1 > > R( getResidual( )->data( ), numUnknowns ); + + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > J( getFlatJacobian( )->data( ), numUnknowns, numUnknowns ); + + RHS.head( numUnknowns ) += ( J.transpose( ) * ( R + J * _dx ) + ( *getMuk( ) ) * _dx ).eval( ); + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( active_constraints[ i ] ){ + + KKTRHSVector[ numUnknowns + i ] = ( *getConstraints( ) )[ i ]; + + for ( unsigned int I = 0; I < numUnknowns; I++ ){ + + KKTRHSVector[ numUnknowns + i ] += ( *getConstraintJacobians( ) )[ numUnknowns * i + I ] * dx[ I ]; + + } + + } + + } + + } + + void hydraBase::solveConstrainedQP( floatVector &dx, const unsigned int kmax ){ + /*! + * Solve the constrained QP problem to estimate the desired step size + * + * \param &dx: The change in the unknown vector + * \param kmax: The maximum number of iterations (defaults to 100) + */ + + const unsigned int numUnknowns = getNumUnknowns( ); + + const unsigned int numConstraints = getNumConstraints( ); + + floatVector K; + + floatVector RHS; + + std::vector< bool > active_constraints( numConstraints, false ); + + assembleKKTRHSVector( dx, RHS, active_constraints ); + + assembleKKTMatrix( K, active_constraints ); + + floatType tol = ( *getRelativeTolerance( ) ) * ( tardigradeVectorTools::l2norm( RHS ) ) + ( *getAbsoluteTolerance( ) ); + + unsigned int k = 0; + + floatVector y( numUnknowns + numConstraints, 0 ); + + floatVector ck( numConstraints, 0 ); + + floatVector ctilde( numConstraints, 0 ); + + floatVector negp( numUnknowns, 0 ); + + floatVector lambda( numConstraints, 0 ); + + floatVector P( numUnknowns + numConstraints, 0 ); + + for ( unsigned int i = 0; i < ( numUnknowns + numConstraints ); i++ ){ + + P[ i ] = 1 / std::max( std::fabs( *std::max_element( K.begin( ) + ( numUnknowns + numConstraints ) * i, + K.begin( ) + ( numUnknowns + numConstraints ) * ( i + 1 ), + [ ]( const floatType &a, const floatType &b ){ return std::fabs( a ) < std::fabs( b ); } ) ), 1e-15 ); + + } + + const floatVector *c0 = getConstraints( ); + + const floatVector *A = getConstraintJacobians( ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > _dx( dx.data( ), numUnknowns ); + + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > _K( K.data( ), numConstraints + numUnknowns, numConstraints + numUnknowns ); + + Eigen::Map< const Eigen::Vector< floatType, -1 > > _RHS( RHS.data( ), numConstraints + numUnknowns ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > _y( y.data( ), numConstraints + numUnknowns ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > _ck( ck.data( ), numConstraints ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > _ctilde( ctilde.data( ), numConstraints ); + + Eigen::Map< Eigen::Vector< floatType, -1 > > _negp( negp.data( ), numUnknowns ); + + Eigen::Map< const Eigen::Vector< floatType, -1 > > _c0( c0->data( ), numConstraints ); + + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > _A( A->data( ), numConstraints, numUnknowns ); + + Eigen::Map< const Eigen::Vector< floatType, -1 > > _P( P.data( ), numUnknowns + numConstraints ); + + tardigradeVectorTools::solverType< floatType > linearSolver; + + while ( k < kmax ){ + + linearSolver = tardigradeVectorTools::solverType< floatType >( _P.asDiagonal( ) * _K ); + + _y = linearSolver.solve( _P.asDiagonal( ) * _RHS ); + + negp = floatVector( y.begin( ), y.begin( ) + numUnknowns ); + + lambda = floatVector( y.begin( ) + numUnknowns, y.end( ) ); + + if ( tardigradeVectorTools::l2norm( negp ) <= tol ){ + + bool negLambda = false; + + floatType minLambda = 1; + + unsigned int imin = 0; + + for ( auto v = std::begin( lambda ); v != std::end( lambda ); v++ ){ + + if ( *v < 0 ){ + + negLambda = true; + + if ( ( *v ) < minLambda ){ + + imin = ( unsigned int )( v - std::begin( lambda ) ); + + minLambda = *v; + + } + + } + + } + + if ( negLambda ){ + + active_constraints[ imin ] = false; + + } + else{ + + return; + + } + + } + else{ + + _ck = ( _c0 + _A * _dx ).eval( ); + + _ctilde = ( _c0 + _A * ( _dx - _negp ) ).eval( ); + + floatType alpha = 1.0; + + unsigned int iblock = 0; + + bool newBlock = false; + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( !active_constraints[ i ] ){ + + if ( ctilde[ i ] < -tol ){ + + floatType alpha_trial = -ck[ i ] / ( ctilde[ i ] - ck[ i ] ); + + if ( alpha_trial <= alpha ){ + + iblock = i; + + alpha = alpha_trial; + + newBlock = true; + + } + + } + + } + + } + + if ( newBlock ){ + + active_constraints[ iblock ] = true; + + } + + dx -= alpha * negp; + + updateKKTMatrix( K, active_constraints ); + + assembleKKTRHSVector( dx, RHS, active_constraints ); + + } + + k++; + + } + + } + + void hydraBase::setConstraints( ){ + /*! + * Set the constraint values + */ + } + + void hydraBase::setConstraintJacobians( ){ + /*! + * Set the constraint Jacobians values + */ + } + errorOut dummyMaterialModel( floatVector &stress, floatVector &statev, floatMatrix &ddsdde, floatType &SSE, floatType &SPD, floatType &SCD, floatType &RPL, floatVector &ddsddt, floatVector &drplde, floatType &DRPLDT, const floatVector &strain, const floatVector &dstrain, const floatVector &time, const floatType &DTIME, const floatType &TEMP, diff --git a/src/cpp/tardigrade_hydra.h b/src/cpp/tardigrade_hydra.h index 09958bc..a038279 100644 --- a/src/cpp/tardigrade_hydra.h +++ b/src/cpp/tardigrade_hydra.h @@ -1097,9 +1097,24 @@ namespace tardigradeHydra{ //! Get a reference to the number of state variables involved in the non-linear solve const unsigned int* getNumNonLinearSolveStateVariables( ){ return &_numNonLinearSolveStateVariables; } - //! Get a reference to the number of terms in the unknown vector + //! Get the number of terms in the unknown vector virtual const unsigned int getNumUnknowns( ){ return ( *getNumConfigurations( ) ) * ( *getConfigurationUnknownCount( ) ) + *getNumNonLinearSolveStateVariables( ); } + //! Get the value of the number of constraint equations + virtual const unsigned int getNumConstraints( ){ + + unsigned int value = 0; + + for ( auto v = getResidualClasses( )->begin( ); v != getResidualClasses( )->end( ); v++ ){ + + value += *( *v )->getNumConstraints( ); + + } + + return value; + + } + //! Get a reference to the dimension constexpr unsigned int getDimension( ){ return _dimension; } @@ -1604,6 +1619,17 @@ namespace tardigradeHydra{ virtual tardigradeHydra::hydraBase::setDataStorageIteration get_setDataStorage_stress( ); + virtual void assembleKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ); + + virtual void updateKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ); + + virtual void assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, std::vector< bool > &active_constraints ); + + virtual void solveConstrainedQP( floatVector &x, const unsigned int kmax=100 ); + + virtual void setConstraints( ); + + virtual void setConstraintJacobians( ); private: @@ -1782,13 +1808,17 @@ namespace tardigradeHydra{ TARDIGRADE_HYDRA_DECLARE_ITERATION_STORAGE( private, dResidualNormdX, floatVector, setdResidualNormdX ) - TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, set_dF1dF, get_dF1dF, dF1dF, secondOrderTensor, setFirstConfigurationJacobians ) + TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, setConstraints, getConstraints, constraints, floatVector, setConstraints ) + + TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, setConstraintJacobians, getConstraintJacobians, constraintJacobians, floatVector, setConstraintJacobians ) + + TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, set_dF1dF, get_dF1dF, dF1dF, secondOrderTensor, setFirstConfigurationJacobians ) - TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, set_dF1dFn, get_dF1dFn, dF1dFn, floatVector, setFirstConfigurationJacobians ) + TARDIGRADE_HYDRA_DECLARE_NAMED_ITERATION_STORAGE( private, set_dF1dFn, get_dF1dFn, dF1dFn, floatVector, setFirstConfigurationJacobians ) - TARDIGRADE_HYDRA_DECLARE_NAMED_PREVIOUS_STORAGE( private, set_previousdF1dF, get_previousdF1dF, previousdF1dF, secondOrderTensor, setPreviousFirstConfigurationJacobians ) + TARDIGRADE_HYDRA_DECLARE_NAMED_PREVIOUS_STORAGE( private, set_previousdF1dF, get_previousdF1dF, previousdF1dF, secondOrderTensor, setPreviousFirstConfigurationJacobians ) - TARDIGRADE_HYDRA_DECLARE_NAMED_PREVIOUS_STORAGE( private, set_previousdF1dFn, get_previousdF1dFn, previousdF1dFn, floatVector, setPreviousFirstConfigurationJacobians ) + TARDIGRADE_HYDRA_DECLARE_NAMED_PREVIOUS_STORAGE( private, set_previousdF1dFn, get_previousdF1dFn, previousdF1dFn, floatVector, setPreviousFirstConfigurationJacobians ) }; From 44f98b0ac1c14c319b45be9734f4b1d84c334b50 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Fri, 30 Aug 2024 15:23:57 -0600 Subject: [PATCH 03/10] FEAT: Added construction of the complete constraint vector and constraint Jacobian --- src/cpp/tardigrade_hydra.cpp | 44 +++++++ src/cpp/tests/test_tardigrade_hydra.cpp | 153 ++++++++++++++++++++++++ 2 files changed, 197 insertions(+) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index f8c4188..f6fbaa7 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2517,12 +2517,56 @@ namespace tardigradeHydra{ /*! * Set the constraint values */ + + const unsigned int numConstraints = getNumConstraints( ); + + auto constraints = get_setDataStorage_constraints( ); + + constraints.zero( numConstraints ); + + unsigned int offset = 0; + + for ( auto v = getResidualClasses( )->begin( ); v != getResidualClasses( )->end( ); v++ ){ + + if ( ( *( *v )->getNumConstraints( ) ) > 0 ){ + + std::copy( ( *v )->getConstraints( )->begin( ), ( *v )->getConstraints( )->end( ), constraints.value->begin( ) + offset ); + + offset += ( *v )->getConstraints( )->size( ); + + } + + } + } void hydraBase::setConstraintJacobians( ){ /*! * Set the constraint Jacobians values */ + + const unsigned int numUnknowns = getNumUnknowns( ); + + const unsigned int numConstraints = getNumConstraints( ); + + auto constraintJacobians = get_setDataStorage_constraintJacobians( ); + + constraintJacobians.zero( numConstraints * numUnknowns ); + + unsigned int offset = 0; + + for ( auto v = getResidualClasses( )->begin( ); v != getResidualClasses( )->end( ); v++ ){ + + if ( ( *( *v )->getNumConstraints( ) ) > 0 ){ + + std::copy( ( *v )->getConstraintJacobians( )->begin( ), ( *v )->getConstraintJacobians( )->end( ), constraintJacobians.value->begin( ) + offset ); + + offset += ( *v )->getConstraintJacobians( )->size( ); + + } + + } + } errorOut dummyMaterialModel( floatVector &stress, floatVector &statev, floatMatrix &ddsdde, floatType &SSE, floatType &SPD, diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index 4bd9b76..d21ddcc 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -6191,3 +6191,156 @@ BOOST_AUTO_TEST_CASE( test_setDataStorageIteration2, * boost::unit_test::toleran BOOST_TEST( !residual._myVectorData.first ); } + +BOOST_AUTO_TEST_CASE( test_hydraBase_setConstraints, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ + + class residualMock : public tardigradeHydra::residualBase{ + + public: + + using tardigradeHydra::residualBase::residualBase; + + residualMock( tardigradeHydra::hydraBase *h, const unsigned int neq, const unsigned int ncon ) : tardigradeHydra::residualBase( h, neq ){ + + setNumConstraints( ncon ); + + } + + protected: + + virtual void setConstraints( ){ + + auto constraints = get_setDataStorage_constraints( ); + + constraints.zero( *getNumConstraints( ) ); + + for ( unsigned int i = 0; i < *getNumConstraints( ); i++ ){ + + ( *constraints.value )[ i ] = ( *getNumConstraints( ) ) + 0.1 * i; + + } + + } + + virtual void setConstraintJacobians( ){ + + const unsigned int numUnknowns = hydra->getNumUnknowns( ); + + auto constraintJacobians = get_setDataStorage_constraintJacobians( ); + + constraintJacobians.zero( ( *getNumConstraints( ) ) * numUnknowns ); + + for ( unsigned int i = 0; i < *getNumConstraints( ); i++ ){ + + for ( unsigned int j = 0; j < numUnknowns; j++ ){ + + ( *constraintJacobians.value )[ numUnknowns * i + j ] = ( *getNumConstraints( ) ) + 0.1 * ( i + j ); + + } + + } + + } + + }; + + class hydraBaseMock : public tardigradeHydra::hydraBase{ + + public: + + residualMock r1; + + residualMock r2; + + residualMock r3; + + unsigned int s1 = 36; + + unsigned int s2 = 2; + + unsigned int s3 = 3; + + unsigned int n1 = 2; + + unsigned int n2 = 0; + + unsigned int n3 = 5; + + using tardigradeHydra::hydraBase::hydraBase; + + using tardigradeHydra::hydraBase::setResidualClasses; + + virtual void setResidualClasses( ){ + + r1 = residualMock( this, s1, n1 ); + + r2 = residualMock( this, s2, n2 ); + + r3 = residualMock( this, s3, n3 ); + + std::vector< tardigradeHydra::residualBase* > residuals( 3 ); + + residuals[ 0 ] = &r1; + + residuals[ 1 ] = &r2; + + residuals[ 2 ] = &r3; + + setResidualClasses( residuals ); + + } + + }; + + 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 ); + + const floatVector constraint_answer = { 2.0, 2.1, 5.0, 5.1, 5.2, 5.3, 5.4 }; + + const floatVector constraintJacobian_answer = { 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, + 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, + 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, + 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, + 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, + 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, 9.3, + 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, 9.3, 9.4 }; + + BOOST_TEST( constraint_answer == *hydra.getConstraints( ), CHECK_PER_ELEMENT ); + + BOOST_TEST( constraintJacobian_answer == *hydra.getConstraintJacobians( ), CHECK_PER_ELEMENT ); + +} From f35a2cb3e53ad30604370aa59c25646e1eef6342 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 08:04:56 -0600 Subject: [PATCH 04/10] TST: Completed tests for the assembly of the KKT matrix --- src/cpp/tardigrade_hydra.cpp | 12 +- src/cpp/tardigrade_hydra.h | 6 +- src/cpp/tests/test_tardigrade_hydra.cpp | 151 ++++++++++++++++++++++++ 3 files changed, 160 insertions(+), 9 deletions(-) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index f6fbaa7..1c81d87 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2203,7 +2203,7 @@ namespace tardigradeHydra{ } - void hydraBase::assembleKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ){ + void hydraBase::assembleKKTMatrix( floatVector &KKTMatrix, const std::vector< bool > &active_constraints ){ /*! * Assemble the Karush-Kuhn-Tucker matrix for an inequality constrained Newton-Raphson solve * @@ -2226,7 +2226,7 @@ namespace tardigradeHydra{ for ( unsigned int I = 0; I < numUnknowns; I++ ){ - KKTMatrix[ ( numUnknowns + numConstraints ) * I + I ] = ( *getMuk( ) ); + KKTMatrix[ ( numUnknowns + numConstraints ) * I + I ] += ( *getMuk( ) ); } @@ -2252,7 +2252,7 @@ namespace tardigradeHydra{ } - void hydraBase::updateKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ){ + void hydraBase::updateKKTMatrix( floatVector &KKTMatrix, const std::vector< bool > &active_constraints ){ /*! * Update the KKTMatrix if the active constraints have changed * @@ -2294,7 +2294,7 @@ namespace tardigradeHydra{ } - void hydraBase::assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, std::vector< bool > &active_constraints ){ + void hydraBase::assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, const std::vector< bool > &active_constraints ){ /*! * Assemble the right hand side vector for the KKT matrix * @@ -2314,10 +2314,10 @@ namespace tardigradeHydra{ Eigen::Map< Eigen::Vector< floatType, -1 > > RHS( KKTRHSVector.data( ), ( numUnknowns + numConstraints ), ( numUnknowns + numConstraints ) ); Eigen::Map< const Eigen::Vector< floatType, -1 > > R( getResidual( )->data( ), numUnknowns ); - + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > J( getFlatJacobian( )->data( ), numUnknowns, numUnknowns ); - RHS.head( numUnknowns ) += ( J.transpose( ) * ( R + J * _dx ) + ( *getMuk( ) ) * _dx ).eval( ); + RHS.head( numUnknowns ) = ( J.transpose( ) * ( R + J * _dx ) + ( *getMuk( ) ) * _dx ).eval( ); for ( unsigned int i = 0; i < numConstraints; i++ ){ diff --git a/src/cpp/tardigrade_hydra.h b/src/cpp/tardigrade_hydra.h index a038279..812ae23 100644 --- a/src/cpp/tardigrade_hydra.h +++ b/src/cpp/tardigrade_hydra.h @@ -1619,11 +1619,11 @@ namespace tardigradeHydra{ virtual tardigradeHydra::hydraBase::setDataStorageIteration get_setDataStorage_stress( ); - virtual void assembleKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ); + virtual void assembleKKTMatrix( floatVector &KKTMatrix, const std::vector< bool > &active_constraints ); - virtual void updateKKTMatrix( floatVector &KKTMatrix, std::vector< bool > &active_constraints ); + virtual void updateKKTMatrix( floatVector &KKTMatrix, const std::vector< bool > &active_constraints ); - virtual void assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, std::vector< bool > &active_constraints ); + virtual void assembleKKTRHSVector( const floatVector &dx, floatVector &KKTRHSVector, const std::vector< bool > &active_constraints ); virtual void solveConstrainedQP( floatVector &x, const unsigned int kmax=100 ); diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index d21ddcc..23bd823 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -6344,3 +6344,154 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_setConstraints, * boost::unit_test::toleran BOOST_TEST( constraintJacobian_answer == *hydra.getConstraintJacobians( ), CHECK_PER_ELEMENT ); } + +BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTMatrix, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ + + class hydraBaseMock : public tardigradeHydra::hydraBase{ + + public: + + using tardigradeHydra::hydraBase::hydraBase; + + using tardigradeHydra::hydraBase::setResidualClasses; + + floatVector initialUnknownVector = { 2, 1 }; + + virtual void public_assembleKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ){ + + assembleKKTMatrix( K, active_constraints ); + + } + + virtual void public_updateKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ){ + + updateKKTMatrix( K, active_constraints ); + + } + + virtual void public_assembleKKTRHSVector( const floatVector &dx, floatVector &RHS, const std::vector< bool > &active_constraints ){ + + assembleKKTRHSVector( dx, RHS, active_constraints ); + + } + + protected: + + virtual void setConstraints( ) override{ + + auto constraints = get_setDataStorage_constraints( ); + + *constraints.value = { 2, 6, 2, 0, 0 }; + + for ( unsigned int i = 0; i < 5; i++ ){ + + for ( unsigned int j = 0; j < 2; j++ ){ + + ( *constraints.value )[ i ] += ( *getConstraintJacobians( ) )[ 2 * i + j ] * initialUnknownVector[ j ]; + + } + + } + + } + + virtual void setConstraintJacobians( ) override{ + + auto constraintJacobians = get_setDataStorage_constraintJacobians( ); + + *constraintJacobians.value = { 1, -2, + -1, -2, + -1, 2, + 1, 0, + 0, 1 }; + + } + + virtual void formNonLinearDerivatives( ) override{ + + floatVector jacobian = { std::pow( 2, 0.5 ), 0, 0, std::pow( 2, 0.5 ) }; + + tardigradeHydra::unit_test::hydraBaseTester::set_flatJacobian( *this, jacobian ); + + } + + virtual const unsigned int getNumUnknowns( ) override{ return initialUnknownVector.size( ); } + + }; + + 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 ); + + hydra.setMuk( 0.1 ); + + floatVector result_KKT; + std::vector< bool > active_constraints( 5, false ); + + floatVector answer1_KKTMatrix = { 2.1, 0., 0., 0., 0., 0., 0., + 0., 2.1, 0., 0., 0., 0., 0., + 0., 0., 1., 0., 0., 0., 0., + 0., 0., 0., 1., 0., 0., 0., + 0., 0., 0., 0., 1., 0., 0., + 0., 0., 0., 0., 0., 1., 0., + 0., 0., 0., 0., 0., 0., 1. }; + + floatVector answer2_KKTMatrix = { 2.1, 0., 0., 0., -1., 0., 0., + 0., 2.1, 0., 0., 2., 0., 1., + 0., 0., 1., 0., 0., 0., 0., + 0., 0., 0., 1., 0., 0., 0., + -1., 2., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 1., 0., + 0., 1., 0., 0., 0., 0., 0. }; + + hydra.public_assembleKKTMatrix( result_KKT, active_constraints ); + + BOOST_TEST( answer1_KKTMatrix == result_KKT, CHECK_PER_ELEMENT ); + + active_constraints[ 2 ] = true; + active_constraints[ 4 ] = true; + + hydra.public_updateKKTMatrix( result_KKT, active_constraints ); + + BOOST_TEST( answer2_KKTMatrix == result_KKT, CHECK_PER_ELEMENT ); + + result_KKT.clear( ); + + hydra.public_assembleKKTMatrix( result_KKT, active_constraints ); + + BOOST_TEST( answer2_KKTMatrix == result_KKT, CHECK_PER_ELEMENT ); + +} From 120f10d926587b96f4892a72648ca586f49805c2 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 10:10:46 -0600 Subject: [PATCH 05/10] FEAT: Added the ability and tests for the construction of the KKT problem RHS vector --- src/cpp/tardigrade_hydra.cpp | 4 +- src/cpp/tests/test_tardigrade_hydra.cpp | 181 +++++++++++++++++++++--- 2 files changed, 167 insertions(+), 18 deletions(-) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index 1c81d87..84865cc 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2213,7 +2213,7 @@ namespace tardigradeHydra{ const unsigned int numUnknowns = getNumUnknowns( ); - const unsigned int numConstraints = active_constraints.size( ); + const unsigned int numConstraints = getNumConstraints( ); KKTMatrix = floatVector( ( numUnknowns + numConstraints ) * ( numUnknowns + numConstraints ), 0 ); @@ -2262,7 +2262,7 @@ namespace tardigradeHydra{ const unsigned int numUnknowns = getNumUnknowns( ); - const unsigned int numConstraints = active_constraints.size( ); + const unsigned int numConstraints = getNumConstraints( ); Eigen::Map< Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > K( KKTMatrix.data( ), ( numUnknowns + numConstraints ), ( numUnknowns + numConstraints ) ); diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index 23bd823..b0d417f 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -6409,7 +6409,7 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTMatrix, * boost::unit_test::tole virtual void formNonLinearDerivatives( ) override{ - floatVector jacobian = { std::pow( 2, 0.5 ), 0, 0, std::pow( 2, 0.5 ) }; + floatVector jacobian = { std::pow( 2, 0.5 ), 0.4, -0.1, std::pow( 2, 0.5 ) }; tardigradeHydra::unit_test::hydraBaseTester::set_flatJacobian( *this, jacobian ); @@ -6417,6 +6417,8 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTMatrix, * boost::unit_test::tole virtual const unsigned int getNumUnknowns( ) override{ return initialUnknownVector.size( ); } + virtual const unsigned int getNumConstraints( ) override{ return 5; } + }; floatType time = 1.1; @@ -6461,21 +6463,21 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTMatrix, * boost::unit_test::tole floatVector result_KKT; std::vector< bool > active_constraints( 5, false ); - floatVector answer1_KKTMatrix = { 2.1, 0., 0., 0., 0., 0., 0., - 0., 2.1, 0., 0., 0., 0., 0., - 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 0., 0., 1. }; - - floatVector answer2_KKTMatrix = { 2.1, 0., 0., 0., -1., 0., 0., - 0., 2.1, 0., 0., 2., 0., 1., - 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 1., 0., 0., 0., - -1., 2., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 1., 0., - 0., 1., 0., 0., 0., 0., 0. }; + floatVector answer1_KKTMatrix = { 2.11 , 0.42426407, 0., 0., 0., 0., 0., + 0.42426407, 2.26 , 0., 0., 0., 0., 0., + 0. , 0. , 1., 0., 0., 0., 0., + 0. , 0. , 0., 1., 0., 0., 0., + 0. , 0. , 0., 0., 1., 0., 0., + 0. , 0. , 0., 0., 0., 1., 0., + 0. , 0. , 0., 0., 0., 0., 1. }; + + floatVector answer2_KKTMatrix = { 2.11 , 0.42426407, 0., 0., -1., 0., 0., + 0.42426407, 2.26 , 0., 0., 2., 0., 1., + 0. , 0. , 1., 0., 0., 0., 0., + 0. , 0. , 0., 1., 0., 0., 0., + -1. , 2. , 0., 0., 0., 0., 0., + 0. , 0. , 0., 0., 0., 1., 0., + 0. , 1. , 0., 0., 0., 0., 0. }; hydra.public_assembleKKTMatrix( result_KKT, active_constraints ); @@ -6495,3 +6497,150 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTMatrix, * boost::unit_test::tole BOOST_TEST( answer2_KKTMatrix == result_KKT, CHECK_PER_ELEMENT ); } + +BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTRHSVector, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ + + class hydraBaseMock : public tardigradeHydra::hydraBase{ + + public: + + using tardigradeHydra::hydraBase::hydraBase; + + using tardigradeHydra::hydraBase::setResidualClasses; + + floatVector initialUnknownVector = { 2, 1 }; + + virtual void public_assembleKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ){ + + assembleKKTMatrix( K, active_constraints ); + + } + + virtual void public_updateKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ){ + + updateKKTMatrix( K, active_constraints ); + + } + + virtual void public_assembleKKTRHSVector( const floatVector &dx, floatVector &RHS, const std::vector< bool > &active_constraints ){ + + assembleKKTRHSVector( dx, RHS, active_constraints ); + + } + + protected: + + virtual void setConstraints( ) override{ + + auto constraints = get_setDataStorage_constraints( ); + + *constraints.value = { 2, 6, 2, 0, 0 }; + + for ( unsigned int i = 0; i < 5; i++ ){ + + for ( unsigned int j = 0; j < 2; j++ ){ + + ( *constraints.value )[ i ] += ( *getConstraintJacobians( ) )[ 2 * i + j ] * initialUnknownVector[ j ]; + + } + + } + + } + + virtual void setConstraintJacobians( ) override{ + + auto constraintJacobians = get_setDataStorage_constraintJacobians( ); + + *constraintJacobians.value = { 1, -2, + -1, -2, + -1, 2, + 1, 0, + 0, 1 }; + + } + + virtual void formNonLinearResidual( ) override{ + + floatVector residual = { 1., 2. }; + + tardigradeHydra::unit_test::hydraBaseTester::set_residual( *this, residual ); + + } + + virtual void formNonLinearDerivatives( ) override{ + + floatVector jacobian = { std::pow( 2, 0.5 ), 0.4, -0.1, std::pow( 2, 0.5 ) }; + + tardigradeHydra::unit_test::hydraBaseTester::set_flatJacobian( *this, jacobian ); + + } + + virtual const unsigned int getNumUnknowns( ) override{ return initialUnknownVector.size( ); } + + 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 ); + + floatVector dx = { -0.2, 1.4 }; + + hydra.setMuk( 0.1 ); + + floatVector result_KKTRHSVector; + std::vector< bool > active_constraints( 5, false ); + + floatVector answer1_KKTRHSVector = { 1.38618326, 6.30757431, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 }; + + floatVector answer2_KKTRHSVector = { 1.38618326, 6.30757431, 0.0 , 0.0 , 5. , 0.0 , 2.4 }; + + hydra.public_assembleKKTRHSVector( dx, result_KKTRHSVector, active_constraints ); + + BOOST_TEST( answer1_KKTRHSVector == result_KKTRHSVector, CHECK_PER_ELEMENT ); + + active_constraints[ 2 ] = true; + active_constraints[ 4 ] = true; + + result_KKTRHSVector.clear( ); + + hydra.public_assembleKKTRHSVector( dx, result_KKTRHSVector, active_constraints ); + + BOOST_TEST( answer2_KKTRHSVector == result_KKTRHSVector, CHECK_PER_ELEMENT ); + +} From c5656f0dc965c339e205bea6f45662786c5062d3 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 13:38:18 -0600 Subject: [PATCH 06/10] TST: Added tests for solving a quadratic problem using the active set method --- src/cpp/tardigrade_hydra.cpp | 71 +++++---- src/cpp/tardigrade_hydra.h | 2 + src/cpp/tests/test_tardigrade_hydra.cpp | 197 ++++++++++++++++++++++++ 3 files changed, 239 insertions(+), 31 deletions(-) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index 84865cc..6c4b9ab 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2337,6 +2337,17 @@ namespace tardigradeHydra{ } + void hydraBase::initializeActiveConstraints( std::vector< bool > &active_constraints ){ + /*! + * Initialize the active constraint vector + * + * \param &active_constraints: The current constraints that are active + */ + + active_constraints = std::vector< bool >( getNumConstraints( ), false ); + + } + void hydraBase::solveConstrainedQP( floatVector &dx, const unsigned int kmax ){ /*! * Solve the constrained QP problem to estimate the desired step size @@ -2353,7 +2364,8 @@ namespace tardigradeHydra{ floatVector RHS; - std::vector< bool > active_constraints( numConstraints, false ); + std::vector< bool > active_constraints; + initializeActiveConstraints( active_constraints ); assembleKKTRHSVector( dx, RHS, active_constraints ); @@ -2365,7 +2377,7 @@ namespace tardigradeHydra{ floatVector y( numUnknowns + numConstraints, 0 ); - floatVector ck( numConstraints, 0 ); + floatVector ck = *getConstraints( ); floatVector ctilde( numConstraints, 0 ); @@ -2383,41 +2395,23 @@ namespace tardigradeHydra{ } - const floatVector *c0 = getConstraints( ); - - const floatVector *A = getConstraintJacobians( ); - - Eigen::Map< Eigen::Vector< floatType, -1 > > _dx( dx.data( ), numUnknowns ); - - Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > _K( K.data( ), numConstraints + numUnknowns, numConstraints + numUnknowns ); - - Eigen::Map< const Eigen::Vector< floatType, -1 > > _RHS( RHS.data( ), numConstraints + numUnknowns ); - - Eigen::Map< Eigen::Vector< floatType, -1 > > _y( y.data( ), numConstraints + numUnknowns ); - - Eigen::Map< Eigen::Vector< floatType, -1 > > _ck( ck.data( ), numConstraints ); - - Eigen::Map< Eigen::Vector< floatType, -1 > > _ctilde( ctilde.data( ), numConstraints ); - - Eigen::Map< Eigen::Vector< floatType, -1 > > _negp( negp.data( ), numUnknowns ); - - Eigen::Map< const Eigen::Vector< floatType, -1 > > _c0( c0->data( ), numConstraints ); - - Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > _A( A->data( ), numConstraints, numUnknowns ); - Eigen::Map< const Eigen::Vector< floatType, -1 > > _P( P.data( ), numUnknowns + numConstraints ); tardigradeVectorTools::solverType< floatType > linearSolver; while ( k < kmax ){ + Eigen::Map< const Eigen::Matrix< floatType, -1, -1, Eigen::RowMajor > > _K( K.data( ), numConstraints + numUnknowns, numConstraints + numUnknowns ); + Eigen::Map< const Eigen::Vector< floatType, -1 > > _RHS( RHS.data( ), numConstraints + numUnknowns ); + Eigen::Map< Eigen::Vector< floatType, -1 > > _y( y.data( ), numConstraints + numUnknowns ); + linearSolver = tardigradeVectorTools::solverType< floatType >( _P.asDiagonal( ) * _K ); _y = linearSolver.solve( _P.asDiagonal( ) * _RHS ); - negp = floatVector( y.begin( ), y.begin( ) + numUnknowns ); + std::copy( y.begin( ), y.begin( ) + numUnknowns, negp.begin( ) ); - lambda = floatVector( y.begin( ) + numUnknowns, y.end( ) ); + std::copy( y.begin( ) + numUnknowns, y.end( ), lambda.begin( ) ); if ( tardigradeVectorTools::l2norm( negp ) <= tol ){ @@ -2459,9 +2453,14 @@ namespace tardigradeHydra{ } else{ - _ck = ( _c0 + _A * _dx ).eval( ); - - _ctilde = ( _c0 + _A * ( _dx - _negp ) ).eval( ); + ck = *getConstraints( ); + ctilde = *getConstraints( ); + for ( unsigned int i = 0; i < numConstraints; i++ ){ + for ( unsigned int j = 0; j < numUnknowns; j++ ){ + ck[ i ] += ( *getConstraintJacobians( ) )[ numUnknowns * i + j ] * dx[ j ]; + ctilde[ i ] += ( *getConstraintJacobians( ) )[ numUnknowns * i + j ] * ( dx[ j ] - negp[ j ] ); + } + } floatType alpha = 1.0; @@ -2501,12 +2500,22 @@ namespace tardigradeHydra{ dx -= alpha * negp; - updateKKTMatrix( K, active_constraints ); + } + + updateKKTMatrix( K, active_constraints ); - assembleKKTRHSVector( dx, RHS, active_constraints ); + assembleKKTRHSVector( dx, RHS, active_constraints ); + + for ( unsigned int i = 0; i < ( numUnknowns + numConstraints ); i++ ){ + + P[ i ] = 1 / std::max( std::fabs( *std::max_element( K.begin( ) + ( numUnknowns + numConstraints ) * i, + K.begin( ) + ( numUnknowns + numConstraints ) * ( i + 1 ), + [ ]( const floatType &a, const floatType &b ){ return std::fabs( a ) < std::fabs( b ); } ) ), 1e-15 ); } + if ( k > 4 ){ throw std::runtime_error("derp"); } + k++; } diff --git a/src/cpp/tardigrade_hydra.h b/src/cpp/tardigrade_hydra.h index 812ae23..297881e 100644 --- a/src/cpp/tardigrade_hydra.h +++ b/src/cpp/tardigrade_hydra.h @@ -1631,6 +1631,8 @@ namespace tardigradeHydra{ virtual void setConstraintJacobians( ); + virtual void initializeActiveConstraints( std::vector< bool > &active_constraints ); + private: // Friend classes diff --git a/src/cpp/tests/test_tardigrade_hydra.cpp b/src/cpp/tests/test_tardigrade_hydra.cpp index b0d417f..9df2c37 100644 --- a/src/cpp/tests/test_tardigrade_hydra.cpp +++ b/src/cpp/tests/test_tardigrade_hydra.cpp @@ -6644,3 +6644,200 @@ BOOST_AUTO_TEST_CASE( test_hydraBase_assembleKKTRHSVector, * boost::unit_test::t BOOST_TEST( answer2_KKTRHSVector == result_KKTRHSVector, CHECK_PER_ELEMENT ); } + +BOOST_AUTO_TEST_CASE( test_hydraBase_solveConstrainedQP, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){ + + class hydraBaseMock : public tardigradeHydra::hydraBase{ + + public: + + using tardigradeHydra::hydraBase::hydraBase; + + using tardigradeHydra::hydraBase::setResidualClasses; + + floatVector initialUnknownVector = { 2, 0 }; + + void public_solveConstrainedQP( floatVector &dx ){ + + solveConstrainedQP( dx ); + + } + + protected: + + virtual void assembleKKTRHSVector( const floatVector &dx, floatVector &RHS, const std::vector< bool > &active_constraints ) override{ + + RHS = floatVector( 7, 0 ); + + RHS[ 0 ] = 2 * ( initialUnknownVector[ 0 ] + dx[ 0 ] - 1.0 ); + RHS[ 1 ] = 2 * ( initialUnknownVector[ 1 ] + dx[ 1 ] - 2.5 ); + + RHS[ 2 + 0 ] = ( initialUnknownVector[ 0 ] + dx[ 0 ] ) - 2 * ( initialUnknownVector[ 1 ] + dx[ 1 ] ) + 2; + RHS[ 2 + 1 ] = -( initialUnknownVector[ 0 ] + dx[ 0 ] ) - 2 * ( initialUnknownVector[ 1 ] + dx[ 1 ] ) + 6; + RHS[ 2 + 2 ] = -( initialUnknownVector[ 0 ] + dx[ 0 ] ) + 2 * ( initialUnknownVector[ 1 ] + dx[ 1 ] ) + 2; + RHS[ 2 + 3 ] = initialUnknownVector[ 0 ] + dx[ 0 ]; + RHS[ 2 + 4 ] = initialUnknownVector[ 1 ] + dx[ 1 ]; + + for ( unsigned int i = 0; i < active_constraints.size( ); i++ ){ + + if ( !active_constraints[ i ] ){ + + RHS[ 2 + i ] = 0; + + } + + } + + } + + virtual void assembleKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ) override{ + + const unsigned int numConstraints = getNumConstraints( ); + + K = floatVector( ( 2 + 5 ) * ( 2 + 5 ), 0 ); + + K[ 7 * 0 + 0 ] = 2; + K[ 7 * 1 + 1 ] = 2; + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( active_constraints[ i ] ){ + + K[ 7 * 0 + i + 2 ] = ( *getConstraintJacobians( ) )[ 2 * i + 0 ]; + K[ 7 * 1 + i + 2 ] = ( *getConstraintJacobians( ) )[ 2 * i + 1 ]; + + K[ 7 * ( i + 2 ) + 0 ] = ( *getConstraintJacobians( ) )[ 2 * i + 0 ]; + K[ 7 * ( i + 2 ) + 1 ] = ( *getConstraintJacobians( ) )[ 2 * i + 1 ]; + + } + else{ + + K[ 7 * ( i + 2 ) + i + 2 ] = 1; + + } + + } + + } + + virtual void updateKKTMatrix( floatVector &K, const std::vector< bool > &active_constraints ) override{ + + const unsigned int numConstraints = getNumConstraints( ); + + for ( unsigned int i = 0; i < numConstraints; i++ ){ + + if ( active_constraints[ i ] ){ + + K[ 7 * 0 + i + 2 ] = ( *getConstraintJacobians( ) )[ 2 * i + 0 ]; + K[ 7 * 1 + i + 2 ] = ( *getConstraintJacobians( ) )[ 2 * i + 1 ]; + + K[ 7 * ( i + 2 ) + 0 ] = ( *getConstraintJacobians( ) )[ 2 * i + 0 ]; + K[ 7 * ( i + 2 ) + 1 ] = ( *getConstraintJacobians( ) )[ 2 * i + 1 ]; + + K[ 7 * ( i + 2 ) + i + 2 ] = 0; + + } + else{ + + K[ 7 * 0 + i + 2 ] = 0; + K[ 7 * 1 + i + 2 ] = 0; + + K[ 7 * ( i + 2 ) + 0 ] = 0; + K[ 7 * ( i + 2 ) + 1 ] = 0; + + K[ 7 * ( i + 2 ) + i + 2 ] = 1; + + } + + } + + } + + virtual void initializeActiveConstraints( std::vector< bool > &active_constraints ) override{ + + active_constraints = { false, false, true, false, true }; + + } + + virtual void setConstraints( ) override{ + + auto constraints = get_setDataStorage_constraints( ); + + *constraints.value = { 2, 6, 2, 0, 0 }; + + for ( unsigned int i = 0; i < 5; i++ ){ + + for ( unsigned int j = 0; j < 2; j++ ){ + + ( *constraints.value )[ i ] += ( *getConstraintJacobians( ) )[ 2 * i + j ] * initialUnknownVector[ j ]; + + } + + } + + } + + virtual void setConstraintJacobians( ) override{ + + auto constraintJacobians = get_setDataStorage_constraintJacobians( ); + + *constraintJacobians.value = { 1, -2, + -1, -2, + -1, 2, + 1, 0, + 0, 1 }; + + } + + virtual const unsigned int getNumUnknowns( ) override{ return initialUnknownVector.size( ); } + + 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 ); + + floatVector result = { 0, 0 }; + + floatVector answer = { 1.4, 1.7 }; + + hydra.public_solveConstrainedQP( result ); + + BOOST_TEST( ( result + hydra.initialUnknownVector ) == answer, CHECK_PER_ELEMENT ); + +} From 9ad9bdafb4179df18202e9ca6f4ddc2f9e8e8ca4 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 13:51:38 -0600 Subject: [PATCH 07/10] MAINT: Updated changelog --- docs/sphinx/changelog.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/sphinx/changelog.rst b/docs/sphinx/changelog.rst index 62bbbb8..aa330b6 100644 --- a/docs/sphinx/changelog.rst +++ b/docs/sphinx/changelog.rst @@ -34,6 +34,7 @@ Internal Changes - Removed set_varname and replaced with setDataStorage for hydra Peryzna isotropic-kenmatic J2 visoplasticity (:pull:`161`). By `Nathan Miller`_. - Removed set_varname and replaced with setDataStorage for hydra mass-change rate (:pull:`162`). By `Nathan Miller`_. - Enabled passing the tests when higher levels of optimization are used (:pull:`163`). By `Nathan Miller`_. +- Added an active set solver for quadratic problems (:pull:`167`). By `Nathan Miller`_. Bug Fixes ========= From d48de5f717ba79b7b574deba6cbcb31b94b59b40 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 14:48:06 -0600 Subject: [PATCH 08/10] MAINT: Removed debugging error --- src/cpp/tardigrade_hydra.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/cpp/tardigrade_hydra.cpp b/src/cpp/tardigrade_hydra.cpp index 6c4b9ab..f022687 100644 --- a/src/cpp/tardigrade_hydra.cpp +++ b/src/cpp/tardigrade_hydra.cpp @@ -2514,8 +2514,6 @@ namespace tardigradeHydra{ } - if ( k > 4 ){ throw std::runtime_error("derp"); } - k++; } From b06b9129556b0f699a38ea08a75344ddbf61ee6c Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Tue, 3 Sep 2024 15:12:04 -0600 Subject: [PATCH 09/10] BUG: Corrected setting the number of constraints --- ...grade_hydraMicromorphicDruckerPragerPlasticityOptimization.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cpp/tardigrade_hydraMicromorphicDruckerPragerPlasticityOptimization.h b/src/cpp/tardigrade_hydraMicromorphicDruckerPragerPlasticityOptimization.h index b1a41f4..f8cadc4 100644 --- a/src/cpp/tardigrade_hydraMicromorphicDruckerPragerPlasticityOptimization.h +++ b/src/cpp/tardigrade_hydraMicromorphicDruckerPragerPlasticityOptimization.h @@ -118,6 +118,8 @@ namespace tardigradeHydra{ setPenaltyIndices( positive_indices ); //!< The indices of values in the unknown vector which should be penalized for being negative + setNumConstraints( 10 ); + } virtual void suggestInitialIterateValues( std::vector< unsigned int > &indices, From 0134b5226f63e3b10a8ce73125c1863790c8faf2 Mon Sep 17 00:00:00 2001 From: Nathan Miller Date: Wed, 4 Sep 2024 14:15:25 -0600 Subject: [PATCH 10/10] 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 ); + +}