Skip to content

Commit 097b4c2

Browse files
zchlrnrzchlrnr
authored
Merge pull request #179 from victorkemp/mass_matrices
Mass matrices
2 parents b8b67d0 + 380d84f commit 097b4c2

6 files changed

Lines changed: 198 additions & 72 deletions

File tree

‎Source/EMG/EMG4/MITC4.f90‎

Lines changed: 4 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
3939
USE SCONTR, ONLY : BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_GAUSS, MAX_STRESS_POINTS, NTSUB, NSUB
4040
USE NONLINEAR_PARAMS, ONLY : LOAD_ISTEP
4141
USE MODEL_STUF, ONLY : NUM_EMG_FATAL_ERRS, PCOMP_PROPS, ELGP, ES, KE, EM, EB, ET, BE1, BE2, BE3, PHI_SQ, &
42-
FCONV, EPROP, PTE, ALPVEC, TREF, DT, PPE, PRESS, MASS_PER_UNIT_AREA, ME, &
42+
FCONV, EPROP, PTE, ALPVEC, TREF, DT, PPE, PRESS, MASS_PER_UNIT_AREA, &
4343
NUM_PLIES, PCOMP_LAM, PLY_NUM, TPLY, STRESS, KED
4444
USE CONSTANTS_1, ONLY : ZERO, ONE, TWO, FOUR
4545

@@ -60,6 +60,7 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
6060
USE CROSS_Interface
6161
USE MITC_SHAPE_FUNCTIONS_Interface
6262
USE MITC_COVARIANT_BASIS_Interface
63+
USE EXPAND_MASS_DOFS_Interface
6364

6465
IMPLICIT NONE
6566

@@ -112,7 +113,7 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
112113
REAL(DOUBLE) :: BMI(6,6*ELGP) ! Strain-displ matrix for membrane for one Gauss point
113114
REAL(DOUBLE) :: BBI(6,6*ELGP) ! Strain-displ matrix for bending for one Gauss point
114115
REAL(DOUBLE) :: BSI(6,6*ELGP) ! Strain-displ matrix for shear for one Gauss point
115-
REAL(DOUBLE) :: M_1DOF(ELGP,ELGP)
116+
REAL(DOUBLE) :: M_1DOF(ELGP,ELGP) ! Consistent mass matrix with 1 DOF per node.
116117
REAL(DOUBLE) :: DENSITY
117118
REAL(DOUBLE) :: FORCEx(IORD_STRESS_Q4*IORD_STRESS_Q4) ! Engineering force in the elem x direction at Gauss points
118119
REAL(DOUBLE) :: FORCEy(IORD_STRESS_Q4*IORD_STRESS_Q4) ! Engineering force in the elem x direction at Gauss points
@@ -132,7 +133,6 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
132133
REAL(DOUBLE) :: DUM14(3,3)
133134
REAL(DOUBLE) :: DUM33(3,3)
134135
REAL(DOUBLE) :: JAC2x2(2,2)
135-
REAL(DOUBLE) :: ROW_SUM
136136

137137
! **********************************************************************************************************************************
138138

@@ -209,7 +209,6 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
209209
! Consistent mass matrix
210210
! ME = ∫ N' ρ N det(J) dv
211211

212-
ME(:,:) = ZERO
213212
M_1DOF(:,:) = ZERO
214213

215214
DENSITY = MASS_PER_UNIT_AREA / EPROP(1)
@@ -240,26 +239,7 @@ SUBROUTINE MITC4 ( OPT, INT_ELEM_ID )
240239
ENDDO
241240
ENDDO
242241

243-
! Row sum to convert to lumped mass Matrix
244-
DO I=1,ELGP
245-
ROW_SUM = ZERO
246-
DO J=1,ELGP
247-
ROW_SUM = ROW_SUM + M_1DOF(I,J)
248-
M_1DOF(I,J) = ZERO
249-
ENDDO
250-
M_1DOF(I,I) = ROW_SUM
251-
ENDDO
252-
253-
! Copy to all 3 translational DOFs of each node.
254-
DO I=1,ELGP
255-
DO J=1,ELGP
256-
KI = (I-1) * 6
257-
KJ = (J-1) * 6
258-
DO L=1,3
259-
ME(KI + L, KJ + L) = M_1DOF(I,J)
260-
ENDDO
261-
ENDDO
262-
ENDDO
242+
CALL EXPAND_MASS_DOFS( M_1DOF )
263243

264244

265245
ENDIF

‎Source/EMG/EMG5/HEXA.f90‎

Lines changed: 25 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -41,16 +41,17 @@ SUBROUTINE HEXA ( OPT, INT_ELEM_ID,IORD, RED_INT_SHEAR, WRITE_WARN )
4141
USE IOUNT1, ONLY : WRT_ERR, WRT_LOG, ERR, F04, F06
4242
USE SCONTR, ONLY : BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_GAUSS, MELDOF, MPLOAD4_3D_DATA, NPLOAD4_3D, NSUB, NTSUB
4343
USE TIMDAT, ONLY : TSEC
44-
USE CONSTANTS_1, ONLY : QUARTER, HALF, ZERO, ONE, EIGHT
44+
USE CONSTANTS_1, ONLY : QUARTER, HALF, ZERO, ONE
4545
USE DEBUG_PARAMETERS, ONLY : DEBUG
4646
USE SUBR_BEGEND_LEVELS, ONLY : HEXA_BEGEND
4747
USE PARAMS, ONLY : EPSIL
4848
USE NONLINEAR_PARAMS, ONLY : LOAD_ISTEP
4949
USE MODEL_STUF, ONLY : AGRID, ALPVEC, BE1, BE2, DT, EID, ELGP, NUM_EMG_FATAL_ERRS, ES, KE, KED, ME, &
5050
NUM_EMG_FATAL_ERRS, PLOAD4_3D_DATA, PPE, PRESS, PTE, RHO, SE1, SE2, STE1, STRESS, TREF, &
5151
TYPE, XEL
52-
5352
USE HEXA_USE_IFs
53+
USE EXPAND_MASS_DOFS_Interface
54+
5455

5556
IMPLICIT NONE
5657

@@ -126,7 +127,6 @@ SUBROUTINE HEXA ( OPT, INT_ELEM_ID,IORD, RED_INT_SHEAR, WRITE_WARN )
126127
REAL(DOUBLE) :: JAC(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian matrix.
127128
REAL(DOUBLE) :: JACI(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian inverse.
128129
REAL(DOUBLE) :: KWW(3,3) ! Portion of differential stiffness matrix
129-
REAL(DOUBLE) :: M0 ! An intermediate variable used in calc elem mass, ME
130130
REAL(DOUBLE) :: PSH(ELGP) ! Output from subr SHP3DH. Shape fcn at Gauss pts SSI, SSJ
131131
REAL(DOUBLE) :: PSIGN !
132132
REAL(DOUBLE) :: SIGxx ! Normal stress in the elem x direction
@@ -145,7 +145,8 @@ SUBROUTINE HEXA ( OPT, INT_ELEM_ID,IORD, RED_INT_SHEAR, WRITE_WARN )
145145
REAL(DOUBLE) :: TGAUSS(1,NTSUB) ! Temp at a Gauss point for a theral subcase
146146
REAL(DOUBLE) :: VOLUME ! 3D element volume
147147
REAL(DOUBLE) :: SSI,SSJ,SSK ! Isoparametric coordinates of a point.
148-
148+
REAL(DOUBLE) :: M_1DOF(ELGP,ELGP) ! Consistent mass matrix with 1 DOF per node.
149+
149150
! **********************************************************************************************************************************
150151
IF (WRT_LOG >= SUBR_BEGEND) THEN
151152
CALL OURTIM
@@ -224,39 +225,32 @@ SUBROUTINE HEXA ( OPT, INT_ELEM_ID,IORD, RED_INT_SHEAR, WRITE_WARN )
224225

225226
IF (OPT(1) == 'Y') THEN
226227

227-
M0 = (RHO(1))*VOLUME/EIGHT
228-
229-
ME( 1 ,1) = M0
230-
ME( 2 ,2) = M0
231-
ME( 3 ,3) = M0
232-
233-
ME( 7 ,7) = M0
234-
ME( 8 ,8) = M0
235-
ME( 9 ,9) = M0
228+
! Consistent mass matrix
229+
! ME = ∫ N' ρ N det(J) dv
236230

237-
ME(13,13) = M0
238-
ME(14,14) = M0
239-
ME(15,15) = M0
231+
M_1DOF(:,:) = ZERO
232+
GAUSS_PT = 0
233+
CALL ORDER_GAUSS ( IORD, SSS, HHH )
234+
DO K=1,IORD
235+
DO J=1,IORD
236+
DO I=1,IORD
237+
GAUSS_PT = GAUSS_PT + 1
240238

241-
ME(19,19) = M0
242-
ME(20,20) = M0
243-
ME(21,21) = M0
239+
CALL SHP3DH ( I, J, K, ELGP, SUBR_NAME, IORD_MSG, IORD, SSS(I), SSS(J), SSS(K), 'N', PSH, DPSHG )
240+
INTFAC = DETJ(GAUSS_PT)*HHH(I)*HHH(J)*HHH(K)
244241

245-
ME(25,25) = M0
246-
ME(26,26) = M0
247-
ME(27,27) = M0
242+
DO L=1,ELGP
243+
DO M=1,ELGP
244+
M_1DOF(L,M) = M_1DOF(L,M) + PSH(L) * PSH(M) * RHO(1) * INTFAC
245+
ENDDO
246+
ENDDO
248247

249-
ME(31,31) = M0
250-
ME(32,32) = M0
251-
ME(33,33) = M0
248+
ENDDO
249+
ENDDO
250+
ENDDO
252251

253-
ME(37,37) = M0
254-
ME(38,38) = M0
255-
ME(39,39) = M0
256252

257-
ME(43,43) = M0
258-
ME(44,44) = M0
259-
ME(45,45) = M0
253+
CALL EXPAND_MASS_DOFS( M_1DOF )
260254

261255
ENDIF
262256

‎Source/EMG/EMG5/PENTA.f90‎

Lines changed: 23 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ SUBROUTINE PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN
4040
USE IOUNT1, ONLY : WRT_ERR, WRT_LOG, ERR, F04, F06
4141
USE SCONTR, ONLY : BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_TRIA, MAX_ORDER_GAUSS, NTSUB
4242
USE TIMDAT, ONLY : TSEC
43-
USE CONSTANTS_1, ONLY : HALF, THIRD, ZERO, SIX
43+
USE CONSTANTS_1, ONLY : HALF, THIRD, ZERO
4444
USE DEBUG_PARAMETERS, ONLY : DEBUG
4545
USE SUBR_BEGEND_LEVELS, ONLY : PENTA_BEGEND
4646
USE PARAMS, ONLY : EPSIL
@@ -49,7 +49,8 @@ SUBROUTINE PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN
4949
SE1, SE2, STE1, STRESS, TREF, TYPE
5050

5151
USE PENTA_USE_IFs
52-
52+
USE EXPAND_MASS_DOFS_Interface
53+
5354
IMPLICIT NONE
5455

5556
CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'PENTA'
@@ -104,7 +105,6 @@ SUBROUTINE PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN
104105
REAL(DOUBLE) :: JAC(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian matrix.
105106
REAL(DOUBLE) :: JACI(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian inverse.
106107
REAL(DOUBLE) :: KWW(3,3) ! Portion of differential stiffness matrix
107-
REAL(DOUBLE) :: M0 ! An intermediate variable used in calc elem mass, ME
108108
REAL(DOUBLE) :: PSH(ELGP) ! Output from subr SHP3DP. Shape fcn at Gauss pts SSI, SSJ
109109
REAL(DOUBLE) :: SIGxx ! Normal stress in the elem x direction
110110
REAL(DOUBLE) :: SIGyy ! Normal stress in the elem y direction
@@ -126,7 +126,7 @@ SUBROUTINE PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN
126126
REAL(DOUBLE) :: TGAUSS(1,NTSUB) ! Temp at a Gauss point for a theral subcase
127127
REAL(DOUBLE) :: VOLUME ! 3D element volume
128128
REAL(DOUBLE) :: SSI,SSJ,SSK ! Isoparametric coordinates of a point.
129-
129+
REAL(DOUBLE) :: M_1DOF(ELGP,ELGP) ! Consistent mass matrix with 1 DOF per node.
130130

131131
! **********************************************************************************************************************************
132132
IF (WRT_LOG >= SUBR_BEGEND) THEN
@@ -205,15 +205,29 @@ SUBROUTINE PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN
205205

206206
IF (OPT(1) == 'Y') THEN
207207

208-
M0 = (RHO(1))*VOLUME/ELGP
208+
! Consistent mass matrix
209+
! ME = ∫ N' ρ N det(J) dv
210+
211+
M_1DOF(:,:) = ZERO
212+
GAUSS_PT = 0
213+
DO K=1,IORD_K
214+
DO IJ=1,IORD_IJ
215+
GAUSS_PT = GAUSS_PT + 1
209216

210-
DO I=1,ELGP
211-
DO J=1,3
212-
K = 6*(I-1) + J
213-
ME(K,K) = M0
217+
CALL SHP3DP ( I, J, K, ELGP, SUBR_NAME, IORD_MSG, IORD_IJ, IORD_K, SS_I(IJ), SS_J(IJ), SS_K(K), 'N', PSH, DPSHG )
218+
INTFAC = DETJ(GAUSS_PT)*HH_IJ(IJ)*HH_K(K)
219+
220+
DO L=1,ELGP
221+
DO M=1,ELGP
222+
M_1DOF(L,M) = M_1DOF(L,M) + PSH(L) * PSH(M) * RHO(1) * INTFAC
223+
ENDDO
224+
ENDDO
214225
ENDDO
215226
ENDDO
216227

228+
229+
CALL EXPAND_MASS_DOFS( M_1DOF )
230+
217231
ENDIF
218232

219233
! **********************************************************************************************************************************

‎Source/EMG/EMG5/TETRA.f90‎

Lines changed: 21 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ SUBROUTINE TETRA ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, WRITE_WARN )
4040
USE IOUNT1, ONLY : WRT_ERR, WRT_LOG, ERR, F04, F06
4141
USE SCONTR, ONLY : BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_TETRA, NTSUB
4242
USE TIMDAT, ONLY : TSEC
43-
USE CONSTANTS_1, ONLY : HALF, QUARTER, ZERO, FOUR
43+
USE CONSTANTS_1, ONLY : HALF, QUARTER, ZERO
4444
USE DEBUG_PARAMETERS, ONLY : DEBUG
4545
USE SUBR_BEGEND_LEVELS, ONLY : TETRA_BEGEND
4646
USE NONLINEAR_PARAMS, ONLY : LOAD_ISTEP
@@ -49,6 +49,7 @@ SUBROUTINE TETRA ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, WRITE_WARN )
4949
SE1, SE2, STE1, STRESS, TREF, TYPE
5050

5151
USE TETRA_USE_IFs
52+
USE EXPAND_MASS_DOFS_Interface
5253

5354
IMPLICIT NONE
5455

@@ -100,7 +101,6 @@ SUBROUTINE TETRA ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, WRITE_WARN )
100101
REAL(DOUBLE) :: JAC(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian matrix.
101102
REAL(DOUBLE) :: JACI(3,3) ! An output from subr JAC3D, called herein. 3 x 3 Jacobian inverse.
102103
REAL(DOUBLE) :: KWW(3,3) ! Portion of differential stiffness matrix
103-
REAL(DOUBLE) :: M0 ! An intermediate variable used in calc elem mass, ME
104104
REAL(DOUBLE) :: PSH(ELGP) ! Output from subr SHP3DT. Shape fcn at Gauss pts SSI, SSJ
105105
REAL(DOUBLE) :: SIGxx ! Normal stress in the elem x direction
106106
REAL(DOUBLE) :: SIGyy ! Normal stress in the elem y direction
@@ -117,7 +117,8 @@ SUBROUTINE TETRA ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, WRITE_WARN )
117117
REAL(DOUBLE) :: TREF1 ! TREF(1)
118118
REAL(DOUBLE) :: VOLUME ! 3D element volume
119119
REAL(DOUBLE) :: SSI,SSJ,SSK ! Isoparametric coordinates of a point.
120-
120+
REAL(DOUBLE) :: M_1DOF(ELGP,ELGP) ! Consistent mass matrix with 1 DOF per node.
121+
121122
! **********************************************************************************************************************************
122123
IF (WRT_LOG >= SUBR_BEGEND) THEN
123124
CALL OURTIM
@@ -190,15 +191,27 @@ SUBROUTINE TETRA ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, WRITE_WARN )
190191

191192
IF (OPT(1) == 'Y') THEN
192193

193-
M0 = (RHO(1))*VOLUME/ELGP
194+
! Consistent mass matrix
195+
! ME = ∫ N' ρ N det(J) dv
196+
197+
M_1DOF(:,:) = ZERO
198+
GAUSS_PT = 0
199+
DO I=1,IORD
200+
GAUSS_PT = GAUSS_PT + 1
201+
202+
CALL SHP3DT ( I, ELGP, SUBR_NAME, IORD_MSG, IORD, SSS_I(I), SSS_J(I), SSS_K(I), 'N', PSH, DPSHG )
203+
INTFAC = DETJ(I)*HHH_IJK(I)
194204

195-
DO I=1,ELGP
196-
DO J=1,3
197-
K = 6*(I-1) + J
198-
ME(K,K) = M0
205+
DO L=1,ELGP
206+
DO M=1,ELGP
207+
M_1DOF(L,M) = M_1DOF(L,M) + PSH(L) * PSH(M) * RHO(1) * INTFAC
208+
ENDDO
199209
ENDDO
200210
ENDDO
201211

212+
213+
CALL EXPAND_MASS_DOFS( M_1DOF )
214+
202215
ENDIF
203216

204217
! **********************************************************************************************************************************

‎Source/EMG/EXPAND_MASS_DOFS.f90‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
! #################################################################################################################################
2+
! Begin MIT license text.
3+
! _______________________________________________________________________________________________________
4+
5+
! Copyright 2022 Dr William R Case, Jr (mystransolver@gmail.com)
6+
7+
! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
8+
! associated documentation files (the "Software"), to deal in the Software without restriction, including
9+
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10+
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
11+
! the following conditions:
12+
13+
! The above copyright notice and this permission notice shall be included in all copies or substantial
14+
! portions of the Software and documentation.
15+
16+
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
17+
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18+
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19+
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20+
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21+
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22+
! THE SOFTWARE.
23+
! _______________________________________________________________________________________________________
24+
25+
! End MIT license text.
26+
SUBROUTINE EXPAND_MASS_DOFS ( M_1DOF )
27+
28+
! Copies per-node mass values to all 3 of each node's translational DOFs in the element mass matrix ME.
29+
30+
USE PENTIUM_II_KIND, ONLY : DOUBLE, LONG
31+
USE MODEL_STUF, ONLY : ELGP, ME
32+
USE CONSTANTS_1, ONLY : ZERO
33+
USE SCONTR, ONLY : SOL_NAME
34+
35+
IMPLICIT NONE
36+
37+
INTEGER(LONG) :: I,J,L
38+
INTEGER(LONG) :: KI, KJ ! For converting grid point number to element DOF number
39+
40+
REAL(DOUBLE), INTENT(INOUT) :: M_1DOF(ELGP,ELGP)
41+
REAL(DOUBLE) :: ROW_SUM
42+
43+
! **********************************************************************************************************************************
44+
45+
! SOL 101, 104, 105, 31(?) must use lumped mass matrix because GRAV and RFORCE assume 6x6 block diagonal structure for MGG.
46+
! SOL 103 should use consistent mass matrix because quadratic elements won't solve with this row sum mass lumping. They
47+
! will solve with equal-mass-per-node lumping but with slightly worse results.
48+
IF ( SOL_NAME(1:5) /= 'MODES') THEN
49+
50+
! Row sum to convert consistent to lumped mass matrix
51+
DO I=1,ELGP
52+
ROW_SUM = ZERO
53+
DO J=1,ELGP
54+
ROW_SUM = ROW_SUM + M_1DOF(I,J)
55+
M_1DOF(I,J) = ZERO
56+
ENDDO
57+
M_1DOF(I,I) = ROW_SUM
58+
ENDDO
59+
60+
ENDIF
61+
62+
ME(:,:) = ZERO
63+
64+
! Copy to all 3 translational DOFs of each node.
65+
DO I=1,ELGP
66+
DO J=1,ELGP
67+
KI = (I-1) * 6
68+
KJ = (J-1) * 6
69+
DO L=1,3
70+
ME(KI + L, KJ + L) = M_1DOF(I,J)
71+
ENDDO
72+
ENDDO
73+
ENDDO
74+
75+
RETURN
76+
77+
78+
! **********************************************************************************************************************************
79+
80+
END SUBROUTINE EXPAND_MASS_DOFS

0 commit comments

Comments
 (0)