9 use hecmw,
only : kint, kreal
14 real(kind=kreal),
private,
parameter :: i33(3,3) = &
15 & reshape( (/1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/) )
23 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
24 time, tincr, u, temperature)
35 integer(kind=kint),
intent(in) :: etype
36 integer(kind=kint),
intent(in) :: nn
37 real(kind=kreal),
intent(in) :: ecoord(3, nn)
39 real(kind=kreal),
intent(out) :: stiff(:,:)
40 integer(kind=kint),
intent(in) :: cdsys_id
41 real(kind=kreal),
intent(inout) :: coords(3, 3)
42 real(kind=kreal),
intent(in) :: time
43 real(kind=kreal),
intent(in) :: tincr
44 real(kind=kreal),
intent(in),
optional :: u(:, :)
45 real(kind=kreal),
intent(in),
optional :: temperature(nn)
49 integer(kind=kint) :: flag
50 integer(kind=kint),
parameter :: ndof = 3
51 real(kind=kreal) :: d(6, 6),b(6, ndof*nn),db(6, ndof*nn)
52 real(kind=kreal) :: gderiv(nn, 3),stress(6),mat(6, 6)
53 real(kind=kreal) :: det, wg, temp, spfunc(nn)
54 integer(kind=kint) :: i, j, p, q, lx, serr
55 real(kind=kreal) :: naturalcoord(3)
56 real(kind=kreal) :: gdispderiv(3, 3)
57 real(kind=kreal) :: b1(6, ndof*nn)
58 real(kind=kreal) :: smat(9, 9), elem(3, nn)
59 real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
60 real(kind=kreal) :: coordsys(3, 3)
62 real(kind=kreal) :: elem0(3,nn), elem1(3, nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3,2)
63 real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
64 real(kind=kreal) :: gderiv2_ave(ndof*nn,ndof*nn)
65 real(kind=kreal) :: fbar(3,3), jratio(8), coeff, sff, dstrain(6), ddfs, fs(3,3), gfs(3,2)
71 flag = gausses(1)%pMaterial%nlgeom_flag
72 if( .not.
present(u) ) flag = infinitesimal
73 elem(:, :) = ecoord(:, :)
74 elem0(:, :) = ecoord(:, :)
75 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
76 elem1(:, :) = ecoord(:, :)+u(:, :)
81 gderiv1_ave(1:nn,1:ndof) = 0.d0
82 gderiv2_ave(1:ndof*nn,1:ndof*nn) = 0.d0
87 if( flag == infinitesimal )
then
89 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv(1:nn, 1:ndof)
91 gdispderiv(1:3, 1:3) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
92 jacob =
determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
93 jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob)
95 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
100 gderiv2_ave(3*p-3+i,3*q-3+j) = gderiv2_ave(3*p-3+i,3*q-3+j) &
101 & + jacob*wg*(gderiv1(p,i)*gderiv1(q,j)-gderiv1(q,i)*gderiv1(p,j))
108 jacob_ave = jacob_ave + jacob*wg
110 if( dabs(v0) > 1.d-12 )
then
111 if( dabs(jacob_ave) < 1.d-12 ) stop
'Error in Update_C3D8Fbar: too small average jacob'
112 jacob_ave = jacob_ave/v0
113 jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8)
114 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
115 gderiv2_ave(1:ndof*nn,1:ndof*nn) = gderiv2_ave(1:ndof*nn,1:ndof*nn)/(v0*jacob_ave)
117 stop
'Error in Update_C3D8Fbar: too small element volume'
125 if( cdsys_id > 0 )
then
127 if( serr == -1 ) stop
"Fail to setup local coordinate"
128 if( serr == -2 )
then
129 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
133 if(
present(temperature) )
then
135 temp = dot_product( temperature, spfunc )
136 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
138 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
141 if( flag == updatelag )
then
142 call geomat_c3( gausses(lx)%stress, mat )
143 d(:, :) = d(:, :)-mat
147 b(1:6, 1:nn*ndof) = 0.0d0
149 b(1, 3*j-2) = gderiv(j, 1)
150 b(2, 3*j-1) = gderiv(j, 2)
151 b(3, 3*j ) = gderiv(j, 3)
152 b(4, 3*j-2) = gderiv(j, 2)
153 b(4, 3*j-1) = gderiv(j, 1)
154 b(5, 3*j-1) = gderiv(j, 3)
155 b(5, 3*j ) = gderiv(j, 2)
156 b(6, 3*j-2) = gderiv(j, 3)
157 b(6, 3*j ) = gderiv(j, 1)
161 if( flag == infinitesimal )
then
162 b2(1:6, 1:nn*ndof) = 0.0d0
164 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
165 b2(1,3*j-2:3*j) = z1(1:3,1)
166 b2(2,3*j-2:3*j) = z1(1:3,1)
167 b2(3,3*j-2:3*j) = z1(1:3,1)
172 b(1:3, j) = b(1:3,j)+b2(1:3,j)
175 else if( flag == totallag )
then
176 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
177 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
178 call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
181 b1(1:6, 1:nn*ndof) = 0.0d0
183 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
184 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
185 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
186 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
187 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
188 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
189 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
190 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
191 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
192 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
193 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
194 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
195 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
196 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
197 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
198 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
199 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
200 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
203 dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
204 dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
205 dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
206 dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
207 dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
208 dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
210 b2(1:6, 1:nn*ndof) = 0.0d0
212 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
213 b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3,1)
214 b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3,1)
215 b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3,1)
216 b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3,1)
217 b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3,1)
218 b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3,1)
223 b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
226 else if( flag == updatelag )
then
227 wg = (jratio(lx)**3.d0)*
getweight(etype, lx)*det
229 b2(1:3, 1:nn*ndof) = 0.0d0
231 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
232 b2(1, 3*j-2:3*j) = z1(1:3,1)
233 b2(2, 3*j-2:3*j) = z1(1:3,1)
234 b2(3, 3*j-2:3*j) = z1(1:3,1)
238 b(1:3, j) = b(1:3,j)+b2(1:3,j)
243 db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
244 forall( i=1:nn*ndof, j=1:nn*ndof )
245 stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
249 if( flag == totallag .or. flag == updatelag )
then
250 stress(1:6) = gausses(lx)%stress
251 if( flag == totallag )
then
253 sff = dot_product(stress(1:6),dstrain(1:6))
254 else if( flag == updatelag )
then
256 sff = stress(1)+stress(2)+stress(3)
258 fbar(1:3,1:3) = i33(1:3,1:3)
261 bn(1:9, 1:nn*ndof) = 0.0d0
263 bn(1, 3*j-2) = coeff*gderiv(j, 1)
264 bn(2, 3*j-1) = coeff*gderiv(j, 1)
265 bn(3, 3*j ) = coeff*gderiv(j, 1)
266 bn(4, 3*j-2) = coeff*gderiv(j, 2)
267 bn(5, 3*j-1) = coeff*gderiv(j, 2)
268 bn(6, 3*j ) = coeff*gderiv(j, 2)
269 bn(7, 3*j-2) = coeff*gderiv(j, 3)
270 bn(8, 3*j-1) = coeff*gderiv(j, 3)
271 bn(9, 3*j ) = coeff*gderiv(j, 3)
272 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
273 bn(1, 3*j-2:3*j) = bn(1, 3*j-2:3*j) + fbar(1,1)*z1(1:3,1)
274 bn(2, 3*j-2:3*j) = bn(2, 3*j-2:3*j) + fbar(2,1)*z1(1:3,1)
275 bn(3, 3*j-2:3*j) = bn(3, 3*j-2:3*j) + fbar(3,1)*z1(1:3,1)
276 bn(4, 3*j-2:3*j) = bn(4, 3*j-2:3*j) + fbar(1,2)*z1(1:3,1)
277 bn(5, 3*j-2:3*j) = bn(5, 3*j-2:3*j) + fbar(2,2)*z1(1:3,1)
278 bn(6, 3*j-2:3*j) = bn(6, 3*j-2:3*j) + fbar(3,2)*z1(1:3,1)
279 bn(7, 3*j-2:3*j) = bn(7, 3*j-2:3*j) + fbar(1,3)*z1(1:3,1)
280 bn(8, 3*j-2:3*j) = bn(8, 3*j-2:3*j) + fbar(2,3)*z1(1:3,1)
281 bn(9, 3*j-2:3*j) = bn(9, 3*j-2:3*j) + fbar(3,3)*z1(1:3,1)
285 smat(j , j ) = stress(1)
286 smat(j , j+3) = stress(4)
287 smat(j , j+6) = stress(6)
288 smat(j+3, j ) = stress(4)
289 smat(j+3, j+3) = stress(2)
290 smat(j+3, j+6) = stress(5)
291 smat(j+6, j ) = stress(6)
292 smat(j+6, j+3) = stress(5)
293 smat(j+6, j+6) = stress(3)
295 sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
296 forall( i=1:nn*ndof, j=1:nn*ndof )
297 stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
301 fs(1,1) = fbar(1,1)*stress(1)+fbar(1,2)*stress(4)+fbar(1,3)*stress(6)
302 fs(1,2) = fbar(1,1)*stress(4)+fbar(1,2)*stress(2)+fbar(1,3)*stress(5)
303 fs(1,3) = fbar(1,1)*stress(6)+fbar(1,2)*stress(5)+fbar(1,3)*stress(3)
304 fs(2,1) = fbar(2,1)*stress(1)+fbar(2,2)*stress(4)+fbar(2,3)*stress(6)
305 fs(2,2) = fbar(2,1)*stress(4)+fbar(2,2)*stress(2)+fbar(2,3)*stress(5)
306 fs(2,3) = fbar(2,1)*stress(6)+fbar(2,2)*stress(5)+fbar(2,3)*stress(3)
307 fs(3,1) = fbar(3,1)*stress(1)+fbar(3,2)*stress(4)+fbar(3,3)*stress(6)
308 fs(3,2) = fbar(3,1)*stress(4)+fbar(3,2)*stress(2)+fbar(3,3)*stress(5)
309 fs(3,3) = fbar(3,1)*stress(6)+fbar(3,2)*stress(5)+fbar(3,3)*stress(3)
311 z1(1:3,1) = (gderiv1_ave(i,1:3)-gderiv1(i,1:3))/3.d0
312 gfs(1:3,1) = coeff*matmul(fs,gderiv(i,1:3))
314 z1(1:3,2) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
315 gfs(1:3,2) = coeff*matmul(fs,gderiv(j,1:3))
318 ddfs = z1(p,1)*z1(q,2)
319 ddfs = ddfs + (gderiv2_ave(3*i-3+p,3*j-3+q)-gderiv1_ave(i,p)*gderiv1_ave(j,q))/3.d0
320 ddfs = ddfs + gderiv1(i,q)*gderiv1(j,p)/3.d0
321 ddfs = sff*ddfs + z1(p,1)*gfs(q,2)+z1(q,2)*gfs(p,1)
322 stiff(3*i-3+p, 3*j-3+q) = stiff(3*i-3+p, 3*j-3+q) + ddfs*wg
338 (etype, nn, ecoord, u, du, cdsys_id, coords, &
339 qf ,gausses, iter, time, tincr, tt,t0, tn )
353 integer(kind=kint),
intent(in) :: etype
354 integer(kind=kint),
intent(in) :: nn
355 real(kind=kreal),
intent(in) :: ecoord(3, nn)
356 real(kind=kreal),
intent(in) :: u(3, nn)
357 real(kind=kreal),
intent(in) :: du(3, nn)
358 integer(kind=kint),
intent(in) :: cdsys_id
359 real(kind=kreal),
intent(inout) :: coords(3, 3)
360 real(kind=kreal),
intent(out) :: qf(nn*3)
362 integer,
intent(in) :: iter
363 real(kind=kreal),
intent(in) :: time
364 real(kind=kreal),
intent(in) :: tincr
365 real(kind=kreal),
intent(in),
optional :: tt(nn)
366 real(kind=kreal),
intent(in),
optional :: t0(nn)
367 real(kind=kreal),
intent(in),
optional :: tn(nn)
371 integer(kind=kint) :: flag
372 integer(kind=kint),
parameter :: ndof = 3
373 real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
374 real(kind=kreal) :: gderiv(nn, 3), gdispderiv(3, 3), det, wg
375 integer(kind=kint) :: j, lx, serr
376 real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
377 real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
378 real(kind=kreal) :: dstrain(6)
379 real(kind=kreal) :: dvol
380 real(kind=kreal) :: ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
381 logical :: ierr, matlaniso
383 real(kind=kreal) :: elem0(3,nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3)
384 real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
385 real(kind=kreal) :: fbar(3,3), jratio(8)
386 real(kind=kreal) :: jacob_ave05, gderiv05_ave(nn,ndof)
392 flag = gausses(1)%pMaterial%nlgeom_flag
393 elem0(1:ndof,1:nn) = ecoord(1:ndof,1:nn)
394 totaldisp(:, :) = u(:, :)+du(:, :)
396 elem(:, :) = ecoord(:, :)
397 elem1(:, :) = ecoord(:, :)
399 elem(:, :) = ecoord(:, :)
400 elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
402 elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
403 elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
404 totaldisp(:, :) = du(:, :)
408 if( cdsys_id > 0 .AND.
present(tt) )
then
410 call fetch_tabledata(
mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
411 if( .not. ierr ) matlaniso = .true.
417 gderiv1_ave(1:nn,1:ndof) = 0.d0
420 gderiv05_ave(1:nn,1:ndof) = 0.d0
428 gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
430 gdispderiv(1:3, 1:3) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
431 jacob =
determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
432 jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob)
434 call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
437 jacob_ave = jacob_ave + jacob*wg
438 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
442 jacob_ave05 = jacob_ave05 + wg
443 gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof) + wg*gderiv1(1:nn, 1:ndof)
446 if( dabs(v0) > 1.d-12 )
then
447 if( dabs(jacob_ave) < 1.d-12 ) stop
'Error in Update_C3D8Fbar: too small average jacob'
448 jacob_ave = jacob_ave/v0
449 jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8)
450 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
451 if( flag ==
updatelag ) gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof)/jacob_ave05
453 stop
'Error in Update_C3D8Fbar: too small element volume'
461 if( cdsys_id > 0 )
then
462 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
463 if( serr == -1 ) stop
"Fail to setup local coordinate"
464 if( serr == -2 )
then
465 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
469 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
477 if(
present(tt) .AND.
present(t0) )
then
479 ttc = dot_product(tt, spfunc)
480 tt0 = dot_product(t0, spfunc)
481 ttn = dot_product(tn, spfunc)
487 dvol = dot_product( totaldisp(1,1:nn), gderiv1_ave(1:nn,1) )
488 dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv1_ave(1:nn,2) )
489 dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv1_ave(1:nn,3) )
490 dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
491 dstrain(1) = gdispderiv(1, 1) + dvol
492 dstrain(2) = gdispderiv(2, 2) + dvol
493 dstrain(3) = gdispderiv(3, 3) + dvol
494 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
495 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
496 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
497 dstrain(:) = dstrain(:)-epsth(:)
499 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
502 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
505 dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
506 dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
507 dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
508 dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
509 dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
510 dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
511 dstrain(:) = dstrain(:)-epsth(:)
513 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
516 dvol = dot_product( totaldisp(1,1:nn), gderiv05_ave(1:nn,1) )
517 dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv05_ave(1:nn,2) )
518 dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv05_ave(1:nn,3) )
519 dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
520 dstrain(1) = gdispderiv(1, 1) + dvol
521 dstrain(2) = gdispderiv(2, 2) + dvol
522 dstrain(3) = gdispderiv(3, 3) + dvol
523 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
524 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
525 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
526 dstrain(:) = dstrain(:)-epsth(:)
529 rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
530 rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
531 rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
533 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
535 call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv1)
536 gdispderiv(1:ndof, 1:ndof) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
537 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
542 if(
present(tt) .AND.
present(t0) )
then
543 call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr, ttc, tt0, ttn )
545 call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr )
552 b(1:6, 1:nn*ndof) = 0.0d0
554 b(1,3*j-2) = gderiv(j, 1)
555 b(2,3*j-1) = gderiv(j, 2)
556 b(3,3*j ) = gderiv(j, 3)
557 b(4,3*j-2) = gderiv(j, 2)
558 b(4,3*j-1) = gderiv(j, 1)
559 b(5,3*j-1) = gderiv(j, 3)
560 b(5,3*j ) = gderiv(j, 2)
561 b(6,3*j-2) = gderiv(j, 3)
562 b(6,3*j ) = gderiv(j, 1)
567 gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
569 b2(1:6, 1:nn*ndof) = 0.0d0
571 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
572 b2(1,3*j-2:3*j) = z1(1:3)
573 b2(2,3*j-2:3*j) = z1(1:3)
574 b2(3,3*j-2:3*j) = z1(1:3)
579 b(:, j) = b(:,j)+b2(:,j)
584 b1(1:6, 1:nn*ndof) = 0.0d0
586 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
587 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
588 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
589 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
590 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
591 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
592 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
593 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
594 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
595 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
596 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
597 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
598 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
599 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
600 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
601 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
602 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
603 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
608 b2(1:6, 1:nn*ndof) = 0.0d0
610 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
611 b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3)
612 b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3)
613 b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3)
614 b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3)
615 b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3)
616 b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3)
621 b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
627 wg = (jratio(lx)**3.d0)*
getweight(etype, lx)*det
628 b(1:6, 1:nn*ndof) = 0.0d0
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)
641 b2(1:3, 1:nn*ndof) = 0.0d0
643 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
644 b2(1, 3*j-2:3*j) = z1(1:3)
645 b2(2, 3*j-2:3*j) = z1(1:3)
646 b2(3, 3*j-2:3*j) = z1(1:3)
650 b(1:3, j) = b(1:3,j)+b2(1:3,j)
656 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
665 (etype, nn, xx, yy, zz, tt, t0, &
666 gausses, vect, cdsys_id, coords)
676 integer(kind=kint),
parameter :: ndof = 3
677 integer(kind=kint),
intent(in) :: etype, nn
679 real(kind=kreal),
intent(in) :: xx(nn), yy(nn), zz(nn)
680 real(kind=kreal),
intent(in) :: tt(nn), t0(nn)
681 real(kind=kreal),
intent(out) :: vect(nn*ndof)
682 integer(kind=kint),
intent(in) :: cdsys_ID
683 real(kind=kreal),
intent(inout) :: coords(3, 3)
687 real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
688 real(kind=kreal) :: z1(3), det, ecoord(3, nn)
689 integer(kind=kint) :: j, LX, serr
690 real(kind=kreal) :: estrain(6), sgm(6), h(nn)
691 real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
692 real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
693 real(kind=kreal) :: tempc, temp0, v0, jacob_ave, thermal_eps, tm(6,6)
694 logical :: ierr, matlaniso
700 if( cdsys_id > 0 )
then
702 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
703 if( .not. ierr ) matlaniso = .true.
715 gderiv1_ave(1:nn,1:ndof) = 0.d0
718 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv)
721 jacob_ave = jacob_ave + wg
722 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + wg*gderiv(1:nn, 1:ndof)
724 if( dabs(v0) > 1.d-12 )
then
725 if( dabs(jacob_ave) < 1.d-12 ) stop
'Error in TLOAD_C3D8Fbar: too small average jacob'
726 jacob_ave = jacob_ave/v0
727 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
729 stop
'Error in TLOAD_C3D8Fbar: too small element volume'
740 call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
741 if( serr == -1 ) stop
"Fail to setup local coordinate"
742 if( serr == -2 )
then
743 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
749 b(1:6, 1:nn*ndof) = 0.0d0
751 b(1,3*j-2) = gderiv(j, 1)
752 b(2,3*j-1) = gderiv(j, 2)
753 b(3,3*j ) = gderiv(j, 3)
754 b(4,3*j-2) = gderiv(j, 2)
755 b(4,3*j-1) = gderiv(j, 1)
756 b(5,3*j-1) = gderiv(j, 3)
757 b(5,3*j ) = gderiv(j, 2)
758 b(6,3*j-2) = gderiv(j, 3)
759 b(6,3*j ) = gderiv(j, 1)
763 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
764 b(1,3*j-2:3*j) = b(1,3*j-2:3*j)+z1(1:3)
765 b(2,3*j-2:3*j) = b(2,3*j-2:3*j)+z1(1:3)
766 b(3,3*j-2:3*j) = b(3,3*j-2:3*j)+z1(1:3)
769 tempc = dot_product( h(1:nn),tt(1:nn) )
770 temp0 = dot_product( h(1:nn),t0(1:nn) )
772 call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
776 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
777 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
779 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
780 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
785 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
786 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
788 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
789 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
802 estrain(:) = matmul( estrain(:), tm )
803 estrain(4:6) = estrain(4:6)*2.0d0
806 estrain(1:3) = thermal_eps
813 sgm(:) = matmul( d(:, :), estrain(:) )
818 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
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.
real(kind=kreal), pointer ref_temp
REFTEMP.
This module manages calculation relates with materials.
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn)
This module contains several strategy to free locking problem in Eight-node hexagonal element.
subroutine update_c3d8fbar(etype, nn, ecoord, u, du, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update Strain stress of this element.
subroutine stf_c3d8fbar(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using F-bar method.
subroutine tload_c3d8fbar(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
This module provides aux functions.
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
subroutine transformation(jacob, tm)
This module provides functions for hyperelastic calculation.
This module summarizes all infomation of material properties.
integer(kind=kint), parameter totallag
character(len=dict_key_length) mc_orthoexp
integer(kind=kint), parameter infinitesimal
integer(kind=kint), parameter updatelag
This modules defines a structure to record history dependent parameter in static analysis.
All data should be recorded in every quadrature points.