AmplePy/ample.py at master · QuackFrog/AmplePy · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
"""
AmplePy - Unofficial Python port of AMPLE: A Material Point Learning Environment
--------------------------------------------------------------------------------

Original Author: William Coombs
Date:   27/08/2020
Ported by: Chen, Yanbin
Date:   08/09/2025
Description:
Large deformation elasto-plastic (EP) material point method (MPM) code
based on an updated Lagrangian (UL) descripition of motion with a
quadrilateral background mesh.

--------------------------------------------------------------------------------
"""

import numpy as np
from scipy.sparse import csr_matrix

# Import core functions
from functions.detMPs import detMPs
from functions.linSolve import linSolve
from functions.detExtForce import detExtForce
from functions.detFDoFs import detFDoFs
from functions.elemMPinfo import elemMPinfo
from functions.updateMPs import updateMPs

# Import setup functions
from setup import setupGrid, setupGrid_beam, setupGrid_collapse

# Import plotting functions
from plotting.postPro import postPro


def main(problem_type='default'):
    """
    Main function to run the AmplePy material point analysis

    Parameters:
    -----------
    problem_type : str
        Type of problem to solve ('default', 'beam', 'collapse')
    """

    # Setup problem based on type
    if problem_type == 'beam':
        lstps, g, mpData, mesh = setupGrid_beam()
    elif problem_type == 'collapse':
        lstps, g, mpData, mesh = setupGrid_collapse()
    else:
        lstps, g, mpData, mesh = setupGrid()

    NRitMax = 10
    tol = 1e-9  # Newton Raphson parameters
    nodes, nD = mesh['coord'].shape  # number of nodes and dimensions
    nels, nen = mesh['etpl'].shape  # number of elements and nodes/element
    nDoF = nodes * nD  # total number of degrees of freedom
    nmp = len(mpData)  # number of material points
    lstp = 0  # zero loadstep counter (for plotting function)
    uvw = np.zeros(nDoF)  # zeros displacements (for plotting function)

    print('AmplePy: A Material Point Learning Environment (Python Version)')
    print('-------------------------------------------------------------')
    print('Problem type:', problem_type)
    print(f'Mesh info: {nodes} nodes, {nD} dimensions, {nels} elements')
    print(f'Material points: {nmp}')
    print(f'Degrees of freedom: {nDoF}')

    # Plotting initial state & mesh
    postPro(lstp, mesh, mpData, problem_type)

    for lstp in range(1, lstps + 1):  # loadstep loop
        # text output to screen (loadstep)
        print(f'\nloadstep {lstp:4d} of {lstps:4d}')

        # Material point - element information
        mesh, mpData = elemMPinfo(mesh, mpData)

        # External force calculation (total)
        fext = detExtForce(nodes, nD, g, mpData)

        # Current external force value
        fext = fext * lstp / lstps

        oobf = fext.copy()  # initial out-of-balance force
        fErr = 1  # initial error
        frct = np.zeros(nDoF)  # zero the reaction forces
        uvw = np.zeros(nDoF)  # zero the displacements
        fd = detFDoFs(mesh)  # free degrees of freedom
        NRit = 0  # zero the iteration counter
        Kt = csr_matrix((nDoF, nDoF))  # zero global stiffness matrix

        # Global equilibrium loop
        while (fErr > tol and NRit < NRitMax) or NRit < 2:
            # Linear solver
            duvw, drct = linSolve(mesh['bc'], Kt, oobf, NRit, fd)

            # Update displacements
            uvw = uvw + duvw

            # Update reaction forces
            frct = frct + drct

            # Global stiffness & internal force
            fint, Kt, mpData = detMPs(uvw, mpData)

            # Out-of-balance force
            oobf = fext - fint + frct

            # Normalised oobf error
            fErr = np.linalg.norm(oobf, 2) / (np.linalg.norm(fext + frct, 2) +
                                              np.finfo(float).eps)

            # Increment the NR counter
            NRit = NRit + 1

            # Text output to screen (NR error)
            print(f'  iteration {NRit:2d} NR error {fErr:8.3e}')

        # Update material points
        mpData = updateMPs(uvw.reshape(-1, 1), mpData)

        # Plotting and post processing
        postPro(lstp, mesh, mpData, problem_type)


if __name__ == "__main__":
    import sys
    if len(sys.argv) > 1:
        main(sys.argv[1])
    else:
        main()