# 9.5.8 Subroutine MINDPB
This subroutine simply reads some additional information required for controlling the convergence check and inserting additional constraints for the Heterosis element.
```asm
SUBROUTINE MINDPB (IFDIS, IFFIX, IFRES, LNODS, MELEM, MTOTV, MIND 1
. NCDIS, NCRES, NELEM, NTYPE) MIND 2
C*****************************************************************************************MIND 3
C MIND 4
C*** READS ADDITIONAL DATA FOR MINDLIN PLATE ANALYSIS MIND 5
C MIND 6
C*****************************************************************************************MIND 7
DIMENSION DERIV(2,9), IFFIX(MTOTV), MIND 8
. LNODS(MELEM,9), NCDIS(4), NCRES(4), SHAPE(9) MIND 9
C MIND 10
C*** READ DATA CONTROLLING CONVERGENCE CHECK MIND 11
C MIND 12
10 READ(5,900) IFDIS, (NCDIS(I), I=1,4) MIND 13
. ,IFRES, (NCRES(I), I=1,4) MIND 14
900 FORMAT(5I1) MIND 15
WRITE(6,901) IFDIS, (NCDIS(I), I=1,4) MIND 16
. ,IFRES, (NCRES(I), I=1,4) MIND 17
901 FORMAT(/,23H CONVERGENCE PARAMETERS,/, MIND 18
. 8H IFDIS =, I2,5X, 8H NCDIS =, 4I1,/, MIND 19
. 8H IFRES =, I2,5X, 8H NCRES =, 4I1, //) MIND 20
C*** INSERT ADDITIONAL CONSTRAINT FOR HETEROSIS ELEMENT MIND 21
IF(NTYPE.NE.5) GO TO 30 MIND 22
DO 20 IELEM=1, NELEM MIND 23
LNODE=LNODS(IELEM,9) MIND 24
NLOCA=LNODE*3-2 MIND 25
20 IFFIX(NLOCA)=-1 MIND 26
30 CONTINUE MIND 27
RETURN MIND 28
END MIND 29
```
# 9.5.9 Subroutine NODEXY
This subroutine evaluates midside nodes for straight sided 8 and 9-node quadrilateral elements. In the original subroutine described in Section 6.4.1 this routine also evaluated the coordinates of the central node. Here, as we are choosing a hierarchical formulation, the values at the central node and the departures from the interpolated Serendipity values are always taken as zero.
Thus the revised subroutine NODEXY is almost identical to its namesake given earlier in Section 6.4.1 and is listed below.
```c
SUBROUTINE NODEXY (COORD, LNODS, MELEM, MPOIN, NDIME, NELEM, NNODE) NODE 1
C***** NODE 2
C***** NODE 3
C***** NODE 4
C*** INTERPOLATES MIDSIDE NODE COORDINATES FOR 8-NODED ELEMENTS NODE 5
C*** INTERPOLATES CENTRAL AND MIDSIDE NODE COORDINATES FOR NODE 6
C*** 9-NODE ELEMENTS PROVIDED THAT THE SIDES ARE STRAIGHT NODE 7
C***** NODE 8
C***** NODE 9
```
```fortran
DIMENSION COORD(MPOIN,2),LNODS(MELEM,9) NODE 10
IF(NNODE.EQ.4) GO TO 60 NODE 11
C NODE 12
C*** LOOP OVER EACH ELEMENT NODE 13
C NODE 14
DO 30 IELEM=1,NELEM NODE 15
C NODE 16
C*** LOOP OVER EACH ELEMENT EDGE NODE 17
C NODE 18
NNOD1=NNODE NODE 19
IF(NNODE.EQ.8).NNOD1=9 NODE 20
DO 20 INODE=1,NNOD1,2 NODE 21
IF(INODE.EQ.9.AND.NNODE.EQ.8) GO TO 30 NODE 22
IF(INODE.EQ.9) GO TO 50 NODE 23
C NODE 24
C*** COMPUTE THE NODE NUMBER OF THE FIRST NODE NODE 25
C NODE 26
NODST=LNODS(IELEM,INODE) NODE 27
IGASH=INODE+2 NODE 28
IF(IGASH.GE.NNODE) IGASH=1 NODE 29
C NODE 30
C*** COMPUTE THE NODE NUMBER OF THE LAST NODE NODE 31
C NODE 32
NODFN=LNODS(IELEM,IGASH) NODE 33
MIDPT=INODE+1 NODE 34
C NODE 35
C*** COMPUTE THE NODE NUMBER OF THE INTERMEDIATE NODE NODE 36
C NODE 37
NODMD=LNODS(IELEM,MIDPT) NODE 38
TOTAL=ABS(COORD(NODMD,1))+ABS(COORD(NODMD,2)) NODE 39
C NODE 40
C*** IF THE COORDINATES OF THE INTERMEDIATE NODE ARE BOTH ZERO NODE 41
C INTERPOLATE BY A STRAIGHT LINE NODE 42
C NODE 43
IF(TOTAL.GT.0.0) GO TO 20 NODE 44
KOUNT=1 NODE 45
10 COORD(NODMD,KOUNT)=(COORD(NODST,KOUNT)+COORD(NODFN,KOUNT))/2.0 NODE 46
KOUNT=KOUNT+1 NODE 47
IF(KOUNT.EQ.2) GO TO 10 NODE 48
20 CONTINUE NODE 49
50 LNODE=LNODS(IELEM,INODE) NODE 50
30 CONTINUE NODE 51
60 CONTINUE NODE 52
RETURN NODE 53
END NODE 54
```
# 9.5.10 Subroutine OUTMP
This subroutine outputs nodal displacements and reactions and also the Gauss point stress resultants.
```csv
SUBROUTINE OUTMP (EPSTN,IITER,MTOTG,MTOTV,MVFIX,NELEM, OUTP 1
. NGAUS,NOFIX,NOUTP,NPOIN,NVFIX,STRSG, OUTP 2
. TDISP,TREAC) OUTP 3
C**********OUTP 4
C OUTP 5
C*** OUTPUT DISPLACEMENTS,REACTIONS AND GAUSS POINT STRESS OUTP 6
C*** RESULTANTS FOR EP MINDLIN PLATE ANALYSIS OUTP 7
C OUTP 8
C**********OUTP 9
```
```csv
DIMENSION EPSTN(MTOTG), GPCOD(2,9), NOFIX(MVFIX), NOUTP(2), OUTP 10
. STRSG(5, MTOTG), TDISP(MTOTV), TREAC(MVFIX, 3) OUTP 11
KOUTP=NOUTP(1) OUTP 12
IF(IITER.GT.1) KOUTP=NOUTP(2) OUTP 13
C
C*** OUTPUT DISPLACEMENTS OUTP 15
C
IF(KOUTP.LT.1) GO TO 10 OUTP 17
WRITE(6,900) OUTP 18
900 FORMAT(1H0,5X,13HDISPLACEMENTS) OUTP 19
WRITE(6,950) OUTP 20
950 FORMAT(1H0,6X,4HNODE,6X,5HDISP.,8X,7HXZ-ROT.,7X,7HYZ-ROT.) OUTP 21
DO 20 IPOIN=1,NPOIN OUTP 22
NGASH=IPOIN*3 OUTP 23
NGISH=NGASH-3+1 OUTP 24
20 WRITE(6,910) IPOIN,(TDISP(IGASH),IGASH=NGISH,NGASH) OUTP 25
910 FORMAT(I10,3E14.6) OUTP 26
10 CONTINUE OUTP 27
C
C*** OUTPUT REACTIONS OUTP 28
C
IF(KOUTP.LT.2) GO TO 30 OUTP 30
WRITE(6,920) OUTP 31
920 FORMAT(1H0,5X,9HREACTIONS) OUTP 32
WRITE(6,960) OUTP 33
960 FORMAT(1H0,6X,4HNODE,6X,5HFORCE,3X,9HXZ-MOMENT,5X,9HYZ-MOMENT) OUTP 34
DO 40 IVFIX=1,NVFIX OUTP 35
40 WRITE(6,910) NOFIX(IVFIX),(TREAC(IVFIX,IDOFN),IDOFN=1,3) OUTP 36
30 CONTINUE OUTP 37
C
C*** OUTPUT STRESSES OUTP 38
C
IF(KOUTP.LT.3) GO TO 50 OUTP 39
REWIND 3 OUTP 40
WRITE(6,970) OUTP 41
970 FORMAT(1H0,5X,8HSTRESSES) OUTP 42
WRITE(6,980) OUTP 43
980 FORMAT(1H0,4HG.P.,2X,8HX-COORD.,2X,8HY-COORD.,3X,8HX-MOMENT,4X,OUTP 44
.8HY-MOMENT,3X,9HXY-MOMENT,3X,OUTP 45
.13HEFF.PL.STRAIN) OUTP 46
KGAUS=0 OUTP 47
DO 60 IELEM=1,NELEM OUTP 48
READ(3)GPCOD OUTP 49
KELGS=0 OUTP 50
WRITE(6,930)IELEM OUTP 51
930 FORMAT(1H0,5X,13HELEMENT NO. =,I5) OUTP 52
DO 60 IGAUS=1,NGAUS OUTP 53
DO 60 JGAUS=1,NGAUS OUTP 54
KGAUS=KGAUS+1 OUTP 55
KELGS=KELGS+1 OUTP 56
WRITE(6,940)KELGS,(GPCOD(IDIME,KELGS),IDIME=1,2), OUTP 57
.(STRSG(ISTRE,KGAUS),ISTRE=1,3),EPSTN(KGAUS) OUTP 58
940 FORMAT(I5,2F10.4,6E12.5) OUTP 59
60 CONTINUE OUTP 60
50 CONTINUE OUTP 61
RETURN OUTP 62
END OUTP 63
OUTP 64
OUTP 65
OUTP 66
```
# 9.5.11 Subroutine RESMP
This subroutine evaluates the residual nodal forces. The structure of this routine is similar to that given in Chapter 7 for the other two dimensional elasto-plastic applications and it is illustrated in Fig. 9.3.
```fortran
SUBROUTINE RESMP (ASDIS,COORD,EFFST,ELOAD,EPSTN,LNODS, RESP 1
MATNO,MELEM,MMATS,MPOIN,MTOTG,MTOTV, RESP 2
NCRIT,NELEM,NEVAB,NGAUS,NNODE,PROPS, RESP 3
STRSG) RESP 4
C**************************RESP 5
C
C*** EVALUATES EQUIVALENT NODAL FORCES FOR THE STRESS RESULTANTS RESP 7
C*** IN MINDLIN PLATES DURING EP ANALYSIS RESP 8
C
C**************************RESP 9
DIMENSION ASDIS(MTOTV),AVECT(5),CARTD(2,9), RESP 11
COORD(MPOIN,2),DERIV(2,9),DESIG(5),DEVIA(4), RESP 12
DVECT(5), RESP 13
EFFST(MTOTG),ELCOD(2,9), RESP 14
ELDIS(3,9),ELOAD(MELEM,27),EPSTN(MTOTG),GPCOD(2,9), RESP 15
LNODS(MELEM,9),MATNO(MELEM),POSGP(4), RESP 16
PROPS(MMATS,8),SGTOT(5),SHAPE(9),SIGMA(5), RESP 17
STRES(5),STRSG(5,MTOTG),WEIGP(4), RESP 18
DFLEX(3,3),DSHER(2,2),BFLEI(3,3),BSHEI(2,3), RESP 19
DUMMY(3,3),FORCE(3),DGRAD(6) RESP 20
NTIME=1 RESP 21
DO 10 IELEM=1,NELEM RESP 22
DO 10 IEVAB=1,NEVAB RESP 23
10 ELOAD(IELEM,IEVAB)=0.0 RESP 24
KGAUS=0 RESP 25
LGAUS=0 RESP 26
DO 20 IELEM=1,NELEM RESP 27
LPROP=MATNO(IELEM) RESP 28
C
C*** COMPUTE COORDINATE AND INCREMENTAL DISPLACEMENTS OF THE RESP 30
ELEMENT NODAL POINTS RESP 31
C
DO 190 INODE =1,NNODE RESP 32
LNODE=IABS(LNODS(IELEM,INODE)) RESP 33
NPOSN=(LNODE-1)*3 RESP 34
DO 30 IDOFN=1,3 RESP 35
NPOSN=NPOSN+1 RESP 36
30 ELDIS(IDOFN,INODE)=ASDIS(NPOSN) RESP 37
DO 180 IDIME=1,2 RESP 38
180 ELCOD(IDIME,INODE)=COORD(LNODE,IDIME) RESP 39
190 CONTINUE RESP 40
KGASP=0 RESP 41
CALL MODPB (DFLEX,DUMMY,DSHER,LPROP,MMATS,PROPS, RESP 42
0, 1, 1) RESP 43
CALL GAUSSQ (NGAUS,POSGP,WEIGP) RESP 44
DO 40 IGAUS=1,NGAUS RESP 45
DO 40 JGAUS=1,NGAUS RESP 46
BRING=1.0 RESP 47
KGAUS=KGAUS+1 RESP 48
EXISP=POSGP(IGAUS) RESP 49
ETASP=POSGP(JGAUS) RESP 50
CALL SFR2 (DERIV,ETASP,EXISP,NNODE,SHAPE) RESP 51
KGASP=KGASP+1 RESP 52
CALL JACOB2 (CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM, RESP 53
KGASP,NNODE,SHAPE) RESP 54
C
```
```csv
DAREA=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)
CALL GRADMP (CARTD,DGRAD,ELDIS, 3,NNODE)
CALL STRMP (CARTD,DFLEX,DGRAD,DSHER,ELDIS,NNODE, SHAPE,STRES, 1, 0)
PREYS=PROPS(LPROP,6)+EPSTN(KGAUS)*PROPS(LPROP,7)
DO 150 ISTRE=1,3
DESIG(ISTRE)=STRES(ISTRE)
150 SIGMA(ISTRE)=STRSG(ISTRE,KGAUS)+STRES(ISTRE)
CALL INVMP (DEVIA,NCRIT,SINT3,STEFF,SIGMA,THETA, VARJ2,YIELD)
ESPRE=EFFST(KGAUS)-PREYS
IF(ESPRE.GE.0.0) GO TO 50
ESCUR=YIELD-PREYS
IF(ESCUR.LE.0.0) GO TO 60
RFACT=ESCUR/(YIELD-EFFST(KGAUS))
GO TO 70
50 ESCUR=YIELD-EFFST(KGAUS)
IF(ESCUR.LE.0.0) GO TO 60
RFACT=1.0
70 MSTEP=ESCUR*8.0/PROPS(LPROP,6)+1.0
ASTEP=MSTEP
REDUC=1.0-RFACT
DO 80 ISTRE=1,3
SGTOT(ISTRE)=STRSG(ISTRE,KGAUS)+REDUC*STRES(ISTRE)
80 STRES(ISTRE)=RFACT*STRES(ISTRE)/ASTEP
DO 90 ISTEP=1,MSTEP
CALL INVMP (DEVIA,NCRIT,SINT3,STEFF,SGTOT,THETA, VARJ2,YIELD)
HARDS=PROPS(LPROP,7)
CALL FLOWMP (ABETA,AVECT,DEVIA,DFLEX,DVECT,HARDS, NCRIT,SINT3,STEFF,THETA,VARJ2)
AGASH=0.0
DO 100 ISTRE=1,3
100 AGASH=AGASH+AVECT(ISTRE)*STRES(ISTRE)
DLAMD=AGASH*ABETA
IF(DLAMD.LT.0.0) DLAMD=0.0
BGASH=0.0
DO 110 ISTRE=1,3
BGASH=BGASH+AVECT(ISTRE)*SGTOT(ISTRE)
110 SGTOT(ISTRE)=SGTOT(ISTRE)+STRES(ISTRE)-DLAMD*DVECT(ISTRE)
90 EPSTN(KGAUS)=EPSTN(KGAUS)+DLAMD*BGASH/YIELD
DO 120 ISTRE=1,3
120 DESIG(ISTRE)=SGTOT(ISTRE)-STRSG(ISTRE,KGAUS)
CALL INVMP (DEVIA,NCRIT,SINT3,STEFF,SGTOT,THETA, VARJ2,YIELD)
CURYS=PROPS(LPROP,6)+EPSTN(KGAUS)*PROPS(LPROP,7)
IF(YIELD.GT.CURYS) BRING=CURYS/YIELD
60 DO 130 ISTRE=1,3
SGTOT(ISTRE)=BRING*(STRSG(ISTRE,KGAUS)+DESIG(ISTRE))
130 STRSG(ISTRE,KGAUS)=SGTOT(ISTRE)
EFFST(KGAUS)=BRING*YIELD
*** CALCULATE THE EQUIVALENT NODAL FORCES AND ASSOCIATE WITH THE ELEMENT NODES
DO 140 INODE=1,NNODE
*** ZERO FORCE VECTOR
CALL VZERO (3,FORCE)
CALL BMATPB (BFLEI,DUMMY,BSHEI,CARTD,INODE,SHAPE, 0, 1, 0)
FORCE(2)=(BFLEI(1,2)*SGTOT(1)+BFLEI(3,2)*SGTOT(3))*DAREA
+FORCE(2)
FORCE(3)=(BFLEI(2,3)*SGTOT(2)+BFLEI(3,3)*SGTOT(3))*DAREA
+FORCE(3)
IPOSN=(INODE-1)*3+1
RESP 56
RESP 57
RESP 58
RESP 59
RESP 60
RESP 61
RESP 62
RESP 63
RESP 64
RESP 65
RESP 66
RESP 67
RESP 68
RESP 69
RESP 70
RESP 71
RESP 72
RESP 73
RESP 74
RESP 75
RESP 76
RESP 77
RESP 78
RESP 79
RESP 80
RESP 81
RESP 82
RESP 83
RESP 84
RESP 85
RESP 86
RESP 87
RESP 88
RESP 89
RESP 90
RESP 91
RESP 92
RESP 93
RESP 94
RESP 95
RESP 96
RESP 97
RESP 98
RESP 99
RESP 100
RESP 101
RESP 102
RESP 103
RESP 104
RESP 105
RESP 106
RESP 107
RESP 108
RESP 109
RESP 110
RESP 111
RESP 112
RESP 113
RESP 114
RESP 115
RESP 116
RESP 117
RESP 118
RESP 119
```
```csv
DO 135 IDOFN=2,3
IPOSN=IPOSN+1
135 ELOAD(IELEM,IPOSN)=ELOAD(IELEM,IPOSN)+FORCE(IDOFN)
140 CONTINUE
40 CONTINUE
C
C*** CALCULATE FORCES ASSOCIATED WITH SHEAR DEFORMATION
C
NGAUM=NGAUS-1
CALL GAUSSQ (NGAUM,POSGP,WEIGP)
C
C*** ENTER LOOPS FOR AREA NUMERICAL INTEGRATION
C
KGASP=0
DO 300 IGAUS=1,NGAUM
DO 300 JGAUS=1,NGAUM
LGAUS=LGAUS+1
EXISP=POSGP(IGAUS)
ETASP=POSGP(JGAUS)
CALL SFR2 (DERIV,ETASP,EXISP,NNODE,SHAPE)
KGASP=KGASP+1
CALL JACOB2 (CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM, KGASP,NNODE,SHAPE)
DAREA=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)
CALL GRADMP (CARTD,DGRAD,ELDIS, 3,NNODE)
CALL STRMP (CARTD,DFLEX,DGRAD,DSHER,ELDIS,NNODE, SHAPE,STRES, 0, 1)
DO 310 ISTRE=4,5
SGTOT(ISTRE)=STRSG(ISTRE,LGAUS)+STRES(ISTRE)
310 STRSG(ISTRE,LGAUS)=SGTOT(ISTRE)
C
C*** CALCULATE THE EQUIVALENT NODAL FORCES
C
DO 320 INODE=1,NNODE
C*** ZERO FORCE VECTOR
CALL VZERO(3,FORCE)
CALL BMATPB (BFLEI,DUMMY,BSHEI,CARTD,INODE,SHAPE, 0, 0, 1)
FORCE(1)=(BSHEI(1,1)*SGTOT(4)+BSHEI(2,1)*SGTOT(5))*DAREA
+FORCE(1)
FORCE(2)=(BSHEI(1,2)*SGTOT(4))*DAREA+FORCE(2)
FORCE(3)=(BSHEI(2,3)*SGTOT(5))*DAREA+FORCE(3)
IPOSN=(INODE-1)*3
DO 315 IDOFN=1,3
IPOSN=IPOSN+1
315 ELOAD(IELEM,IPOSN)=ELOAD(IELEM,IPOSN)+FORCE(IDOFN)
320 CONTINUE
300 CONTINUE
20 CONTINUE
RETURN
END
RESP 120
RESP 121
RESP 122
RESP 123
RESP 124
RESP 125
RESP 126
RESP 127
RESP 128
RESP 129
RESP 130
RESP 131
RESP 132
RESP 133
RESP 134
RESP 135
RESP 136
RESP 137
RESP 138
RESP 139
RESP 140
RESP 141
RESP 142
RESP 143
RESP 144
RESP 145
RESP 146
RESP 147
RESP 148
RESP 149
RESP 150
RESP 151
RESP 152
RESP 153
RESP 154
RESP 155
RESP 156
RESP 157
RESP 158
RESP 159
RESP 160
RESP 161
RESP 162
RESP 163
RESP 164
RESP 165
RESP 166
RESP 167
RESP 168
RESP 169
RESP 170
```
# 9.5.12 Subroutine SFR2
This subroutine evaluates the shape functions and their derivatives for 4, 8 and 9-node quadrilateral isoparametric elements. The 9-node element is treated as a hierarchical element as described in Section 9.3.2. This enables the Heterosis element to be easily incorporated.
flowchart
```mermaid
graph TD
A["START"] --> B["Set to zero ELOAD ( . . . )"]
B --> C["Extract local element material property set number, displacements and coordinates"]
C --> D["Call MODPB to evaluate Df, Ds"]
D --> E["Call GAUSSQ to evaluate n-point Gauss-Legendre sampling positions and weights"]
E --> F["Call SFR2, JACOB2, GRADMP and STRMP to evaluate elastic stress increment dσf"]
F --> G["Calculate the effective stress necessary for yielding to occur"]
G --> H["Calculate the total bending moments at the current Gauss points"]
H --> I["If the current bending moments are outside of the yield surface bring them back to the yield surface taking into account unloading if it has taken place"]
I --> J["A"]
K["ELEMENT LOOP"] --> L["GAUSS LOOPS"]
L --> M["Calculate the effective stress necessary for yielding to occur"]
M --> N["Calculate the total bending moments at the current Gauss points"]
N --> O["If the current bending moments are outside of the yield surface bring them back to the yield surface taking into account unloading if it has taken place"]
```
Fig. 9.3 Overall structure of subroutine RESMP.
flowchart
```mermaid
graph TD
A["A"] --> B["Evaluate [Bf"]T σ̂f × Gauss weights × det J and add into ELOAD ( , ). Use routines VZERO and BMATPB]
B --> C["Call GAUSSQ to evaluate (n - 1) point Gauss-Legendre sampling positions and weights"]
C --> D["Call SFR2, JACOB2, GRADMP and STRMP to evaluate elastic stress increment dσ̂s"]
D --> E["Evaluate [Bs"]T σ̂s × Gauss weights × det J and add into ELOAD ( , ). Use routines VZERO and BMATPB]
E --> F["RETURN"]
G["GAUSS LOOPS"] --> D
G --> E
```
Fig. 9.3 Overall structure of subroutine RESMP (continued).
Subroutine SFR2 is identical to its namesake given earlier in Section 6.4.3 except that SFR2 72–118 are replaced by SFRH 67–73.
IF(NNODE.EQ.8) RETURN
SFR2
67
C*** BUBBLE FUNCTION FOR HIERARCHICAL AND HETEROSIS ELEMENTS
SFRH
68
SHAPE(9)=(1.0-SS)*(1.0-TT)
SFRH
69
DERIV(1,9)=-S2*(1.0-TT)
SFRH
70
DERIV(2,9)=-T2*(1.0-SS)
SFRH
71
RETURN
SFRH
72
END
SFRH
73
# 9.5.13 Subroutine STIFMP
This routine evaluates the stiffness matrix for the nonlayered elasto-plastic Mindlin plate elements. The overall structure is shown in Fig. 9.4.
flowchart
```mermaid
graph TD
A["START"] --> B["Rewind tapes 1 and 3"]
B --> C["Extract local element material property set number and coordinates"]
C --> D["Initialise array used to store element stiffness matrices"]
D --> E["Call GAUSSQ to evaluate n-point Gauss-Legendre sampling positions and weights"]
E --> F["Call SFR2 and JACOB2 to evaluate N_i^(e) / ∂N_i^(e)/∂x, ∂N_i^(e)/∂y and det J^(e)"]
F --> G["Call MODPB to evaluate D_f"]
G --> H{Is this the first load increment?}
H -->|Yes| I["GAUSS LOOPS"]
H -->|No| J{Was this Gauss point plastic in the last load increment?}
J -->|No| I
J -->|Yes| K["A"]
```
flowchart
```mermaid
graph TD
A["A"] --> B["Call INVMP and FLOWMP to evaluate a' and d_D and hence calculate D_ep"]
B --> C["Call BMATPB and SUBMP to add [B_f^(e)"]^T D_f B_f^(e) det J^(e) × Gauss weights into K_ij^(e)]
C --> D["Call GAUSSQ to evaluate (n-1)-point Gauss-Legendre sampling positions and weights"]
D --> E["Call SFR2 and JACOB2 to evaluate N_i^(e), ∂N_i^(e)/∂x, ∂N_i^(e)/∂y and det J^(e)"]
E --> F["Call MODPB to evaluate D_s"]
F --> G["Call BMATPB and SUBMP to add [B_δi^(e)"]^T D_δ B_δj^(e) det J × Gauss weights into K_ij^(e)]
G --> H["Store stiffness matrix K^(e) and Gauss point coordinates on files 1 and 3 respectively"]
H --> I["RETURN"]
I --> J["GAUSS LOOPS"]
```
Fig. 9.4 Overall structure of subroutine STIFMP (continued).