Material Point Solver for Finite Strain Materials

Overview

The Material Point Solver (Finite Strain) provides a standalone driver for simulating the response of a single material point using finite strain constitutive models from the MAteRialMOdellingToolbox (marmot) framework.

Unlike the hypo-elastic solver, this class handles large deformations using full \(3\times3\) tensor kinematics. It solves for the Displacement Gradient (\(\nabla \mathbf{u}\)) and the Kirchhoff Stress (\(\boldsymbol{\tau}\)).

It is designed for:

• Testing and verifying finite strain material models (e.g., hyperelasticity, finite plasticity),

• Performing mixed displacement-gradient/stress-controlled loading paths,

• Recording full tensor histories including material tangents and internal state variables.

The solver automatically handles time stepping, adaptive cutbacks upon non-convergence, and output formatting.

Main Features

• Finite Strain Kinematics: Solves for \(\mathbf{F} = \mathbf{I} + \nabla \mathbf{u}\).

• Mixed Control: Supports independent control of displacement gradient or stress for each of the 9 tensor components.

• Adaptive Time Stepping: Automatically reduces time step size (\(\Delta t\)) if convergence fails.

• Newton–Raphson Solver: robust iterative solver with configurable residual and correction tolerances.

• History Recording: Stores time, stress, deformation gradient, tangent moduli, and state variables.

• CSV Export: Exports simulation history for post-processing.

Solver Structure

The solver is organized around a few key data structures:

• SolverOptions — Defines numerical tolerances (residual/correction) and iteration limits.

• Step — Represents a loading phase. It defines target increments for \(\nabla \mathbf{u}\) and \(\boldsymbol{\tau}\), along with boolean flags determining which variable controls which tensor component.

• Increment — A sub-step automatically generated within each Step based on time discretization.

• HistoryEntry — Records the converged state (Time, \(\boldsymbol{\tau}\), \(\mathbf{F}\), \(\mathbb{C}\), state variables) at specific points.

Usage Workflow

A typical workflow consists of:

1. Initialize the solver with a material model name and properties.

2. Set the initial state (initial stress and internal variables).

3. Define one or more loading Step objects with control flags and targets.

4. Add the steps to the solver.

5. Solve the problem.

6. Inspect or export the resulting history.

Example Usage

Here is a simple example demonstrating how to set up and run the finite strain material point solver using a compressible Neo-Hookean material model. It can be found in the examples/finite-strain-mp-solver/CompressibleNeoHooke directory of the Marmot repository together with the appropriate CMakeLists.txt.

Finite Strain Solver Example

#include "Marmot/MarmotFastorTensorBasics.h"
#include "Marmot/MarmotMaterialPointSolverFiniteStrain.h"

int main()
{
  using namespace Marmot::Solvers;
  using namespace Marmot::FastorStandardTensors;

  // 1) Define the material model
  std::string materialName = "COMPRESSIBLENEOHOOKE";
  double      properties[] = { 3500, 1500 }; // Example: Bulk and Shear moduli
  int         nProps       = 2;

  // 2) Configure solver options
  MarmotMaterialPointSolverFiniteStrain::SolverOptions options;
  options.maxIterations       = 25;
  options.residualTolerance   = 1e-10;
  options.correctionTolerance = 1e-10;

  // 3) Create solver instance
  MarmotMaterialPointSolverFiniteStrain solver( materialName, properties, nProps, options );

  // 4) Set the initial state
  int             nSV       = solver.getNumberOfStateVariables();
  Eigen::VectorXd initialSV = Eigen::VectorXd::Zero( nSV );

  // Initial stress is zero
  solver.setInitialState( Tensor33d( 0.0 ), initialSV );

  // 5) Define a loading step (Uniaxial extension in 11-direction)
  MarmotMaterialPointSolverFiniteStrain::Step step;
  step.timeStart = 0.0;
  step.timeEnd   = 1.0;
  step.dTStart   = 0.1;
  step.dTMin     = 1e-6;
  step.dTMax     = 0.5;

  // Initialize targets and flags
  step.gradUIncrementTarget        = Tensor33d( 0.0 );
  step.stressIncrementTarget       = Tensor33d( 0.0 );
  step.isGradUComponentControlled  = Tensor33t< bool >( false );
  step.isStressComponentControlled = Tensor33t< bool >( true );

  // Setup Control:
  // - Control gradU_11 (stretch to 50%)
  step.gradUIncrementTarget( 0, 0 )       = 0.5;
  step.isGradUComponentControlled( 0, 0 ) = true;

  // - Unset control on tau_11
  step.isStressComponentControlled( 0, 0 ) = false;

  // 6) Run the solver
  solver.addStep( step );
  solver.solve();

  // 7) Output the results
  solver.exportHistoryToCSV( "finite_strain_history.csv" );

  return 0;
}

To compile and run the example, you can use cmake with the configuration CMakeLists.txt. Ensure that the paths to Marmot, Eigen, and Fastor match your system configuration.

cmake configuration for Finite Strain Material Point Solver Example

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

Typical Output

Running the solver prints detailed convergence information to the console:

Solving step from 0 to 1
+------------------------------------------------------------------------------+
  Solving increment 1, time: 0 to 0.1, dT: 0.1
    Iteration 0, ||ddGradU||: 0.00e+00, ||R||: 5.00e-02
    Iteration 1, ||ddGradU||: 1.25e-03, ||R||: 4.12e-05
    Iteration 2, ||ddGradU||: 2.10e-07, ||R||: 1.10e-11
    Converged after 3 iterations.
  Material state for time: 1.000000e-01
  tau:
   [2.500000e+01, 0.000000e+00, 0.000000e+00
    0.000000e+00, 0.000000e+00, 0.000000e+00
    0.000000e+00, 0.000000e+00, 0.000000e+00]
  ...

Guidelines and Notes

• Tensors: The solver uses the Fastor library. Inputs usually expect Tensor33d (3x3 double) or Tensor9d (flattened 9x1 double).

• Control Flags: For every component \((i,j)\) of the \(3\times3\) tensor, exactly one of isGradUComponentControlled(i,j) or isStressComponentControlled(i,j) must be set to true. The Step::checkControl() method enforces this.

• Singularity Handling: When using mixed control, the tangent matrix is modified to ensure invertibility. Rows corresponding to displacement-controlled components are zeroed and the diagonal set to 1.

• Stress Measure: The solver works with Kirchhoff stress (\(\boldsymbol{\tau}\)). If the material model computes Cauchy stress, it must be converted internally or by the model.

CSV Export

The history can be exported to CSV. The columns are formatted as follows:

• Time: Simulation time.

• tau[9]: Kirchhoff stress components (11, 12, 13, 21, 22, 23, 31, 32, 33).

• F[9]: Deformation gradient components (11, 12, 13, 21, 22, 23, 31, 32, 33).

• SV[N]: Internal state variables (SV1, SV2, …).

class MarmotMaterialPointSolverFiniteStrain

Collaboration diagram for Marmot::Solvers::MarmotMaterialPointSolverFiniteStrain:

Solver for material point problems with finite strain materials.

This class implements a solver for material point problems using finite strain material models. It supports loading steps with controlled displacement gradient and stress components, adaptive time stepping, and history recording.

Public Functions

MarmotMaterialPointSolverFiniteStrain(const std::string &materialName, const double *materialProperties, const int nMaterialProperties, const SolverOptions &options)

Constructor for the MarmotMaterialPointSolverFiniteStrain class.

Parameters:

• materialName – Name of the finite strain material model

• materialProperties – Array of material properties

• nMaterialProperties – Number of material properties

• options – Solver options controlling nonlinear iteration settings

void addStep(const Step &step)

Add a loading step to the solver.

Parameters:

step – The Step to be added

inline std::vector<Step> getSteps() const

Get the list of added loading steps.

Returns:

A vector of Step containing the added steps

inline void clearSteps()

Clear all added loading steps.

void setInitialState(const FastorStandardTensors::Tensor33d &initialStress, const Eigen::VectorXd &initialStateVars)

Set the initial state of the material model.

Parameters:

• initialStress – The initial stress in Voigt notation

• initialStateVars – The initial state variables

inline int getNumberOfStateVariables() const

Get the number of state variables in the material model.

Returns:

The number of state variables

void resetToInitialState()

Reset the solver to the initial state.

This function resets the initial stress and state variables of the material model.

void solve()

Solve the material point problem for all added steps.

inline const std::vector<HistoryEntry> &getHistory() const

Get the recorded history of the simulation.

Returns:

A vector of HistoryEntry containing the recorded history

inline void clearHistory()

Clear the recorded history.

void printHistory()

Print the recorded history to the console.

void exportHistoryToCSV(const std::string &filename)

Export the recorded history to a CSV file.

Parameters:

filename – The name of the CSV file to export to

Private Functions

void solveStep(const Step &step)

Solve a single loading step.

This function iterates over the increments defined in the step, calling solveIncrement for each increment. It manages time stepping and ensures that the entire step is covered.

Parameters:

step – The Step to be solved

void solveIncrement(const Increment &increment)

Solve a single increment within a loading step.

This function implements a Newton-Raphson iterative solver to compute the stress and deformation state for the given increment. It updates the material state variables and records the history after convergence.

Parameters:

increment – The Increment to be solved

Throws:

std::runtime_error – if the solver does not converge

FastorStandardTensors::Tensor9d computeResidual(const FastorStandardTensors::Tensor9d &stressIncrement, const FastorStandardTensors::Tensor9d &target, const Increment &increment)

Compute the residual for the current increment.

This function calculates the residual vector based on the difference between the computed stress increment and the target increment, taking into account which components are controlled by strain or stress.

Parameters:

• stressIncrement – The computed stress increment

• target – The target (mixed) stress/strain increment

• increment – The Increment containing control information

Returns:

The computed residual vector

void modifyTangent(FastorStandardTensors::Tensor99d &tangent, const Increment &increment)

Modify the material tangent matrix based on control type.

This function adjusts the tangent matrix to account for the components that are controlled by displacement gradient or stress, ensuring that the solver correctly handles mixed control scenarios. This is done by zeroing out rows corresponding to displacement gradient controlled components and setting their diagonal entries to one.

Parameters:

• tangent – The material tangent matrix to be modified

• increment – The Increment containing control information

Private Members

std::unique_ptr<MarmotMaterialFiniteStrain> material

The finite strain material model.

int nStateVars

Number of state variables in the material model.

Eigen::VectorXd stateVars

Current state variables.

Eigen::VectorXd _initialStateVars

Initial state variables.

Eigen::VectorXd stateVarsTemp

Temporary state variables for computations.

FastorStandardTensors::Tensor33d stress = FastorStandardTensors::Tensor33d(0.0)

The Kirchhoff stress.

FastorStandardTensors::Tensor33d _initialStress = FastorStandardTensors::Tensor33d(0.0)

The initial Kirchhoff stress.

FastorStandardTensors::Tensor33d gradU = FastorStandardTensors::Tensor33d(0.0)

The displacement gradient.

FastorStandardTensors::Tensor3333d dTau_dF = FastorStandardTensors::Tensor3333d(0.0)

The material tangent dTau/dF.

std::vector<Step> steps

List of loading steps.

std::vector<HistoryEntry> history = std::vector<HistoryEntry>()

History of the simulation.

const SolverOptions options

Solver options.

struct HistoryEntry

Struct to record the history of the simulation.

Each entry contains time, stress, deformation, and state variables.

Public Functions

inline void print() const

Print the history entry to the console.

inline HistoryEntry(double time, const FastorStandardTensors::Tensor33d &stress, const FastorStandardTensors::Tensor33d &F, const FastorStandardTensors::Tensor3333d &dTau_dF, const Eigen::VectorXd &stateVars)

Constructor for HistoryEntry.

Parameters:

• time – Time at the history entry

• stress – Stress tensor at the history entry

• F – Deformation gradient at the history entry

• dTau_dF – Material tangent at the history entry

• stateVars – State variables at the history entry

inline HistoryEntry()

Public Members

double time

Time at the history entry.

FastorStandardTensors::Tensor33d stress

Stress at the history entry.

FastorStandardTensors::Tensor33d F

deformation gradient at the history entry

FastorStandardTensors::Tensor3333d dTau_dF

Material tangent at the history entry.

Eigen::VectorXd stateVars

State variables at the history entry.

struct Increment

Struct to define a loading increment.

Each increment contains displacement gradient and stress increments, control flags, time information, and iteration limits.

Public Members

FastorStandardTensors::Tensor9d gradUIncrement

Target displacement gradient increment for the step.

FastorStandardTensors::Tensor9d stressIncrement

Target Kirchhoff stress increment for the step.

FastorStandardTensors::Tensor9t<bool> isGradUComponentControlled

Flags to indicate which displecement gradient components are controlled

FastorStandardTensors::Tensor9t<bool> isStressComponentControlled

Flags to indicate which stress components are controlled.

double timeOld

Old time at the beginning of the increment.

double dT

Time step size for the increment.

struct SolverOptions

Struct to define solver options.

Contains parameters for controlling the solver’s behavior.

Public Members

int maxIterations = 25

Maximum number of iterations per increment (default: 25)

double residualTolerance = 1e-10

Convergence tolerance (default: 1e-10)

double correctionTolerance = 1e-10

Correction tolerance (default: 1e-10)

struct Step

Struct to define a loading step.

Each step contains targets for displacement gradient and stress states, time information and time step control parameters.

Public Functions

inline void checkControl() const

Check that for each component, either deformation or stress is controlled.

Throws:

std::runtime_error – if the condition is not met

Public Members

FastorStandardTensors::Tensor33d gradUIncrementTarget

Target displacement gradient increment for the step.

FastorStandardTensors::Tensor33d stressIncrementTarget

Target Kirchhoff stress increment for the step.

FastorStandardTensors::Tensor33t<bool> isGradUComponentControlled

Flags to indicate which displacement gradient components are controlled

FastorStandardTensors::Tensor33t<bool> isStressComponentControlled

Flags to indicate which stress components are controlled.

double timeStart = 0.0

Start time of the step.

double timeEnd = 1.0

End time of the step.

double dTStart = 0.1

Initial time step size.

double dTMin = 1e-6

Minimum time step size.

double dTMax = 0.5

Maximum time step size.

int maxIncrements = 100

Maximum number of increments in the step.