FEAT: Added option for the constrained incrementation solver

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
 ================

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 ){

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( );

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 );
+
+}