FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_3dIC.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
10
11 use hecmw, only : kint, kreal
12 use m_utilities
13 use elementinfo
14
15 implicit none
16
17contains
18
20 !----------------------------------------------------------------------*
21 subroutine stf_c3d8ic &
22 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
23 time, tincr, u, aux, temperature)
24 !----------------------------------------------------------------------*
25
26 use mmechgauss
27 use m_matmatrix
29 use m_static_lib_3d, only: geomat_c3
30
31 !---------------------------------------------------------------------
32
33 integer(kind=kint), intent(in) :: etype
34 integer(kind=kint), intent(in) :: nn
35 real(kind=kreal), intent(in) :: ecoord(3, nn)
36 type(tgaussstatus), intent(in) :: gausses(:)
37 real(kind=kreal), intent(out) :: stiff(:, :)
38 integer(kind=kint), intent(in) :: cdsys_id
39 real(kind=kreal), intent(inout) :: coords(3, 3)
40 real(kind=kreal), intent(in) :: time
41 real(kind=kreal), intent(in) :: tincr
42 real(kind=kreal), intent(in), optional :: u(3, nn)
43 real(kind=kreal), intent(in), optional :: aux(3, 3)
44 real(kind=kreal), intent(in), optional :: temperature(nn)
45
46 !---------------------------------------------------------------------
47
48 integer(kind=kint) :: flag
49 integer(kind=kint), parameter :: ndof = 3
50 real(kind=kreal) :: d(6, 6), b(6, ndof*(nn+3)), db(6, ndof*(nn+3))
51 real(kind=kreal) :: gderiv(nn+3, 3), stress(6)
52 real(kind=kreal) :: xj(9, 9), jacobian(3, 3), inverse(3, 3)
53 real(kind=kreal) :: tmpstiff((nn+3)*3, (nn+3)*3), tmpk(nn*3, 9)
54 real(kind=kreal) :: det, wg, elem(3, nn), mat(6, 6)
55 integer(kind=kint) :: i, j, lx
56 integer(kind=kint) :: serr
57 real(kind=kreal) :: naturalcoord(3), unode(3, nn+3)
58 real(kind=kreal) :: gdispderiv(3, 3), coordsys(3, 3)
59 real(kind=kreal) :: b1(6, ndof*(nn+3))
60 real(kind=kreal) :: smat(9, 9)
61 real(kind=kreal) :: bn(9, ndof*(nn+3)), sbn(9, ndof*(nn+3))
62 real(kind=kreal) :: spfunc(nn), temp
63
64 !---------------------------------------------------------------------
65
66 if( present(u) .AND. present(aux) ) then
67 unode(:, 1:nn) = u(:, :)
68 unode(:, nn+1:nn+3) = aux(:, :)
69 end if
70
71 ! we suppose the same material type in the element
72 flag = gausses(1)%pMaterial%nlgeom_flag
73 if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
74 elem(:, :) = ecoord(:, :)
75 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+unode(:, 1:nn)
76
77 ! --- Inverse of Jacobian at elemental center
78 naturalcoord(:) = 0.0d0
79 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
80 inverse(:, :)= inverse(:, :)*det
81 ! ---- We now calculate stiff matrix include imcompatible mode
82 ! [ Kdd Kda ]
83 ! [ Kad Kaa ]
84 tmpstiff(:, :) = 0.0d0
85 b1(1:6, 1:(nn+3)*ndof) = 0.0d0
86 bn(1:9, 1:(nn+3)*ndof) = 0.0d0
87
88 do lx = 1, numofquadpoints(etype)
89
90 call getquadpoint(etype, lx, naturalcoord)
91 call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
92
93 if( cdsys_id > 0 ) then
94 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
95 if( serr == -1 ) stop "Fail to setup local coordinate"
96 if( serr == -2 ) then
97 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
98 end if
99 end if
100
101 if( present(temperature) ) then
102 call getshapefunc( etype, naturalcoord, spfunc )
103 temp = dot_product( temperature, spfunc )
104 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
105 else
106 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
107 end if
108
109 if( flag == updatelag ) then
110 call geomat_c3( gausses(lx)%stress, mat )
111 d(:, :) = d(:, :)-mat
112 endif
113
114 ! -- Derivative of shape function of imcompatible mode --
115 ! [ -2*a 0, 0 ]
116 ! [ 0, -2*b, 0 ]
117 ! [ 0, 0, -2*c ]
118 ! we don't call getShapeDeriv but use above value directly for computation efficiency
119 gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
120 gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
121 gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
122
123 wg = getweight(etype, lx)*det
124 b(1:6, 1:(nn+3)*ndof)=0.0d0
125 do j = 1, nn+3
126 b(1, 3*j-2) = gderiv(j, 1)
127 b(2, 3*j-1) = gderiv(j, 2)
128 b(3, 3*j ) = gderiv(j, 3)
129 b(4, 3*j-2) = gderiv(j, 2)
130 b(4, 3*j-1) = gderiv(j, 1)
131 b(5, 3*j-1) = gderiv(j, 3)
132 b(5, 3*j ) = gderiv(j, 2)
133 b(6, 3*j-2) = gderiv(j, 3)
134 b(6, 3*j ) = gderiv(j, 1)
135 end do
136
137 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
138 if( flag == totallag ) then
139 ! ---dudx(i,j) ==> gdispderiv(i,j)
140 gdispderiv(1:ndof,1:ndof) = matmul( unode(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
141 do j = 1, nn+3
142
143 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
144 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
145 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
146 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
147 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
148 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
149 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
150 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
151 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
152 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
153 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
154 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
155 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
156 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
157 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
158 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
159 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
160 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
161
162 end do
163 ! ---BL = BL0 + BL1
164 do j = 1, (nn+3)*ndof
165 b(:, j) = b(:, j)+b1(:, j)
166 end do
167
168 end if
169
170 db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
171 forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof )
172 tmpstiff(i, j) = tmpstiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
173 end forall
174
175 ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
176 if( flag == totallag .OR. flag == updatelag ) then
177 stress(1:6) = gausses(lx)%stress
178 do j = 1, nn+3
179 bn(1, 3*j-2) = gderiv(j, 1)
180 bn(2, 3*j-1) = gderiv(j, 1)
181 bn(3, 3*j ) = gderiv(j, 1)
182 bn(4, 3*j-2) = gderiv(j, 2)
183 bn(5, 3*j-1) = gderiv(j, 2)
184 bn(6, 3*j ) = gderiv(j, 2)
185 bn(7, 3*j-2) = gderiv(j, 3)
186 bn(8, 3*j-1) = gderiv(j, 3)
187 bn(9, 3*j ) = gderiv(j, 3)
188 end do
189 smat(:, :) = 0.0d0
190 do j = 1, 3
191 smat(j , j ) = stress(1)
192 smat(j , j+3) = stress(4)
193 smat(j , j+6) = stress(6)
194 smat(j+3, j ) = stress(4)
195 smat(j+3, j+3) = stress(2)
196 smat(j+3, j+6) = stress(5)
197 smat(j+6, j ) = stress(6)
198 smat(j+6, j+3) = stress(5)
199 smat(j+6, j+6) = stress(3)
200 end do
201 sbn(1:9, 1:(nn+3)*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:(nn+3)*ndof) )
202 forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof )
203 tmpstiff(i, j) = tmpstiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
204 end forall
205 end if
206
207 end do
208
209 ! -----Condense tmpstiff to stiff
210 xj(1:9, 1:9)= tmpstiff(nn*ndof+1:(nn+3)*ndof, nn*ndof+1:(nn+3)*ndof)
211 call calinverse(9, xj)
212 tmpk = matmul( tmpstiff( 1:nn*ndof, nn*ndof+1:(nn+3)*ndof ), xj )
213 stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)-matmul( tmpk, tmpstiff(nn*ndof+1:(nn+3)*ndof, 1:nn*ndof) )
214
215 end subroutine stf_c3d8ic
216
217
219 !---------------------------------------------------------------------*
220 subroutine update_c3d8ic &
221 (etype,nn,ecoord, u, du, ddu, cdsys_id, coords, &
222 qf, gausses, iter, time, tincr, aux, ddaux, tt, t0, tn )
223 !---------------------------------------------------------------------*
224
225 use m_fstr
226 use mmaterial
227 use mmechgauss
228 use m_matmatrix
230 use m_utilities
232
233 integer(kind=kint), intent(in) :: etype
234 integer(kind=kint), intent(in) :: nn
235 real(kind=kreal), intent(in) :: ecoord(3, nn)
236 real(kind=kreal), intent(in) :: u(3, nn)
237 real(kind=kreal), intent(in) :: du(3, nn)
238 real(kind=kreal), intent(in) :: ddu(3, nn)
239 integer(kind=kint), intent(in) :: cdsys_id
240 real(kind=kreal), intent(inout) :: coords(3, 3)
241 real(kind=kreal), intent(out) :: qf(nn*3)
242 type(tgaussstatus), intent(inout) :: gausses(:)
243 integer, intent(in) :: iter
244 real(kind=kreal), intent(in) :: time
245 real(kind=kreal), intent(in) :: tincr
246 real(kind=kreal), intent(inout) :: aux(3, 3)
247 real(kind=kreal), intent(out) :: ddaux(3, 3)
248 real(kind=kreal), intent(in), optional :: tt(nn)
249 real(kind=kreal), intent(in), optional :: t0(nn)
250 real(kind=kreal), intent(in), optional :: tn(nn)
251
252 ! LCOAL VARIAVLES
253 integer(kind=kint) :: flag
254 integer(kind=kint), parameter :: ndof = 3
255 real(kind=kreal) :: d(6,6), b(6,ndof*(nn+3)), b1(6,ndof*(nn+3)), spfunc(nn), ina(1)
256 real(kind=kreal) :: bn(9,ndof*(nn+3)), db(6, ndof*(nn+3)), stress(6), smat(9, 9), sbn(9, ndof*(nn+3))
257 real(kind=kreal) :: gderiv(nn+3,3), gderiv0(nn+3,3), gdispderiv(3,3), f(3,3), det, det0, wg, ttc, tt0, ttn
258 integer(kind=kint) :: i, j, lx, mtype, serr
259 real(kind=kreal) :: naturalcoord(3), rot(3,3), mat(6,6), epsth(6)
260 real(kind=kreal) :: totaldisp(3,nn+3), elem(3,nn), elem1(3,nn), coordsys(3,3)
261 real(kind=kreal) :: dstrain(6)
262 real(kind=kreal) :: alpo(3)
263 logical :: ierr, matlaniso
264 real(kind=kreal) :: stiff_ad(9, nn*3), stiff_aa(9, 9), qf_a(9)
265 real(kind=kreal) :: xj(9, 9)
266 real(kind=kreal) :: tmpforce(9), tmpdisp(9), tmp_d(nn*3), tmp_a(9)
267 real(kind=kreal) :: jacobian(3, 3), inverse(3, 3), inverse1(3, 3), inverse0(3, 3)
268
269 !---------------------------------------------------------------------
270
271 totaldisp(:, 1:nn) = u(:, :) + du(:, :) - ddu(:, :)
272 totaldisp(:, nn+1:nn+3) = aux(:, :)
273
274 ! we suppose the same material type in the element
275 flag = gausses(1)%pMaterial%nlgeom_flag
276 elem(:, :) = ecoord(:, :)
277 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+totaldisp(:, 1:nn)
278
279 ! --- Inverse of Jacobian at elemental center
280 naturalcoord(:) = 0.0d0
281 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
282 inverse(:, :)= inverse(:, :)*det
283 ! ---- We now calculate stiff matrix include imcompatible mode
284 ! [ Kdd Kda ]
285 ! [ Kad Kaa ]
286 stiff_ad(:, :) = 0.0d0
287 stiff_aa(:, :) = 0.0d0
288 b1(1:6, 1:(nn+3)*ndof) = 0.0d0
289 bn(1:9, 1:(nn+3)*ndof) = 0.0d0
290
291 do lx = 1, numofquadpoints(etype)
292
293 call getquadpoint(etype, lx, naturalcoord)
294 call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
295
296 if( cdsys_id > 0 ) then
297 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
298 if( serr == -1 ) stop "Fail to setup local coordinate"
299 if( serr == -2 ) then
300 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
301 end if
302 end if
303
304 if( present(tt) ) then
305 call getshapefunc( etype, naturalcoord, spfunc )
306 ttc = dot_product( tt, spfunc )
307 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, ttc )
308 else
309 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
310 end if
311
312 if( flag == updatelag ) then
313 call geomat_c3( gausses(lx)%stress, mat )
314 d(:, :) = d(:, :)-mat
315 endif
316
317 ! -- Derivative of shape function of imcompatible mode --
318 ! [ -2*a 0, 0 ]
319 ! [ 0, -2*b, 0 ]
320 ! [ 0, 0, -2*c ]
321 ! we don't call getShapeDeriv but use above value directly for computation efficiency
322 gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
323 gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
324 gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
325
326 wg = getweight(etype, lx)*det
327 b(1:6, 1:(nn+3)*ndof)=0.0d0
328 do j = 1, nn+3
329 b(1, 3*j-2) = gderiv(j, 1)
330 b(2, 3*j-1) = gderiv(j, 2)
331 b(3, 3*j ) = gderiv(j, 3)
332 b(4, 3*j-2) = gderiv(j, 2)
333 b(4, 3*j-1) = gderiv(j, 1)
334 b(5, 3*j-1) = gderiv(j, 3)
335 b(5, 3*j ) = gderiv(j, 2)
336 b(6, 3*j-2) = gderiv(j, 3)
337 b(6, 3*j ) = gderiv(j, 1)
338 end do
339
340 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
341 if( flag == totallag ) then
342 ! ---dudx(i,j) ==> gdispderiv(i,j)
343 gdispderiv(1:ndof,1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
344 do j = 1, nn+3
345
346 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
347 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
348 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
349 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
350 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
351 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
352 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
353 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
354 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
355 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
356 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
357 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
358 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
359 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
360 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
361 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
362 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
363 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
364
365 end do
366 ! ---BL = BL0 + BL1
367 do j = 1, (nn+3)*ndof
368 b(:, j) = b(:, j)+b1(:, j)
369 end do
370
371 end if
372
373 db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
374 forall( i=1:3*ndof, j=1:nn*ndof )
375 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
376 end forall
377 forall( i=1:3*ndof, j=1:3*ndof )
378 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
379 end forall
380
381 ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
382 if( flag == totallag .OR. flag == updatelag ) then
383 stress(1:6) = gausses(lx)%stress
384 do j = 1, nn+3
385 bn(1, 3*j-2) = gderiv(j, 1)
386 bn(2, 3*j-1) = gderiv(j, 1)
387 bn(3, 3*j ) = gderiv(j, 1)
388 bn(4, 3*j-2) = gderiv(j, 2)
389 bn(5, 3*j-1) = gderiv(j, 2)
390 bn(6, 3*j ) = gderiv(j, 2)
391 bn(7, 3*j-2) = gderiv(j, 3)
392 bn(8, 3*j-1) = gderiv(j, 3)
393 bn(9, 3*j ) = gderiv(j, 3)
394 end do
395 smat(:, :) = 0.0d0
396 do j = 1, 3
397 smat(j , j ) = stress(1)
398 smat(j , j+3) = stress(4)
399 smat(j , j+6) = stress(6)
400 smat(j+3, j ) = stress(4)
401 smat(j+3, j+3) = stress(2)
402 smat(j+3, j+6) = stress(5)
403 smat(j+6, j ) = stress(6)
404 smat(j+6, j+3) = stress(5)
405 smat(j+6, j+6) = stress(3)
406 end do
407 sbn(1:9, 1:(nn+3)*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:(nn+3)*ndof) )
408 forall( i=1:3*ndof, j=1:nn*ndof )
409 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, j) )*wg
410 end forall
411 forall( i=1:3*ndof, j=1:3*ndof )
412 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, nn*ndof+j) )*wg
413 end forall
414 end if
415
416 end do
417
418 ! -----recover incompatible dof of nodal displacement
419 xj(1:3*ndof, 1:3*ndof)= stiff_aa(1:3*ndof, 1:3*ndof)
420 call calinverse(3*ndof, xj)
421 ! --- [Kad] * ddu
422 do i=1,nn
423 tmp_d((i-1)*ndof+1:i*ndof) = ddu(1:ndof, i)
424 enddo
425 tmpforce(1:3*ndof) = matmul( stiff_ad(1:3*ndof, 1:nn*ndof), tmp_d(1:nn*ndof) )
426 ! --- ddaux = -[Kaa]-1 * ([Kad] * ddu)
427 tmp_a(1:3*ndof) = -matmul( xj(1:3*ndof, 1:3*ndof), tmpforce(1:3*ndof) )
428 do i=1,3
429 ddaux(1:ndof, i) = tmp_a((i-1)*ndof+1:i*ndof)
430 enddo
431
432
433 !---------------------------------------------------------------------
434
435 totaldisp(:,1:nn) = u(:,:)+ du(:,:)
436 totaldisp(:,nn+1:nn+3) = aux(:,:) + ddaux(:,:)
437
438 !---------------------------------------------------------------------
439
440 qf(:) = 0.0d0
441 qf_a(:) = 0.0d0
442
443 stiff_ad(:, :) = 0.0d0
444 stiff_aa(:, :) = 0.0d0
445
446 !---------------------------------------------------------------------
447
448 if( flag == updatelag ) then
449 elem(:,:) = (0.5d0*du(:,:)+u(:,:) ) +ecoord(:,:)
450 elem1(:,:) = (du(:,:)+u(:,:) ) +ecoord(:,:)
451 ! elem = elem1
452 totaldisp(:,1:nn) = du(:,:)
453 if( iter == 1 ) aux(:,:) = 0.0d0 !--- is this correct???
454 totaldisp(:,nn+1:nn+3) = aux(:,:)
455 end if
456
457 matlaniso = .false.
458 if( present(tt) .and. cdsys_id > 0 ) then
459 ina = tt(1)
460 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
461 if( .not. ierr ) matlaniso = .true.
462 end if
463
464 ! --- Inverse of Jacobian at elemental center
465 naturalcoord(:) = 0.0d0
466 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
467 inverse(:, :)= inverse(:, :)*det
468 if( flag == updatelag ) then
469 call getjacobian(etype, nn, naturalcoord, elem1, det, jacobian, inverse1)
470 inverse1(:, :)= inverse1(:, :)*det
471 call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse0)
472 inverse0(:, :)= inverse0(:, :)*det !for deformation gradient F
473 endif
474
475 do lx = 1, numofquadpoints(etype)
476
477 mtype = gausses(lx)%pMaterial%mtype
478
479 call getquadpoint( etype, lx, naturalcoord(:) )
480 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
481
482 if( cdsys_id > 0 ) then
483 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
484 if( serr == -1 ) stop "Fail to setup local coordinate"
485 if( serr == -2 ) then
486 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
487 end if
488 end if
489
490 ! ========================================================
491 ! UPDATE STRAIN and STRESS
492 ! ========================================================
493
494 ! Thermal Strain
495 epsth = 0.0d0
496 if( present(tt) .AND. present(t0) ) then
497 call getshapefunc(etype, naturalcoord, spfunc)
498 ttc = dot_product(tt, spfunc)
499 tt0 = dot_product(t0, spfunc)
500 ttn = dot_product(tn, spfunc)
501 call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
502 end if
503
504 ! -- Derivative of shape function of imcompatible mode --
505 ! [ -2*a 0, 0 ]
506 ! [ 0, -2*b, 0 ]
507 ! [ 0, 0, -2*c ]
508 ! we don't call getShapeDeriv but use above value directly for computation efficiency
509 gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
510 gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
511 gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
512
513 ! Small strain
514 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
515 dstrain(1) = gdispderiv(1, 1)
516 dstrain(2) = gdispderiv(2, 2)
517 dstrain(3) = gdispderiv(3, 3)
518 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
519 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
520 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
521 dstrain(:) = dstrain(:)-epsth(:) ! allright?
522
523 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
524 if( flag == infinitesimal ) then
525 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
526 f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
527
528 else if( flag == totallag ) then
529 ! Green-Lagrange strain
530 dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
531 dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
532 dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
533 dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
534 +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
535 dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
536 +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
537 dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
538 +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
539
540 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
541
542 else if( flag == updatelag ) then
543 ! CALL GEOMAT_C3( gausses(LX)%stress, mat )
544 ! D(:, :) = D(:, :)+mat(:, :)
545 rot = 0.0d0
546 rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
547 rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
548 rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
549
550 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+ dstrain(1:6)+epsth(:)
551 call getglobalderiv(etype, nn, naturalcoord, ecoord, det0, gderiv0)
552 gderiv0(nn+1, :) = -2.0d0*naturalcoord(1)*inverse0(1, :)/det0
553 gderiv0(nn+2, :) = -2.0d0*naturalcoord(2)*inverse0(2, :)/det0
554 gderiv0(nn+3, :) = -2.0d0*naturalcoord(3)*inverse0(3, :)/det0
555 f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn+3)+du(1:ndof, 1:nn+3), gderiv0(1:nn+3, 1:ndof) )
556
557 end if
558
559 ! Update stress
560 if( present(tt) .AND. present(t0) ) then
561 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
562 else
563 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
564 end if
565
566 ! ========================================================
567 ! calculate the internal force ( equivalent nodal force )
568 ! ========================================================
569 ! Small strain
570 b(1:6, 1:(nn+3)*ndof) = 0.0d0
571 do j=1,nn+3
572 b(1,3*j-2) = gderiv(j, 1)
573 b(2,3*j-1) = gderiv(j, 2)
574 b(3,3*j ) = gderiv(j, 3)
575 b(4,3*j-2) = gderiv(j, 2)
576 b(4,3*j-1) = gderiv(j, 1)
577 b(5,3*j-1) = gderiv(j, 3)
578 b(5,3*j ) = gderiv(j, 2)
579 b(6,3*j-2) = gderiv(j, 3)
580 b(6,3*j ) = gderiv(j, 1)
581 end do
582
583 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
584 if( flag == infinitesimal ) then
585
586 else if( flag == totallag ) then
587
588 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
589 b1(1:6, 1:(nn+3)*ndof)=0.0d0
590 do j = 1,nn+3
591 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
592 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
593 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
594 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
595 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
596 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
597 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
598 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
599 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
600 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
601 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
602 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
603 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
604 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
605 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
606 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
607 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
608 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
609 end do
610 ! BL = BL0 + BL1
611 do j=1,(nn+3)*ndof
612 b(:,j) = b(:,j)+b1(:,j)
613 end do
614
615 else if( flag == updatelag ) then
616
617 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv(1:nn,1:3))
618
619 ! -- Derivative of shape function of imcompatible mode --
620 ! [ -2*a 0, 0 ]
621 ! [ 0, -2*b, 0 ]
622 ! [ 0, 0, -2*c ]
623 ! we don't call getShapeDeriv but use above value directly for computation efficiency
624 gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse1(1, :)/det
625 gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse1(2, :)/det
626 gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse1(3, :)/det
627
628 b(1:6, 1:(nn+3)*ndof) = 0.0d0
629 do j = 1, nn+3
630 b(1, 3*j-2) = gderiv(j, 1)
631 b(2, 3*j-1) = gderiv(j, 2)
632 b(3, 3*j ) = gderiv(j, 3)
633 b(4, 3*j-2) = gderiv(j, 2)
634 b(4, 3*j-1) = gderiv(j, 1)
635 b(5, 3*j-1) = gderiv(j, 3)
636 b(5, 3*j ) = gderiv(j, 2)
637 b(6, 3*j-2) = gderiv(j, 3)
638 b(6, 3*j ) = gderiv(j, 1)
639 end do
640
641 end if
642
643 wg=getweight( etype, lx )*det
644
645 db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
646 forall( i=1:3*ndof, j=1:nn*ndof )
647 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
648 end forall
649 forall( i=1:3*ndof, j=1:3*ndof )
650 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
651 end forall
652
653 ! calculate the Internal Force
654 qf(1:nn*ndof) &
655 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
656 qf_a(1:3*ndof) &
657 = qf_a(1:3*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,nn*ndof+1:(nn+3)*ndof) )*wg
658 end do
659
660 ! condence ( qf - [Kda] * [Kaa]-1 * qf_a )
661 xj(1:9, 1:9)= stiff_aa(1:3*ndof, 1:3*ndof)
662 call calinverse(9, xj)
663 tmpdisp(:) = matmul( xj(:,:), qf_a(:) )
664 do i=1,nn*ndof
665 qf(i) = qf(i) - dot_product( stiff_ad(:,i), tmpdisp(:) )
666 enddo
667 end subroutine update_c3d8ic
668
669
671 !----------------------------------------------------------------------*
672 subroutine tload_c3d8ic &
673 (etype, nn, xx, yy, zz, tt, t0, &
674 gausses, vect, cdsys_id, coords)
675 !----------------------------------------------------------------------*
676
677 use m_fstr
678 use mmechgauss
679 use m_matmatrix
681
682 !---------------------------------------------------------------------
683
684 integer(kind=kint), parameter :: ndof=3
685 integer(kind=kint), intent(in) :: etype
686 integer(kind=kint), intent(in) :: nn
687 real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
688 real(kind=kreal), intent(in) :: tt(nn), t0(nn)
689 type(tgaussstatus), intent(inout) :: gausses(:)
690 real(kind=kreal), intent(out) :: vect(nn*ndof)
691 integer(kind=kint), intent(in) :: cdsys_id
692 real(kind=kreal), intent(inout) :: coords(3, 3)
693
694 !---------------------------------------------------------------------
695
696 real(kind=kreal) :: alp, alp0
697 real(kind=kreal) :: d(6, 6), b(6, ndof*(nn+3)), db_a(6, ndof*3)
698 real(kind=kreal) :: det, wg, ecoord(3, nn)
699 integer(kind=kint) :: i, j, ic, serr
700 real(kind=kreal) :: epsth(6), sgm(6), h(nn), alpo(3), alpo0(3), coordsys(3, 3), xj(9, 9)
701 real(kind=kreal) :: naturalcoord(3), gderiv(nn+3, 3)
702 real(kind=kreal) :: jacobian(3, 3),inverse(3, 3)
703 real(kind=kreal) :: stiff_da(nn*3, 3*3), stiff_aa(3*3, 3*3)
704 real(kind=kreal) :: vect_a(3*3)
705 real(kind=kreal) :: icdisp(9)
706 real(kind=kreal) :: tempc, temp0, thermal_eps, tm(6, 6), outa(1), ina(1)
707 logical :: ierr, matlaniso
708
709 !---------------------------------------------------------------------
710
711 matlaniso = .false.
712 if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
713 ina = tt(1)
714 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
715 if( .not. ierr ) matlaniso = .true.
716 end if
717
718 ecoord(1, :) = xx(:)
719 ecoord(2, :) = yy(:)
720 ecoord(3, :) = zz(:)
721 ! ---- Calculate enhanced displacement at first
722 ! --- Inverse of Jacobian at elemental center
723 naturalcoord(:) = 0.0d0
724 call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse)
725 inverse(:, :) = inverse(:, :)*det
726 ! ---- We now calculate stiff matrix include imcompatible mode
727 ! [ Kdd Kda ]
728 ! [ Kad Kaa ]
729 stiff_da(:, :) = 0.0d0
730 stiff_aa(:, :) = 0.0d0
731 b(1:6, 1:(nn+3)*ndof) = 0.0d0
732 vect(1:nn*ndof) = 0.0d0
733 vect_a(1:3*ndof) = 0.0d0
734
735 do ic = 1, numofquadpoints(etype)
736
737 call getquadpoint(etype, ic, naturalcoord)
738 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv(1:nn, 1:3) )
739
740 if( cdsys_id > 0 ) then
741 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
742 if( serr == -1 ) stop "Fail to setup local coordinate"
743 if( serr == -2 ) then
744 write(*,*) "WARNING! Cannot setup local coordinate, it is modified automatically"
745 end if
746 end if
747
748 call getshapefunc( etype, naturalcoord, h(1:nn) )
749 tempc = dot_product( h(1:nn), tt(1:nn) )
750 temp0 = dot_product( h(1:nn), t0(1:nn) )
751 call matlmatrix( gausses(ic), d3, d, 1.d0, 1.0d0, coordsys, tempc )
752
753 ! -- Derivative of shape function of imcompatible mode --
754 ! [ -2*a 0, 0 ]
755 ! [ 0, -2*b, 0 ]
756 ! [ 0, 0, -2*c ]
757 ! we don't call getShapeDeriv but use above value directly for computation efficiency
758 gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
759 gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
760 gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
761
762 wg = getweight(etype, ic)*det
763 do j = 1, nn+3
764 b(1, 3*j-2) = gderiv(j, 1)
765 b(2, 3*j-1) = gderiv(j, 2)
766 b(3, 3*j ) = gderiv(j, 3)
767 b(4, 3*j-2) = gderiv(j, 2)
768 b(4, 3*j-1) = gderiv(j, 1)
769 b(5, 3*j-1) = gderiv(j, 3)
770 b(5, 3*j ) = gderiv(j, 2)
771 b(6, 3*j-2) = gderiv(j, 3)
772 b(6, 3*j ) = gderiv(j, 1)
773 end do
774
775 db_a(1:6, 1:3*ndof) = matmul( d, b(1:6, nn*ndof+1:(nn+3)*ndof) )
776 forall( i=1:nn*ndof, j=1:3*ndof )
777 stiff_da(i, j) = stiff_da(i, j) + dot_product( b(1:6, i), db_a(1:6, j) ) * wg
778 end forall
779 forall( i=1:3*ndof, j=1:3*ndof )
780 stiff_aa(i, j) = stiff_aa(i, j) + dot_product( b(1:6, nn*ndof+i), db_a(1:6, j) ) * wg
781 end forall
782
783 ina(1) = tempc
784 if( matlaniso ) then
785 call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo(:), ierr, ina )
786 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
787 else
788 call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
789 if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
790 alp = outa(1)
791 end if
792 ina(1) = temp0
793 if( matlaniso ) then
794 call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo0(:), ierr, ina )
795 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
796 else
797 call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
798 if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
799 alp0 = outa(1)
800 end if
801
802 !**
803 !** THERMAL strain
804 !**
805 if( matlaniso ) then
806 do j = 1, 3
807 epsth(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
808 end do
809 epsth(4:6) = 0.0d0
810 call transformation(coordsys, tm)
811 epsth(:) = matmul( epsth(:), tm ) ! to global coord
812 epsth(4:6) = epsth(4:6)*2.0d0
813 else
814 thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
815 epsth(1:3) = thermal_eps
816 epsth(4:6) = 0.0d0
817 end if
818
819 !**
820 !** SET SGM {s}=[D]{e}
821 !**
822 sgm(:) = matmul( d(:, :), epsth(:) )
823
824 !**
825 !** CALCULATE LOAD {F}=[B]T{e}
826 !**
827 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(1:6), b(1:6, 1:nn*ndof) )*wg
828 vect_a(1:3*ndof) = vect_a(1:3*ndof)+matmul( sgm(1:6), b(1:6, nn*ndof+1:(nn+3)*ndof) )*wg
829
830 end do
831
832 ! --- condense vect
833 xj(1:9,1:9)= stiff_aa(1:9, 1:9)
834 call calinverse(9, xj)
835 ! --- [Kaa]-1 * fa
836 icdisp(1:9) = matmul( xj(:, :), vect_a(1:3*ndof) )
837 ! --- fd - [Kda] * [Kaa]-1 * fa
838 vect(1:nn*ndof) = vect(1:nn*ndof) - matmul( stiff_da(1:nn*ndof, 1:3*ndof), icdisp(1:3*ndof) )
839
840 end subroutine tload_c3d8ic
841
842
843end module m_static_lib_3dic
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getjacobian(fetype, nn, localcoord, elecoord, det, jacobian, inverse)
calculate Jacobian matrix, its determinant and inverse
Definition: element.f90:813
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, secttype, matrix, time, dtime, cdsys, temperature, isep)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine update_stress3d(flag, gauss, rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn)
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, epsth)
Eight-node hexagonal element with imcompatible mode.
subroutine update_c3d8ic(etype, nn, ecoord, u, du, ddu, cdsys_id, coords, qf, gausses, iter, time, tincr, aux, ddaux, tt, t0, tn)
Update strain and stress inside element.
subroutine stf_c3d8ic(etype, nn, ecoord, gausses, stiff, cdsys_id, coords, time, tincr, u, aux, temperature)
CALCULATION STIFF Matrix for C3D8IC ELEMENT.
subroutine tload_c3d8ic(etype, nn, xx, yy, zz, tt, t0, gausses, vect, cdsys_id, coords)
This subroutine calculatess thermal loading.
This module provides aux functions.
Definition: utilities.f90:6
subroutine calinverse(nn, a)
calculate inverse of matrix a
Definition: utilities.f90:258
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13