11 use hecmw,
only : kint, kreal
22 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
23 time, tincr, u, aux, temperature)
33 integer(kind=kint),
intent(in) :: etype
34 integer(kind=kint),
intent(in) :: nn
35 real(kind=kreal),
intent(in) :: ecoord(3, nn)
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)
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
66 if(
present(u) .AND.
present(aux) )
then
67 unode(:, 1:nn) = u(:, :)
68 unode(:, nn+1:nn+3) = aux(:, :)
72 flag = gausses(1)%pMaterial%nlgeom_flag
73 if( .not.
present(u) ) flag = infinitesimal
74 elem(:, :) = ecoord(:, :)
75 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+unode(:, 1:nn)
78 naturalcoord(:) = 0.0d0
79 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
80 inverse(:, :)= inverse(:, :)*det
84 tmpstiff(:, :) = 0.0d0
85 b1(1:6, 1:(nn+3)*ndof) = 0.0d0
86 bn(1:9, 1:(nn+3)*ndof) = 0.0d0
91 call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
93 if( cdsys_id > 0 )
then
95 if( serr == -1 ) stop
"Fail to setup local coordinate"
97 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
101 if(
present(temperature) )
then
103 temp = dot_product( temperature, spfunc )
104 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
106 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
109 if( flag == updatelag )
then
110 call geomat_c3( gausses(lx)%stress, mat )
111 d(:, :) = d(:, :)-mat
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
124 b(1:6, 1:(nn+3)*ndof)=0.0d0
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)
138 if( flag == totallag )
then
140 gdispderiv(1:ndof,1:ndof) = matmul( unode(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
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)
164 do j = 1, (nn+3)*ndof
165 b(:, j) = b(:, j)+b1(:, j)
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
176 if( flag == totallag .OR. flag == updatelag )
then
177 stress(1:6) = gausses(lx)%stress
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)
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)
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
210 xj(1:9, 1:9)= tmpstiff(nn*ndof+1:(nn+3)*ndof, nn*ndof+1:(nn+3)*ndof)
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) )
221 (etype,nn,ecoord, u, du, ddu, cdsys_id, coords, &
222 qf, gausses, iter, time, tincr, aux, ddaux, tt, t0, tn )
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)
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)
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)
271 totaldisp(:, 1:nn) = u(:, :) + du(:, :) - ddu(:, :)
272 totaldisp(:, nn+1:nn+3) = aux(:, :)
275 flag = gausses(1)%pMaterial%nlgeom_flag
276 elem(:, :) = ecoord(:, :)
277 if( flag ==
updatelag ) elem(:, :) = ecoord(:, :)+totaldisp(:, 1:nn)
280 naturalcoord(:) = 0.0d0
281 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
282 inverse(:, :)= inverse(:, :)*det
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
294 call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
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"
304 if(
present(tt) )
then
306 ttc = dot_product( tt, spfunc )
307 call matlmatrix( gausses(lx),
d3, d, time, tincr, coordsys, ttc )
309 call matlmatrix( gausses(lx),
d3, d, time, tincr, coordsys )
313 call geomat_c3( gausses(lx)%stress, mat )
314 d(:, :) = d(:, :)-mat
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
327 b(1:6, 1:(nn+3)*ndof)=0.0d0
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)
343 gdispderiv(1:ndof,1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
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)
367 do j = 1, (nn+3)*ndof
368 b(:, j) = b(:, j)+b1(:, j)
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
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
383 stress(1:6) = gausses(lx)%stress
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)
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)
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
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
419 xj(1:3*ndof, 1:3*ndof)= stiff_aa(1:3*ndof, 1:3*ndof)
423 tmp_d((i-1)*ndof+1:i*ndof) = ddu(1:ndof, i)
425 tmpforce(1:3*ndof) = matmul( stiff_ad(1:3*ndof, 1:nn*ndof), tmp_d(1:nn*ndof) )
427 tmp_a(1:3*ndof) = -matmul( xj(1:3*ndof, 1:3*ndof), tmpforce(1:3*ndof) )
429 ddaux(1:ndof, i) = tmp_a((i-1)*ndof+1:i*ndof)
435 totaldisp(:,1:nn) = u(:,:)+ du(:,:)
436 totaldisp(:,nn+1:nn+3) = aux(:,:) + ddaux(:,:)
443 stiff_ad(:, :) = 0.0d0
444 stiff_aa(:, :) = 0.0d0
449 elem(:,:) = (0.5d0*du(:,:)+u(:,:) ) +ecoord(:,:)
450 elem1(:,:) = (du(:,:)+u(:,:) ) +ecoord(:,:)
452 totaldisp(:,1:nn) = du(:,:)
453 if( iter == 1 ) aux(:,:) = 0.0d0
454 totaldisp(:,nn+1:nn+3) = aux(:,:)
458 if(
present(tt) .and. cdsys_id > 0 )
then
460 call fetch_tabledata(
mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
461 if( .not. ierr ) matlaniso = .true.
465 naturalcoord(:) = 0.0d0
466 call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
467 inverse(:, :)= inverse(:, :)*det
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
477 mtype = gausses(lx)%pMaterial%mtype
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"
496 if(
present(tt) .AND.
present(t0) )
then
498 ttc = dot_product(tt, spfunc)
499 tt0 = dot_product(t0, spfunc)
500 ttn = dot_product(tn, spfunc)
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
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(:)
523 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0;
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)
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) )
540 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
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)
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) )
560 if(
present(tt) .AND.
present(t0) )
then
561 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
563 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
570 b(1:6, 1:(nn+3)*ndof) = 0.0d0
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)
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
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)
612 b(:,j) = b(:,j)+b1(:,j)
617 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv(1:nn,1:3))
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
628 b(1:6, 1:(nn+3)*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)
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
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
655 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
657 = qf_a(1:3*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,nn*ndof+1:(nn+3)*ndof) )*wg
661 xj(1:9, 1:9)= stiff_aa(1:3*ndof, 1:3*ndof)
663 tmpdisp(:) = matmul( xj(:,:), qf_a(:) )
665 qf(i) = qf(i) - dot_product( stiff_ad(:,i), tmpdisp(:) )
673 (etype, nn, xx, yy, zz, tt, t0, &
674 gausses, vect, cdsys_id, coords)
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)
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)
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
712 if( cdsys_id > 0 )
then
714 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
715 if( .not. ierr ) matlaniso = .true.
723 naturalcoord(:) = 0.0d0
724 call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse)
725 inverse(:, :) = inverse(:, :)*det
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
738 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv(1:nn, 1:3) )
740 if( cdsys_id > 0 )
then
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"
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 )
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
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)
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
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
785 call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo(:), ierr, ina )
786 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
788 call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
789 if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
794 call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo0(:), ierr, ina )
795 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
797 call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
798 if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
810 call transformation(coordsys, tm)
811 epsth(:) = matmul( epsth(:), tm )
812 epsth(4:6) = epsth(4:6)*2.0d0
815 epsth(1:3) = thermal_eps
822 sgm(:) = matmul( d(:, :), epsth(:) )
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
833 xj(1:9,1:9)= stiff_aa(1:9, 1:9)
834 call calinverse(9, xj)
836 icdisp(1:9) = matmul( xj(:, :), vect_a(1:3*ndof) )
838 vect(1:nn*ndof) = vect(1:nn*ndof) - matmul( stiff_da(1:nn*ndof, 1:3*ndof), icdisp(1:3*ndof) )
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 getjacobian(fetype, nn, localcoord, elecoord, det, jacobian, inverse)
calculate Jacobian matrix, its determinant and inverse
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 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.
subroutine calinverse(nn, a)
calculate inverse of matrix a
This module summarizes all infomation of material properties.
integer(kind=kint), parameter d3
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.