8 use hecmw,
only : kint, kreal
13 real(kind=kreal),
parameter,
private :: i3(3,3) = reshape((/ &
14 & 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/))
21 real(kind=kreal),
intent(in) :: stress(6)
22 real(kind=kreal),
intent(out) :: mat(6, 6)
26 mat(1, 1) = 2.0d0*stress(1); mat(1, 2) = 0.0d0; mat(1, 3) = 0.0d0
27 mat(1, 4) = stress(4); mat(1, 5) = 0.0d0; mat(1, 6) = stress(6)
28 mat(2, 1) = mat(1, 2); mat(2, 2) = 2.0d0*stress(2); mat(2, 3) = 0.0d0
29 mat(2, 4) = stress(4); mat(2, 5) = stress(5); mat(2, 6) = 0.0d0
30 mat(3, 1) = mat(1, 3); mat(3, 2) = mat(2, 3); mat(3, 3) = 2.0d0*stress(3)
31 mat(3, 4) = 0.0d0; mat(3, 5) = stress(5); mat(3, 6) = stress(6)
33 mat(4, 1) = mat(1, 4); mat(4, 2) = mat(2, 4); mat(4, 3) = mat(3, 4)
34 mat(4, 4) = 0.5d0*( stress(1)+stress(2) ); mat(4, 5) = 0.5d0*stress(6); mat(4, 6) = 0.5d0*stress(5)
35 mat(5, 1) = mat(1, 5); mat(5, 2) = mat(2, 5); mat(5, 3) = mat(3, 5)
36 mat(5, 4) = mat(4, 5); mat(5, 5) = 0.5d0*( stress(3)+stress(2) ); mat(5, 6) = 0.5d0*stress(4)
37 mat(6, 1) = mat(1, 6); mat(6, 2) = mat(2, 6); mat(6, 3) = mat(3, 6)
38 mat(6, 4) = mat(4, 6); mat(6, 5) = mat(5, 6); mat(6, 6) = 0.5d0*( stress(1)+stress(3) );
51 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
52 time, tincr, u ,temperature)
61 integer(kind=kint),
intent(in) :: etype
62 integer(kind=kint),
intent(in) :: nn
63 real(kind=kreal),
intent(in) :: ecoord(3,nn)
65 real(kind=kreal),
intent(out) :: stiff(:,:)
66 integer(kind=kint),
intent(in) :: cdsys_id
67 real(kind=kreal),
intent(inout) :: coords(3,3)
68 real(kind=kreal),
intent(in) :: time
69 real(kind=kreal),
intent(in) :: tincr
70 real(kind=kreal),
intent(in) :: temperature(nn)
71 real(kind=kreal),
intent(in),
optional :: u(:,:)
75 integer(kind=kint) :: flag
76 integer(kind=kint),
parameter :: ndof = 3
77 real(kind=kreal) :: d(6, 6), b(6, ndof*nn), db(6, ndof*nn)
78 real(kind=kreal) :: gderiv(nn, 3), stress(6), mat(6, 6)
79 real(kind=kreal) :: det, wg
80 integer(kind=kint) :: i, j, lx, serr
81 real(kind=kreal) :: temp, naturalcoord(3)
82 real(kind=kreal) :: spfunc(nn), gdispderiv(3, 3)
83 real(kind=kreal) :: b1(6, ndof*nn), coordsys(3, 3)
84 real(kind=kreal) :: smat(9, 9), elem(3, nn)
85 real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
91 flag = gausses(1)%pMaterial%nlgeom_flag
92 if( .not.
present(u) ) flag = infinitesimal
93 elem(:, :) = ecoord(:, :)
94 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
101 if( cdsys_id > 0 )
then
103 if( serr == -1 ) stop
"Fail to setup local coordinate"
104 if( serr == -2 )
then
105 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
110 temp = dot_product(temperature, spfunc)
111 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
113 if( flag == updatelag )
then
114 call geomat_c3( gausses(lx)%stress, mat )
115 d(:, :) = d(:, :)-mat
119 b(1:6, 1:nn*ndof) = 0.0d0
121 b(1, 3*j-2)=gderiv(j, 1)
122 b(2, 3*j-1)=gderiv(j, 2)
123 b(3, 3*j )=gderiv(j, 3)
124 b(4, 3*j-2)=gderiv(j, 2)
125 b(4, 3*j-1)=gderiv(j, 1)
126 b(5, 3*j-1)=gderiv(j, 3)
127 b(5, 3*j )=gderiv(j, 2)
128 b(6, 3*j-2)=gderiv(j, 3)
129 b(6, 3*j )=gderiv(j, 1)
133 if( flag == totallag )
then
135 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
136 b1(1:6, 1:nn*ndof)=0.0d0
138 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
139 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
140 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
141 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
142 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
143 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
144 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
145 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
146 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
147 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
148 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
149 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
150 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
151 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
152 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
153 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
154 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
155 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
159 b(:, j) = b(:, j)+b1(:, j)
163 db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
166 stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
171 if( flag == totallag .OR. flag==updatelag )
then
172 stress(1:6) = gausses(lx)%stress
173 bn(1:9, 1:nn*ndof) = 0.0d0
175 bn(1, 3*j-2) = gderiv(j, 1)
176 bn(2, 3*j-1) = gderiv(j, 1)
177 bn(3, 3*j ) = gderiv(j, 1)
178 bn(4, 3*j-2) = gderiv(j, 2)
179 bn(5, 3*j-1) = gderiv(j, 2)
180 bn(6, 3*j ) = gderiv(j, 2)
181 bn(7, 3*j-2) = gderiv(j, 3)
182 bn(8, 3*j-1) = gderiv(j, 3)
183 bn(9, 3*j ) = gderiv(j, 3)
187 smat(j , j ) = stress(1)
188 smat(j , j+3) = stress(4)
189 smat(j , j+6) = stress(6)
190 smat(j+3, j ) = stress(4)
191 smat(j+3, j+3) = stress(2)
192 smat(j+3, j+6) = stress(5)
193 smat(j+6, j ) = stress(6)
194 smat(j+6, j+3) = stress(5)
195 smat(j+6, j+6) = stress(3)
197 sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
200 stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
213 subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
230 integer(kind=kint),
intent(in) :: etype, nn
231 real(kind=kreal),
intent(in) :: xx(:),yy(:),zz(:)
232 real(kind=kreal),
intent(in) :: params(0:6)
233 real(kind=kreal),
intent(inout) :: vect(:)
235 integer(kind=kint) ltype,nsize
238 integer(kind=kint) ndof
240 real(kind=kreal) h(nn)
241 real(kind=kreal) plx(nn), ply(nn), plz(nn)
242 real(kind=kreal) xj(3, 3), det, wg
243 integer(kind=kint) ivol, isuf
244 integer(kind=kint) nod(nn)
245 integer(kind=kint) IG2, LX, I , SURTYPE, NSUR
246 real(kind=kreal) vx, vy, vz, xcod, ycod, zcod
247 real(kind=kreal) ax, ay, az, rx, ry, rz, hx, hy, hz, val
248 real(kind=kreal) phx, phy, phz
249 real(kind=kreal) coefx, coefy, coefz
250 real(kind=kreal) normal(3), localcoord(3), elecoord(3, nn), deriv(nn, 3)
252 ax = 0.0d0; ay = 0.0d0; az = 0.0d0; rx = 0.0d0; ry = 0.0d0; rz = 0.0d0;
262 if( ltype.LT.10 )
then
264 if( ltype.EQ.5 )
then
272 else if( ltype.GE.10 )
then
274 call getsubface( etype, ltype/10, surtype, nod )
284 elecoord(1,i)=xx(nod(i))
285 elecoord(2,i)=yy(nod(i))
286 elecoord(3,i)=zz(nod(i))
290 call getshapefunc( surtype, localcoord(1:2), h(1:nsur) )
293 normal=
surfacenormal( surtype, nsur, localcoord(1:2), elecoord(:,1:nsur) )
295 vect(3*nod(i)-2)=vect(3*nod(i)-2)+val*wg*h(i)*normal(1)
296 vect(3*nod(i)-1)=vect(3*nod(i)-1)+val*wg*h(i)*normal(2)
297 vect(3*nod(i) )=vect(3*nod(i) )+val*wg*h(i)*normal(3)
312 xj(1,1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
313 xj(2,1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
314 xj(3,1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
316 det=xj(1,1)*xj(2,2)*xj(3,3) &
317 +xj(2,1)*xj(3,2)*xj(1,3) &
318 +xj(3,1)*xj(1,2)*xj(2,3) &
319 -xj(3,1)*xj(2,2)*xj(1,3) &
320 -xj(2,1)*xj(1,2)*xj(3,3) &
321 -xj(1,1)*xj(3,2)*xj(2,3)
328 xcod=dot_product( h(1:nn),xx(1:nn) )
329 ycod=dot_product( h(1:nn),yy(1:nn) )
330 zcod=dot_product( h(1:nn),zz(1:nn) )
331 hx=ax+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*rx
332 hy=ay+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*ry
333 hz=az+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*rz
337 coefx=rho*val*dabs(val)*phx
338 coefy=rho*val*dabs(val)*phy
339 coefz=rho*val*dabs(val)*phz
344 plx(i)=plx(i)+h(i)*wg*coefx
345 ply(i)=ply(i)+h(i)*wg*coefy
346 plz(i)=plz(i)+h(i)*wg*coefz
351 vect(3*i-2)=val*plx(i)
353 else if( ltype.EQ.2 )
then
355 vect(3*i-1)=val*ply(i)
357 else if( ltype.EQ.3 )
then
359 vect(3*i )=val*plz(i)
361 else if( ltype.EQ.4 )
then
365 vx=vx/sqrt(params(1)**2+params(2)**2+params(3)**2)
366 vy=vy/sqrt(params(1)**2+params(2)**2+params(3)**2)
367 vz=vz/sqrt(params(1)**2+params(2)**2+params(3)**2)
369 vect(3*i-2)=val*plx(i)*rho*vx
370 vect(3*i-1)=val*ply(i)*rho*vy
371 vect(3*i )=val*plz(i)*rho*vz
373 else if( ltype.EQ.5 )
then
387 (etype, nn, xx, yy, zz, tt, t0, gausses, &
388 vect, cdsys_id, coords)
398 integer(kind=kint),
parameter :: ndof = 3
399 integer(kind=kint),
intent(in) :: etype, nn
401 real(kind=kreal),
intent(in) :: xx(nn), yy(nn), zz(nn)
402 real(kind=kreal),
intent(in) :: tt(nn),t0(nn)
403 real(kind=kreal),
intent(out) :: vect(nn*ndof)
404 integer(kind=kint),
intent(in) :: cdsys_ID
405 real(kind=kreal),
intent(inout) :: coords(3, 3)
409 real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
410 real(kind=kreal) :: det, ecoord(3, nn)
411 integer(kind=kint) :: j, LX, serr
412 real(kind=kreal) :: epsth(6),sgm(6), h(nn), alpo(3), alpo0(3), coordsys(3, 3)
413 real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
414 real(kind=kreal) :: wg, outa(1), ina(1)
415 real(kind=kreal) :: tempc, temp0, thermal_eps, tm(6, 6)
416 logical :: ierr, matlaniso
422 if( cdsys_id > 0 )
then
424 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
425 if( .not. ierr ) matlaniso = .true.
439 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv )
441 if( cdsys_id > 0 )
then
442 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys, serr )
443 if( serr == -1 ) stop
"Fail to setup local coordinate"
444 if( serr == -2 )
then
445 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
451 b(1:6,1:nn*ndof)=0.0d0
453 b(1,3*j-2)=gderiv(j,1)
454 b(2,3*j-1)=gderiv(j,2)
455 b(3,3*j )=gderiv(j,3)
456 b(4,3*j-2)=gderiv(j,2)
457 b(4,3*j-1)=gderiv(j,1)
458 b(5,3*j-1)=gderiv(j,3)
459 b(5,3*j )=gderiv(j,2)
460 b(6,3*j-2)=gderiv(j,3)
461 b(6,3*j )=gderiv(j,1)
464 tempc = dot_product( h(1:nn), tt(1:nn) )
465 temp0 = dot_product( h(1:nn), t0(1:nn) )
467 call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
471 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
472 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
474 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
475 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
480 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
481 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
483 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
484 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
497 epsth(:) = matmul( epsth(:), tm )
498 epsth(4:6) = epsth(4:6)*2.0d0
501 epsth(1:3) = thermal_eps
508 sgm(:) = matmul( d(:, :), epsth(:) )
512 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:),b(:, :) )*wg
521 real(kind=kreal),
INTENT(IN) :: tt0
522 real(kind=kreal),
INTENT(IN) :: ttc
523 type( tmaterial ),
INTENT(IN) :: material
524 real(kind=kreal),
INTENT(IN) :: coordsys(3,3)
525 logical,
INTENT(IN) :: matlaniso
526 real(kind=kreal),
INTENT(OUT) :: epsth(6)
528 integer(kind=kint) :: j
529 real(kind=kreal) :: ina(1), outa(1)
531 real(kind=kreal) :: alp, alp0, alpo(3), alpo0(3), tm(6,6)
534 if( dabs(ttc-tt0) < 1.d-14 )
return
538 call fetch_tabledata( mc_orthoexp, material%dict, alpo(:), ierr, ina )
539 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
541 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
542 if( ierr ) outa(1) = material%variables(m_exapnsion)
547 call fetch_tabledata( mc_orthoexp, material%dict, alpo0(:), ierr, ina )
548 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
550 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
551 if( ierr ) outa(1) = material%variables(m_exapnsion)
559 epsth(:) = matmul( epsth(:), tm )
560 epsth(4:6) = epsth(4:6)*2.0d0
571 real(kind=kreal),
intent(in) :: rot(3,3)
572 real(kind=kreal),
intent(in) :: stress_in(3,3)
573 real(kind=kreal),
intent(out) :: stress_out(3,3)
575 real(kind=kreal) :: dr(3,3)
580 dr = matmul(dr,i3+0.5d0*rot)
582 stress_out = matmul(dr,stress_in)
583 stress_out = matmul(stress_out,transpose(dr))
586 subroutine update_stress3d( flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn, hdflag )
593 integer(kind=kint),
intent(in) :: flag
594 real(kind=kreal),
intent(in) :: rot(3,3)
595 real(kind=kreal),
intent(in) :: dstrain(6)
596 real(kind=kreal),
intent(in) :: f(3,3)
597 real(kind=kreal),
intent(in) :: coordsys(3,3)
598 real(kind=kreal),
intent(in) :: time
599 real(kind=kreal),
intent(in) :: tincr
600 real(kind=kreal),
intent(in) :: ttc
601 real(kind=kreal),
intent(in) :: tt0
602 real(kind=kreal),
intent(in) :: ttn
603 integer(kind=kint),
intent(in),
optional :: hdflag
605 integer(kind=kint) :: mtype, i, j, k
606 integer(kind=kint) :: isep, hdflag_in
607 real(kind=kreal) :: d(6,6), dstress(6), dumstress(3,3), dum(3,3), trd, det
608 real(kind=kreal) :: tensor(6)
609 real(kind=kreal) :: eigval(3)
610 real(kind=kreal) :: princ(3,3), norm
612 mtype = gauss%pMaterial%mtype
614 if( iselastoplastic(mtype) .OR. mtype == norton )
then
621 if(
present(hdflag) ) hdflag_in = hdflag
623 call matlmatrix( gauss, d3, d, time, tincr, coordsys, ttc, isep, hdflag=hdflag_in )
625 if( flag == infinitesimal )
then
627 dstress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6)-gauss%strain_bak(1:6) )
628 gauss%stress(1:6) = gauss%stress_bak(1:6) + dstress(1:6)
629 if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 )
then
630 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
631 gauss%stress = real(gauss%stress)
633 gauss%strain_energy = gauss%strain_energy_bak+dot_product(gauss%stress_bak(1:6),dstrain(1:6)-gauss%strain_bak(1:6))
634 gauss%strain_energy = gauss%strain_energy+0.5d0*dot_product(dstress,dstrain(1:6)-gauss%strain_bak(1:6))
637 else if( flag == totallag )
then
639 if( ishyperelastic(mtype) .OR. mtype == userelastic .OR. mtype == usermaterial )
then
640 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, temp=0.d0, tempn=0.d0, hdflag=hdflag_in )
641 else if( ( isviscoelastic(mtype) .OR. mtype == norton ) .AND. tincr /= 0.0d0 )
then
642 gauss%pMaterial%mtype=mtype
643 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
645 gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
646 gauss%strain_energy = 0.5d0*dot_product(gauss%stress(1:6),dstrain(1:6))
649 else if( flag == updatelag )
then
653 if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 )
then
654 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, tt0, hdflag=hdflag_in )
656 dstress = real( matmul( d(1:6,1:6), dstrain(1:6) ) )
657 dumstress(1,1) = gauss%stress_bak(1)
658 dumstress(2,2) = gauss%stress_bak(2)
659 dumstress(3,3) = gauss%stress_bak(3)
660 dumstress(1,2) = gauss%stress_bak(4); dumstress(2,1)=dumstress(1,2)
661 dumstress(2,3) = gauss%stress_bak(5); dumstress(3,2)=dumstress(2,3)
662 dumstress(3,1) = gauss%stress_bak(6); dumstress(1,3)=dumstress(3,1)
665 trd = dstrain(1)+dstrain(2)+dstrain(3)
666 dum(:,:) = dumstress + matmul( rot,dumstress ) - matmul( dumstress, rot ) - dumstress*trd
669 gauss%stress(1) = dum(1,1) + dstress(1)
670 gauss%stress(2) = dum(2,2) + dstress(2)
671 gauss%stress(3) = dum(3,3) + dstress(3)
672 gauss%stress(4) = dum(1,2) + dstress(4)
673 gauss%stress(5) = dum(2,3) + dstress(5)
674 gauss%stress(6) = dum(3,1) + dstress(6)
676 if( mtype == usermaterial )
then
677 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, temp=0.d0, tempn=0.d0, hdflag=hdflag_in )
678 else if( mtype == norton )
then
680 if( tincr /= 0.0d0 .AND. any( gauss%stress /= 0.0d0 ) )
then
682 call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
686 gauss%strain_energy = gauss%strain_energy_bak+dot_product(gauss%stress(1:6),dstrain(1:6))
687 gauss%strain_energy = gauss%strain_energy-0.5d0*dot_product(dstress(1:6),dstrain(1:6))
691 if( iselastoplastic(mtype) )
then
692 call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
693 gauss%istatus(1), gauss%fstatus, gauss%plpotential, ttc )
694 gauss%strain_energy = gauss%strain_energy+gauss%plpotential
700 if( flag == infinitesimal )
then
701 gauss%stress_out(1:6) = gauss%stress(1:6)
702 gauss%strain_out(1:6) = gauss%strain(1:6)
705 if( flag == totallag )
then
706 dumstress(1,1) = gauss%stress(1)
707 dumstress(2,2) = gauss%stress(2)
708 dumstress(3,3) = gauss%stress(3)
709 dumstress(1,2) = gauss%stress(4); dumstress(2,1)=dumstress(1,2)
710 dumstress(2,3) = gauss%stress(5); dumstress(3,2)=dumstress(2,3)
711 dumstress(3,1) = gauss%stress(6); dumstress(1,3)=dumstress(3,1)
714 if( det == 0.d0 ) stop
"Fail to convert stress: detF=0"
716 dumstress(1:3,1:3) = matmul(dumstress(1:3,1:3),transpose(f(1:3,1:3)))
717 dumstress(1:3,1:3) = (1.d0/det)*matmul(f(1:3,1:3),dumstress(1:3,1:3))
719 gauss%stress_out(1) = dumstress(1,1)
720 gauss%stress_out(2) = dumstress(2,2)
721 gauss%stress_out(3) = dumstress(3,3)
722 gauss%stress_out(4) = dumstress(1,2)
723 gauss%stress_out(5) = dumstress(2,3)
724 gauss%stress_out(6) = dumstress(3,1)
725 else if( flag == updatelag )
then
726 gauss%stress_out(1:6) = gauss%stress(1:6)
730 dum(1:3,1:3) = matmul(f(1:3,1:3),transpose(f(1:3,1:3)))
740 if( eigval(k) <= 0.d0 ) stop
"Fail to calc log strain: stretch<0"
741 eigval(k) = 0.5d0*dlog(eigval(k))
742 norm = dsqrt(dot_product(princ(1:3,k),princ(1:3,k)))
743 if( norm <= 0.d0 ) stop
"Fail to calc log strain: stretch direction vector=0"
744 princ(1:3,k) = princ(1:3,k)/norm
750 dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
754 gauss%strain_out(1) = dum(1,1)
755 gauss%strain_out(2) = dum(2,2)
756 gauss%strain_out(3) = dum(3,3)
757 gauss%strain_out(4) = 2.d0*dum(1,2)
758 gauss%strain_out(5) = 2.d0*dum(2,3)
759 gauss%strain_out(6) = 2.d0*dum(3,1)
763 gauss%stress_out(1:6) = gauss%stress(1:6)
764 gauss%strain_out(1:6) = gauss%strain(1:6)
772 (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
773 gausses, iter, time, tincr, tt, t0, tn)
783 integer(kind=kint),
intent(in) :: etype
784 integer(kind=kint),
intent(in) :: nn
785 real(kind=kreal),
intent(in) :: ecoord(3, nn)
786 real(kind=kreal),
intent(in) :: u(3, nn)
787 real(kind=kreal),
intent(in) :: ddu(3, nn)
788 integer(kind=kint),
intent(in) :: cdsys_id
789 real(kind=kreal),
intent(inout) :: coords(3, 3)
790 real(kind=kreal),
intent(out) :: qf(nn*3)
792 integer,
intent(in) :: iter
793 real(kind=kreal),
intent(in) :: time
794 real(kind=kreal),
intent(in) :: tincr
795 real(kind=kreal),
intent(in) :: tt(nn)
796 real(kind=kreal),
intent(in) :: t0(nn)
797 real(kind=kreal),
intent(in) :: tn(nn)
800 integer(kind=kint) :: flag
801 integer(kind=kint),
parameter :: ndof = 3
802 real(kind=kreal) :: b(6,ndof*nn), b1(6,ndof*nn), spfunc(nn), ina(1)
803 real(kind=kreal) :: gderiv(nn,3), gderiv1(nn,3), gdispderiv(3,3), f(3,3), det, det1, wg, ttc, tt0, ttn
804 integer(kind=kint) :: j, lx, serr
805 real(kind=kreal) :: naturalcoord(3), rot(3,3), epsth(6)
806 real(kind=kreal) :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3)
807 real(kind=kreal) :: dstrain(6)
808 real(kind=kreal) :: alpo(3)
809 logical :: ierr, matlaniso
813 flag = gausses(1)%pMaterial%nlgeom_flag
814 elem(:,:) = ecoord(:,:)
815 totaldisp(:,:) = u(:,:)+ ddu(:,:)
817 elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
818 elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
820 totaldisp(:,:) = ddu(:,:)
825 call fetch_tabledata(
mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
826 if( .not. ierr ) matlaniso = .true.
833 if( cdsys_id > 0 )
then
834 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
835 if( serr == -1 ) stop
"Fail to setup local coordinate"
836 if( serr == -2 )
then
837 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
848 ttc = dot_product(tt, spfunc)
849 tt0 = dot_product(t0, spfunc)
850 ttn = dot_product(tn, spfunc)
855 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
856 dstrain(1) = gdispderiv(1, 1)
857 dstrain(2) = gdispderiv(2, 2)
858 dstrain(3) = gdispderiv(3, 3)
859 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
860 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
861 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
862 dstrain(:) = dstrain(:)-epsth(:)
864 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0;
866 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
870 dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
871 dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
872 dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
873 dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
874 +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
875 dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
876 +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
877 dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
878 +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
880 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
881 f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
885 rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
886 rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
887 rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
889 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
891 call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
892 f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn)+ddu(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
897 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
903 b(1:6, 1:nn*ndof) = 0.0d0
905 b(1,3*j-2) = gderiv(j, 1)
906 b(2,3*j-1) = gderiv(j, 2)
907 b(3,3*j ) = gderiv(j, 3)
908 b(4,3*j-2) = gderiv(j, 2)
909 b(4,3*j-1) = gderiv(j, 1)
910 b(5,3*j-1) = gderiv(j, 3)
911 b(5,3*j ) = gderiv(j, 2)
912 b(6,3*j-2) = gderiv(j, 3)
913 b(6,3*j ) = gderiv(j, 1)
921 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
922 b1(1:6, 1:nn*ndof)=0.0d0
924 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
925 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
926 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
927 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
928 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
929 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
930 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
931 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
932 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
933 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
934 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
935 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
936 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
937 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
938 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
939 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
940 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
941 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
945 b(:,j) = b(:,j)+b1(:,j)
951 b(1:6, 1:nn*ndof) = 0.0d0
953 b(1, 3*j-2) = gderiv(j, 1)
954 b(2, 3*j-1) = gderiv(j, 2)
955 b(3, 3*j ) = gderiv(j, 3)
956 b(4, 3*j-2) = gderiv(j, 2)
957 b(4, 3*j-1) = gderiv(j, 1)
958 b(5, 3*j-1) = gderiv(j, 3)
959 b(5, 3*j ) = gderiv(j, 2)
960 b(6, 3*j-2) = gderiv(j, 3)
961 b(6, 3*j ) = gderiv(j, 1)
969 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
972 gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
987 integer(kind=kint),
intent(in) :: etype, nn
989 real(kind=kreal),
intent(out) :: ndstrain(nn,6)
990 real(kind=kreal),
intent(out) :: ndstress(nn,6)
995 real(kind=kreal) :: temp(12)
1004 temp(1:6) = temp(1:6) +gausses(i)%strain_out(1:6)
1005 temp(7:12) = temp(7:12)+gausses(i)%stress_out(1:6)
1008 temp(1:12) = temp(1:12)/ic
1011 ndstrain(i, 1:6) = temp(1:6)
1012 ndstress(i, 1:6) = temp(7:12)
1028 integer(kind=kint),
intent(in) :: etype
1030 real(kind=kreal),
intent(out) :: strain(6)
1031 real(kind=kreal),
intent(out) :: stress(6)
1039 strain(:) = 0.0d0; stress(:) = 0.0d0
1044 strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1045 stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1048 strain(:) = strain(:)/ic
1049 stress(:) = stress(:)/ic
1059 integer(kind=kint),
intent(in) :: etype, nn
1060 real(kind=kreal),
intent(in) :: xx(:), yy(:), zz(:)
1064 real(kind=kreal) :: xj(3, 3), det
1065 integer(kind=kint) :: lx
1066 real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1079 xj(1, 1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
1080 xj(2, 1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
1081 xj(3, 1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
1084 det = xj(1, 1)*xj(2, 2)*xj(3, 3) &
1085 +xj(2, 1)*xj(3, 2)*xj(1, 3) &
1086 +xj(3, 1)*xj(1, 2)*xj(2, 3) &
1087 -xj(3, 1)*xj(2, 2)*xj(1, 3) &
1088 -xj(2, 1)*xj(1, 2)*xj(3, 3) &
1089 -xj(1, 1)*xj(3, 2)*xj(2, 3)
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
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.
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural 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.
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
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 defines common data and basic structures for analysis.
integer(kind=kint), parameter kopss_solution
integer(kind=kint) opsstype
real(kind=kreal), pointer ref_temp
REFTEMP.
This module manages calculation relates with materials.
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
subroutine stressupdate(gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn, hdflag)
Update strain and stress for elastic and hyperelastic materials.
This module provides common functions of Solid elements.
subroutine elementstress_c3(etype, gausses, strain, stress)
subroutine tload_c3(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
Distributed external load.
subroutine hughes_winget_rotation_3d(rot, stress_in, stress_out)
subroutine geomat_c3(stress, mat)
subroutine update_c3(etype, nn, ecoord, u, ddu, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update strain and stress inside element.
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, hdflag)
subroutine stf_c3(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix of general solid elements.
real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
Volume of element.
subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
This module provides aux functions.
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
subroutine get_principal(tensor, eigval, princmatrix)
subroutine transformation(jacob, tm)
subroutine calinverse(NN, A)
calculate inverse of matrix a
This module summarizes all information 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.