diff --git a/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/CMakeLists.txt b/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/CMakeLists.txt new file mode 100644 index 00000000..a1c437d0 --- /dev/null +++ b/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/CMakeLists.txt @@ -0,0 +1,37 @@ +cmake_minimum_required(VERSION 3.16) +project(MarmotSolverExample LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +find_package(Python3 COMPONENTS Interpreter REQUIRED) + +execute_process( + COMMAND "${Python3_EXECUTABLE}" -c "import sys; print(sys.prefix, end='')" + OUTPUT_VARIABLE PYTHON_PREFIX + OUTPUT_STRIP_TRAILING_WHITESPACE +) + +set(MARMOT_ROOT "${PYTHON_PREFIX}") +set(EIGEN_DIR "${PYTHON_PREFIX}/include/eigen3") +set(FASTOR_DIR "${PYTHON_PREFIX}/include/Fastor") + +add_executable(solver_example example.cpp) + +target_include_directories(solver_example PRIVATE + "${MARMOT_ROOT}/include" + "${EIGEN_DIR}" + "${FASTOR_DIR}" +) + +target_link_directories(solver_example PRIVATE "${MARMOT_ROOT}/lib") + +target_link_libraries(solver_example PRIVATE Marmot) + +target_compile_options(solver_example PRIVATE -Wall -Wextra -fPIC) + +set_target_properties(solver_example PROPERTIES + BUILD_RPATH "${MARMOT_ROOT}/lib" + INSTALL_RPATH "${MARMOT_ROOT}/lib" +) diff --git a/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/example.cpp b/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/example.cpp new file mode 100644 index 00000000..e1707bdf --- /dev/null +++ b/examples/finite-strain-mp-solver/FiniteStrainJ2PlasticitySubstepped/example.cpp @@ -0,0 +1,94 @@ +#include "Marmot/MarmotFastorTensorBasics.h" +#include "Marmot/MarmotMaterialPointSolverFiniteStrain.h" +// Ensure the registration unit is linked or included if using static registration +// #include "Marmot/MaterialRegistration.cpp" + +int main() +{ + using namespace Marmot::Solvers; + using namespace Marmot::FastorStandardTensors; + + // 1) Define the material model + // Must match the key used in MarmotMaterialFiniteStrainFactory::registerMaterial + std::string materialName = "FINITESTRAINJ2PLASTICITY_SUBSTEPPED"; + + // Define Properties + // [0]: nSubsteps (Specific to SubsteppingMaterial) + // [1...8]: Standard FiniteStrainJ2Plasticity properties + double properties[] = { + 20, // nSubsteps (20 substeps per global increment) + 35, // K (Bulk Modulus) + 15, // G (Shear Modulus) + 0.45, // fy (Initial Yield Stress) + 0.75, // fyInf (Saturated Yield Stress) + 0.5, // eta (Saturation rate) + 0.1, // H (Linear Hardening) + 1.0, // implementationType (1 = Full Analytical Return Mapping) + 0.0 // density + }; + int nProps = 9; + + // 2) Configure solver options + // With a consistent analytical tangent from the substepper, + // we expect quadratic convergence (typically < 5 iterations). + MarmotMaterialPointSolverFiniteStrain::SolverOptions options; + options.maxIterations = 15; + options.residualTolerance = 1e-8; + options.correctionTolerance = 1e-8; + + // 3) Create solver instance + // The factory will instantiate SubsteppingMaterial, which internally instantiates FiniteStrainJ2Plasticity + MarmotMaterialPointSolverFiniteStrain solver( materialName, properties, nProps, options ); + + // 5) Define Loading Steps + + // --- Step 1: Loading (Uniaxial extension to 10%) --- + MarmotMaterialPointSolverFiniteStrain::Step step1; + step1.timeStart = 0.0; + step1.timeEnd = 1.0; + step1.dTStart = 1e-0; // 10 global increments + step1.dTMin = 1e-1; + step1.dTMax = 0.5; + + step1.gradUIncrementTarget = Tensor33d( 0.0 ); + step1.stressIncrementTarget = Tensor33d( 0.0 ); + step1.isGradUComponentControlled = Tensor33t< bool >( false ); + step1.isStressComponentControlled = Tensor33t< bool >( true ); + + // Control gradU_11 (stretch to 10%) + step1.gradUIncrementTarget( 0, 0 ) = 1; + step1.isGradUComponentControlled( 0, 0 ) = true; + step1.isStressComponentControlled( 0, 0 ) = false; + + // Mixed control for lateral contraction (Plane Stress approximation via Solver) + // We keep stress 22 and 33 controlled to 0.0 (default true) + + solver.addStep( step1 ); + + // --- Step 2: Unloading (Return to 0% global strain) --- + // This validates that F_n is correctly updated and stored between steps + MarmotMaterialPointSolverFiniteStrain::Step step2 = step1; + step2.timeStart = 1.0; + step2.timeEnd = 2.0; + + // Reverse direction: target increment is -0.10 to go back to 0 + step2.gradUIncrementTarget( 0, 0 ) = -1; + + solver.addStep( step2 ); + + // 6) Run the solver + try { + solver.solve(); + } + catch ( const std::exception& e ) { + // Catch exceptions thrown by the substepper (e.g., if inner Newton fails) + std::cerr << "\n!!! SOLVER FAILED !!!\nError: " << e.what() << std::endl; + return 1; + } + + // 7) Output the results + // The CSV will contain extra columns for the substepper state (Substepping_F_n components) + solver.exportHistoryToCSV( "j2_substepped_results.csv" ); + + return 0; +} diff --git a/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h index c20c2ea9..2a84ffce 100644 --- a/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h +++ b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h @@ -125,6 +125,7 @@ class MarmotMaterialFiniteStrain { AlgorithmicModuli< 3 >& tangents, const Deformation< 3 >&, const TimeIncrement& ) const = 0; + /** * @brief Computes the Kirchhoff stress given the deformation, time increment, and eigen deformation. * @param[inout] response ConstitutiveResponse instance diff --git a/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h new file mode 100644 index 00000000..df6cb55d --- /dev/null +++ b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrainSubstepped.h @@ -0,0 +1,288 @@ +/* --------------------------------------------------------------------- + * _ + * _ __ ___ __ _ _ __ _ __ ___ ___ | |_ + * | '_ ` _ \ / _` | '__| '_ ` _ \ / _ \| __| + * | | | | | | (_| | | | | | | | | (_) | |_ + * |_| |_| |_|\__,_|_| |_| |_| |_|\___/ \__| + * + * Unit of Strength of Materials and Structural Analysis + * University of Innsbruck, + * 2020 - today + * + * festigkeitslehre@uibk.ac.at + * + * Alexander Dummer alexander.dummer@uibk.ac.at + * + * This file is part of the MAteRialMOdellingToolbox (marmot). + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of marmot. + * --------------------------------------------------------------------- + */ + +#pragma once + +#include "Marmot/MarmotFastorTensorBasics.h" +#include "Marmot/MarmotMaterialFiniteStrain.h" +#include "Marmot/MarmotNumericalDifferentiation.h" +#include +#include +#include +#include + +namespace Marmot::Materials { + + using namespace Fastor; + + /** + * @class MarmotMaterialFiniteStrainSubstepped + * @brief A decorator that applies time substepping with analytical tangent accumulation. + * @tparam BaseMaterialType The concrete material class to wrap. + * Properties: [nSubsteps, Prop1, Prop2, ...] + */ + template < typename BaseMaterialType > + class MarmotMaterialFiniteStrainSubstepped : public MarmotMaterialFiniteStrain { + protected: + std::unique_ptr< BaseMaterialType > baseMaterial; + int nSubsteps; + + public: + /** + * @struct StateSensitivities + * @brief Sensitivity matrices required for analytical substepping. + * All matrices represent derivatives of flattened arrays and are stored in Eigen's column-major format. + */ + struct StateSensitivities { + /** @brief Jacobian of State_new w.r.t Deformation Gradient F_new. + * Size: nStateVars x 9 + */ + Eigen::MatrixXd dState_dF; + + /** @brief Jacobian of State_new w.r.t State_old. + * Size: nStateVars x nStateVars + */ + Eigen::MatrixXd dState_dStateOld; + + /** @brief Jacobian of Stress w.r.t State_old. + * Size: 9 x nStateVars + */ + Eigen::MatrixXd dStress_dStateOld; + }; + + MarmotMaterialFiniteStrainSubstepped( const double* matProperties_, int nMaterialProperties_, int materialNumber_ ) + : MarmotMaterialFiniteStrain( matProperties_, nMaterialProperties_, materialNumber_ ) + { + // 1. Parse Substepping parameters + if ( nMaterialProperties_ < 1 ) { + throw std::invalid_argument( "MarmotMaterialFiniteStrainSubstepped requires at least 1 property (nSubsteps)" ); + } + this->nSubsteps = static_cast< int >( matProperties_[0] ); + if ( this->nSubsteps < 1 ) + this->nSubsteps = 1; + + baseMaterial = std::make_unique< BaseMaterialType >( matProperties_ + 1, + nMaterialProperties_ - 1, + materialNumber_ ); + + initializeStateLayout(); + } + + virtual ~MarmotMaterialFiniteStrainSubstepped() = default; + + void initializeStateLayout() override + { + + // Add our state variable: Deformation Gradient at start of GLOBAL step (F_n) + this->stateLayout.add( "Substepping_F_n", 9 ); + + // Add Base Material state variables + this->stateLayout.add( "materialstate", baseMaterial->getNumberOfRequiredStateVars() ); + + this->stateLayout.finalize(); + } + + void initializeYourself( double* stateVars, int nStateVars ) override + { + // 1. Initialize Base Material part + int baseVarsCount = baseMaterial->getNumberOfRequiredStateVars(); + + baseMaterial->initializeYourself( stateLayout.getPtr( stateVars, "materialstate" ), baseVarsCount ); + + FastorStandardTensors::Tensor33d& + Fn = this->stateLayout.getAs< FastorStandardTensors::Tensor33d& >( stateVars, "Substepping_F_n" ); + Fn.eye(); + } + + /** + * @brief Extended stress update computing sensitivities for substepping. + * * The default implementation uses Forward Finite Differences. + * @param response Constitutive response to be computed. + * @param tangents Algorithmic moduli to be computed. + * @param sensitivities Sensitivity matrices to be computed. + * @param deformation Current deformation state. + * @param timeIncrement Time increment information. + * + * @note This function can be overridden in derived classes to provide analytical sensitivities. + * @note The base material's computeStress function is called within this function. + */ + virtual void computeStressWithSensitivities( ConstitutiveResponse< 3 >& response, + AlgorithmicModuli< 3 >& tangents, + StateSensitivities& sensitivities, + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement ) const + { + using namespace Eigen; + using namespace Marmot::NumericalAlgorithms::Differentiation; + + int nState = baseMaterial->getNumberOfRequiredStateVars(); + + std::vector< double > stateOld( nState ); + std::memcpy( stateOld.data(), response.stateVars, nState * sizeof( double ) ); + + baseMaterial->computeStress( response, tangents, deformation, timeIncrement ); + + if ( nState == 0 ) + return; + + sensitivities.dState_dF.resize( nState, 9 ); + sensitivities.dState_dStateOld.resize( nState, nState ); + sensitivities.dStress_dStateOld.resize( 9, nState ); + + 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; + respPert.stateVars = stateTemp.data(); + + AlgorithmicModuli< 3 > tanPert; + + baseMaterial->computeStress( respPert, tanPert, defPerturbed, timeIncrement ); + + Map< VectorXd > res( respPert.stateVars, nState ); + return res; + }; + + auto func_dStressAndState_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 + baseMaterial->computeStress( respPert, tanPert, deformation, timeIncrement ); + + // return [stress; state_new] as a single vector + VectorXd res( 9 + nState ); + // Fill Stress part + std::memcpy( res.data(), respPert.tau.data(), 9 * sizeof( double ) ); + // Fill State_new part + std::memcpy( res.data() + 9, respPert.stateVars, nState * sizeof( double ) ); + + return res; + }; + + // Compute dState/dF_new + sensitivities.dState_dF = forwardDifference( func_dState_dF, Map< const VectorXd >( deformation.F.data(), 9 ) ); + + // Compute [dStress/dStateOld; dState/dStateOld] + const MatrixXd dStressAndState_dStateOld = forwardDifference( func_dStressAndState_dStateOld, + Map< const VectorXd >( stateOld.data(), nState ) ); + + // Extract dState/dStateOld + sensitivities.dState_dStateOld = dStressAndState_dStateOld.block( 9, 0, nState, nState ); + + // Extract dStress/dStateOld + sensitivities.dStress_dStateOld = dStressAndState_dStateOld.block( 0, 0, 9, nState ); + } + + /** + * @brief Compute stress with time substepping and analytical tangent accumulation. + * @param response Constitutive response to be computed. + * @param tangents Algorithmic moduli to be computed. + * @param deformation Current deformation state. + * @param timeIncrement Time increment information. + * + * @details The deformation increment is divided into nSubsteps sub-increments. + * For each sub-increment, the base material's stress update is called, + * and the sensitivities are used to accumulate the overall tangent. + */ + void computeStress( ConstitutiveResponse< 3 >& response, + AlgorithmicModuli< 3 >& tangents, + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement ) const override + { + using namespace Eigen; + using namespace FastorStandardTensors; + + Tensor33d& Fn_ref = this->stateLayout.getAs< Tensor33d& >( response.stateVars, "Substepping_F_n" ); + const Tensor33d Fn = Fn_ref; + const Tensor33d Fn1 = deformation.F; + + int nBaseState = baseMaterial->getNumberOfRequiredStateVars(); + + // J_accum = d(State_i)/d(F_global). Size: nState x 9 + MatrixXd J_accum = MatrixXd::Zero( nBaseState, 9 ); + + double dt_sub = timeIncrement.dT / static_cast< double >( nSubsteps ); + double t_curr = timeIncrement.time; + + ConstitutiveResponse< 3 > subResponse; + subResponse.stateVars = this->stateLayout.getPtr( response.stateVars, "materialstate" ); + + AlgorithmicModuli< 3 > subTangents; + StateSensitivities subSensitivities; + Deformation< 3 > subDef; + + for ( int i = 1; i <= nSubsteps; ++i ) { + double xi = static_cast< double >( i ) / static_cast< double >( nSubsteps ); + + // Linear Interpolation: F(xi) = (1-xi)*F_n + xi*F_n1 + subDef.F = ( 1.0 - xi ) * Fn + xi * Fn1; + const double dFsub_dFglobal = xi; + + TimeIncrement subTime = { t_curr, dt_sub }; + + computeStressWithSensitivities( subResponse, subTangents, subSensitivities, subDef, subTime ); + + // Update accumulated Jacobian for history + if ( i < nSubsteps ) { + if ( nBaseState > 0 ) { + J_accum = subSensitivities.dState_dF * dFsub_dFglobal + subSensitivities.dState_dStateOld * J_accum; + } + } + else { + // Final Step: Compute overall Tangent dTau/dF_global + MatrixXd C_hist = MatrixXd::Zero( 9, 9 ); + if ( nBaseState > 0 ) { + C_hist = subSensitivities.dStress_dStateOld * J_accum; + } + + // build global tangent + // transpose because second term is from Eigen's world with column-major layout + Fastor::Tensor< double, 9, 9 > C_global = subTangents.dTau_dF * dFsub_dFglobal + + transpose( Fastor::Tensor< double, 9, 9 >( C_hist.data() ) ); + + std::memcpy( tangents.dTau_dF.data(), C_global.data(), 81 * sizeof( double ) ); + } + t_curr += dt_sub; + } + + // set final response + response.tau = subResponse.tau; + response.elasticEnergyDensity = subResponse.elasticEnergyDensity; + + // Update stored F_n to F_n1 + Fn_ref = deformation.F; + } + }; + +} // namespace Marmot::Materials diff --git a/modules/core/MarmotMathCore/module.cmake b/modules/core/MarmotMathCore/module.cmake index 0cf8f44e..fea13995 100644 --- a/modules/core/MarmotMathCore/module.cmake +++ b/modules/core/MarmotMathCore/module.cmake @@ -6,6 +6,7 @@ list(APPEND publicheaders "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotMath.h" "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotConstants.h" "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotFastorTensorBasics.h" + "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotNumericalDifferentiation.h" "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotTensor.h" "${CMAKE_CURRENT_LIST_DIR}/include/Marmot/MarmotTypedefs.h" ) diff --git a/modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp b/modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp index c66446d8..621a16fa 100644 --- a/modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp +++ b/modules/core/MarmotMathCore/src/MarmotNumericalDifferentiation.cpp @@ -48,33 +48,40 @@ namespace Marmot { MatrixXd centralDifference( const vector_to_vector_function_type& F, const VectorXd& X ) { - const auto xSize = X.rows(); - MatrixXd J( xSize, xSize ); + + MatrixXd J; VectorXd leftX( xSize ); VectorXd rightX( xSize ); for ( auto i = 0; i < xSize; i++ ) { double volatile h = std::max( 1.0, std::abs( X( i ) ) ) * Marmot::Constants::CubicRootEps; - // clang-format off + leftX = X; rightX = X; leftX( i ) -= h; rightX( i ) += h; - J.col( i ) = ( F( rightX ) - F( leftX ) ) - / //------------------------------------ - ( 2. * h ); + VectorXd yRight = F( rightX ); + VectorXd yLeft = F( leftX ); + + if ( i == 0 ) { + J.resize( yRight.size(), xSize ); + } + + // clang-format off + J.col( i ) = ( yRight - yLeft ) + / //----------------------- + ( 2. * h ); // clang-format on } return J; } - namespace Complex { /* - * Implementation of Numerical Differantiation using Complex Step Approximations + * Implementation of Numerical Differentiation using Complex Step Approximations * * further Information can be found in * - Martins et al. (2003) The Complex-Step Derivative Approximation diff --git a/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2Plasticity.cpp b/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2Plasticity.cpp index db0e53c7..27505fa8 100644 --- a/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2Plasticity.cpp +++ b/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2Plasticity.cpp @@ -684,6 +684,6 @@ namespace Marmot::Materials { } Tensor33d& Fp = stateLayout.getAs< Tensor33d& >( stateVars, "Fp" ); - Fp.eye(); + memcpy( Fp.data(), Spatial3D::I.data(), 9 * sizeof( double ) ); } } // namespace Marmot::Materials diff --git a/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2PlasticityRegistration.cpp b/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2PlasticityRegistration.cpp index 992fc5fa..c934f906 100644 --- a/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2PlasticityRegistration.cpp +++ b/modules/materials/FiniteStrainJ2Plasticity/src/FiniteStrainJ2PlasticityRegistration.cpp @@ -1,5 +1,6 @@ #include "Marmot/FiniteStrainJ2Plasticity.h" #include "Marmot/MarmotMaterialFiniteStrainFactory.h" +#include "Marmot/MarmotMaterialFiniteStrainSubstepped.h" namespace Marmot::Materials { @@ -9,5 +10,13 @@ namespace Marmot::Materials { const static bool FiniteStrainJ2PlasticityRegistered = MarmotMaterialFiniteStrainFactory::registerMaterial< FiniteStrainJ2Plasticity >( "FINITESTRAINJ2PLASTICITY" ); + + // Register the SUBSTEPPED J2 model + // This allows you to use "FINITESTRAINJ2PLASTICITY_SUBSTEPPED" in your input file. + // The properties array must start with [nSubsteps, K, G, fy, ...] + const static bool FiniteStrainJ2PlasticitySubsteppedRegistered = MarmotMaterialFiniteStrainFactory:: + registerMaterial< MarmotMaterialFiniteStrainSubstepped< FiniteStrainJ2Plasticity > >( + "FINITESTRAINJ2PLASTICITY_SUBSTEPPED" ); + } // namespace Registration } // namespace Marmot::Materials