feat(materials): add substepper for finite strain materials#32
feat(materials): add substepper for finite strain materials#32alexdummer wants to merge 13 commits intoMAteRialMOdelingToolbox:next_v26.05from
Conversation
There was a problem hiding this comment.
Pull request overview
This WIP pull request introduces a substepping capability for finite strain material models through a decorator pattern, enabling improved numerical stability and accuracy for large deformation increments through analytical tangent accumulation.
Changes:
- Adds a new template decorator class
MarmotMaterialFiniteStrainSubsteppedthat wraps existing material models with time substepping functionality - Implements
computeStressWithSensitivitiesmethod with default numerical differentiation for sensitivity matrices required by substepping - Fixes a bug in
centralDifferencewhere the Jacobian matrix was allocated before knowing the function output size
Reviewed changes
Copilot reviewed 8 out of 8 changed files in this pull request and generated 19 comments.
Show a summary per file
| File | Description |
|---|---|
| modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h | New template class implementing substepping decorator with analytical tangent accumulation |
| modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h | Adds StateSensitivities struct and computeStressWithSensitivities method with numerical differentiation |
| modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp | Fixes bug in centralDifference by deferring matrix allocation until output size is known |
| modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2PlasticityRegistration.cpp | Registers substepped variant of J2 plasticity model |
| modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2Plasticity.cpp | Changes initialization method from .eye() to memcpy for consistency |
| modules/core/MarmotMathCore/module.cmake | Adds public headers for Fastor tensor basics and numerical differentiation |
| examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/example.cpp | Demonstrates usage of substepped material with mixed loading controls |
| examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/Makefile | Build configuration for the example |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| * @brief Extended stress update computing sensitivities for substepping. | ||
| * * The default implementation uses Central Finite Differences. |
There was a problem hiding this comment.
...s/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h
Outdated
Show resolved
Hide resolved
| /* std::cout << "deformation.F: \n" << deformation.F << std::endl; */ | ||
| /* std::cout << "stateOld: \n" << Map<const VectorXd>(stateOld.data(), nState) << std::endl; */ | ||
| // Run the actual update once to get Stress and New State (at t_n+1) | ||
| computeStress( response, tangents, deformation, timeIncrement ); | ||
|
|
||
| /* std::cout << "Base Stress after update: \n" << response.tau << std::endl; */ | ||
|
|
||
| // If there are no state variables (elastic), sensitivities are zero/empty | ||
| if ( nState == 0 ) | ||
| return; | ||
|
|
||
| // --- Numerical Differentiation Setup --- | ||
|
|
||
| // Prepare sensitivity matrices | ||
| sensitivities.dState_dF.resize( nState, 9 ); | ||
| sensitivities.dState_dStateOld.resize( nState, nState ); | ||
| sensitivities.dStress_dStateOld.resize( 9, nState ); | ||
|
|
||
| // We need temporary storage for the FD routines to avoid corrupting the actual results | ||
| // or the input pointers. | ||
|
|
||
| // 2. Compute d(State_new) / d(F_new) | ||
| // We perturb F, run computeStress (starting from fixed StateOld), and measure State_new. | ||
| auto func_dState_dF = [&]( const VectorXd& F_vec ) -> VectorXd { | ||
| Deformation< 3 > defPerturbed; | ||
| std::memcpy( defPerturbed.F.data(), F_vec.data(), 9 * sizeof( double ) ); | ||
|
|
||
| ConstitutiveResponse< 3 > respPert; | ||
| std::vector< double > stateTemp = stateOld; // Copy from OLD | ||
| respPert.stateVars = stateTemp.data(); | ||
|
|
||
| AlgorithmicModuli< 3 > tanPert; // dummy | ||
|
|
||
| // Run update | ||
| computeStress( respPert, tanPert, defPerturbed, timeIncrement ); | ||
|
|
||
| // Return NEW state | ||
| Map< VectorXd > res( respPert.stateVars, nState ); | ||
| return res; | ||
| }; | ||
|
|
||
| Map< const VectorXd > F_map( deformation.F.data(), 9 ); | ||
| sensitivities.dState_dF = forwardDifference( func_dState_dF, F_map ); | ||
|
|
||
| // print dState_dF if det is nan | ||
| if ( sensitivities.dState_dF.hasNaN() ) { | ||
| std::cout << "Computed dState/dF has NaN!" << std::endl; | ||
| std::cout << "F: \n" << deformation.F << std::endl; | ||
| std::cout << "StateOld: \n" << Map< const VectorXd >( stateOld.data(), nState ) << std::endl; | ||
| } | ||
| /* std::cout << "Computed dState/dF via FD: \n" << sensitivities.dState_dF << std::endl; */ | ||
|
|
||
| // 3. Compute d(State_new) / d(State_old) AND d(Stress) / d(State_old) | ||
| // We perturb State_old, run computeStress (fixed F), and measure State_new and Stress. | ||
|
|
||
| // Since centralDifference returns one matrix, we do this in two passes or combine outputs. | ||
| // Let's do two passes for clarity, though it doubles the cost of this specific block. | ||
| // Given the prompt constraints, we use the provided interface. | ||
|
|
||
| // A. d(State_new) / d(State_old) | ||
| auto func_dState_dStateOld = [&]( const VectorXd& StateOld_vec ) -> VectorXd { | ||
| ConstitutiveResponse< 3 > respPert; | ||
| std::vector< double > stateTemp( nState ); | ||
| std::memcpy( stateTemp.data(), StateOld_vec.data(), nState * sizeof( double ) ); | ||
| respPert.stateVars = stateTemp.data(); | ||
|
|
||
| AlgorithmicModuli< 3 > tanPert; // dummy | ||
|
|
||
| // Run update with perturbed old state | ||
| computeStress( respPert, tanPert, deformation, timeIncrement ); | ||
|
|
||
| Map< VectorXd > res( respPert.stateVars, nState ); | ||
| return res; | ||
| }; | ||
|
|
||
| Map< const VectorXd > StateOld_map( stateOld.data(), nState ); | ||
| sensitivities.dState_dStateOld = forwardDifference( func_dState_dStateOld, StateOld_map ); | ||
|
|
||
| // print dState_dF if det is nan | ||
| if ( sensitivities.dState_dStateOld.hasNaN() ) { | ||
| std::cout << "Computed dState/dStateOld has NaN!" << std::endl; | ||
| std::cout << "F: \n" << deformation.F << std::endl; | ||
| std::cout << "StateOld: \n" << Map< const VectorXd >( stateOld.data(), nState ) << std::endl; | ||
| } | ||
| // B. d(Stress) / d(State_old) | ||
| auto func_dStress_dStateOld = [&]( const VectorXd& StateOld_vec ) -> VectorXd { | ||
| ConstitutiveResponse< 3 > respPert; | ||
| std::vector< double > stateTemp( nState ); | ||
| std::memcpy( stateTemp.data(), StateOld_vec.data(), nState * sizeof( double ) ); | ||
| respPert.stateVars = stateTemp.data(); | ||
|
|
||
| AlgorithmicModuli< 3 > tanPert; // dummy | ||
|
|
||
| computeStress( respPert, tanPert, deformation, timeIncrement ); | ||
|
|
||
| // Return Stress (flattened 9x1) | ||
| Map< VectorXd > res( respPert.tau.data(), 9 ); | ||
| return res; | ||
| }; | ||
|
|
||
| sensitivities.dStress_dStateOld = forwardDifference( func_dStress_dStateOld, StateOld_map ); | ||
| if ( sensitivities.dStress_dStateOld.hasNaN() ) { | ||
| std::cout << "Computed dStress/dStateOld has NaN!" << std::endl; | ||
| std::cout << "F: \n" << deformation.F << std::endl; | ||
| std::cout << "StateOld: \n" << Map< const VectorXd >( stateOld.data(), nState ) << std::endl; | ||
| } |
There was a problem hiding this comment.
Multiple debug print statements are left commented out in production code (lines 110-111, 155-163, 221-227, 255-260, 278-282). While this is a WIP PR, these should either be removed or replaced with a proper logging framework before merging to maintain code cleanliness.
modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h
Outdated
Show resolved
Hide resolved
modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp
Outdated
Show resolved
Hide resolved
| virtual void computeStressWithSensitivities( ConstitutiveResponse< 3 >& response, | ||
| AlgorithmicModuli< 3 >& tangents, | ||
| StateSensitivities& sensitivities, | ||
| const Deformation< 3 >& deformation, | ||
| const TimeIncrement& timeIncrement ) const |
There was a problem hiding this comment.
The computeStressWithSensitivities method modifies response.stateVars as a side effect while also computing sensitivities. This could be confusing for users who might expect the method to be const (like computeStress) or who might not expect state updates when computing sensitivities. Consider documenting this behavior clearly in the method documentation, explicitly stating that the method updates the state variables in addition to computing sensitivities.
...s/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h
Outdated
Show resolved
Hide resolved
...s/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h
Outdated
Show resolved
Hide resolved
...s/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h
Outdated
Show resolved
Hide resolved
examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/example.cpp
Outdated
Show resolved
Hide resolved
…rmotMaterialFiniteStrainSubstepped.h Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…rmotMaterialFiniteStrain.h Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
….cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…rmotMaterialFiniteStrainSubstepped.h Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…epped/example.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…rmotMaterialFiniteStrainSubstepped.h Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…rmotMaterialFiniteStrainSubstepped.h Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…epped/example.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/Makefile
Outdated
Show resolved
Hide resolved
| /** | ||
| * @struct StateSensitivities | ||
| * @brief Sensitivity matrices required for analytical substepping. | ||
| * All matrices represent derivatives of flattened arrays (Row-Major Fastor layout). |
There was a problem hiding this comment.
Attention: Eigen itself is not Row-Major;
No description provided.