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),
optional :: 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"
109 if(
present(temperature) )
then
111 temp = dot_product(temperature, spfunc)
112 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
114 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
117 if( flag == updatelag )
then
118 call geomat_c3( gausses(lx)%stress, mat )
119 d(:, :) = d(:, :)-mat
123 b(1:6, 1:nn*ndof) = 0.0d0
125 b(1, 3*j-2)=gderiv(j, 1)
126 b(2, 3*j-1)=gderiv(j, 2)
127 b(3, 3*j )=gderiv(j, 3)
128 b(4, 3*j-2)=gderiv(j, 2)
129 b(4, 3*j-1)=gderiv(j, 1)
130 b(5, 3*j-1)=gderiv(j, 3)
131 b(5, 3*j )=gderiv(j, 2)
132 b(6, 3*j-2)=gderiv(j, 3)
133 b(6, 3*j )=gderiv(j, 1)
137 if( flag == totallag )
then
139 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
140 b1(1:6, 1:nn*ndof)=0.0d0
142 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
143 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
144 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
145 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
146 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
147 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
148 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
149 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
150 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
151 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
152 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
153 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
154 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
155 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
156 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
157 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
158 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
159 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
163 b(:, j) = b(:, j)+b1(:, j)
167 db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
168 forall( i=1:nn*ndof, j=1:nn*ndof )
169 stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
173 if( flag == totallag .OR. flag==updatelag )
then
174 stress(1:6) = gausses(lx)%stress
175 bn(1:9, 1:nn*ndof) = 0.0d0
177 bn(1, 3*j-2) = gderiv(j, 1)
178 bn(2, 3*j-1) = gderiv(j, 1)
179 bn(3, 3*j ) = gderiv(j, 1)
180 bn(4, 3*j-2) = gderiv(j, 2)
181 bn(5, 3*j-1) = gderiv(j, 2)
182 bn(6, 3*j ) = gderiv(j, 2)
183 bn(7, 3*j-2) = gderiv(j, 3)
184 bn(8, 3*j-1) = gderiv(j, 3)
185 bn(9, 3*j ) = gderiv(j, 3)
189 smat(j , j ) = stress(1)
190 smat(j , j+3) = stress(4)
191 smat(j , j+6) = stress(6)
192 smat(j+3, j ) = stress(4)
193 smat(j+3, j+3) = stress(2)
194 smat(j+3, j+6) = stress(5)
195 smat(j+6, j ) = stress(6)
196 smat(j+6, j+3) = stress(5)
197 smat(j+6, j+6) = stress(3)
199 sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
200 forall( i=1:nn*ndof, j=1:nn*ndof )
201 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*val*phx
338 coefy=rho*val*val*phy
339 coefz=rho*val*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 )
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)
535 call fetch_tabledata( mc_orthoexp, material%dict, alpo(:), ierr, ina )
536 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
538 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
539 if( ierr ) outa(1) = material%variables(m_exapnsion)
544 call fetch_tabledata( mc_orthoexp, material%dict, alpo0(:), ierr, ina )
545 if( ierr ) stop
"Fails in fetching orthotropic expansion coefficient!"
547 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
548 if( ierr ) outa(1) = material%variables(m_exapnsion)
556 epsth(:) = matmul( epsth(:), tm )
557 epsth(4:6) = epsth(4:6)*2.0d0
568 real(kind=kreal),
intent(in) :: rot(3,3)
569 real(kind=kreal),
intent(in) :: stress_in(3,3)
570 real(kind=kreal),
intent(out) :: stress_out(3,3)
572 real(kind=kreal) :: dr(3,3)
577 dr = matmul(dr,i3+0.5d0*rot)
579 stress_out = matmul(dr,stress_in)
580 stress_out = matmul(stress_out,transpose(dr))
583 subroutine update_stress3d( flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn )
590 integer(kind=kint),
intent(in) :: flag
591 real(kind=kreal),
intent(in) :: rot(3,3)
592 real(kind=kreal),
intent(in) :: dstrain(6)
593 real(kind=kreal),
intent(in) :: f(3,3)
594 real(kind=kreal),
intent(in) :: coordsys(3,3)
595 real(kind=kreal),
intent(in) :: time
596 real(kind=kreal),
intent(in) :: tincr
597 real(kind=kreal),
intent(in),
optional :: ttc
598 real(kind=kreal),
intent(in),
optional :: tt0
599 real(kind=kreal),
intent(in),
optional :: ttn
601 integer(kind=kint) :: mtype, i, j, k
602 integer(kind=kint) :: isep
603 real(kind=kreal) :: d(6,6), dstress(6), dumstress(3,3), dum(3,3), trd, det
604 real(kind=kreal) :: tensor(6)
605 real(kind=kreal) :: eigval(3)
606 real(kind=kreal) :: princ(3,3), norm
608 mtype = gauss%pMaterial%mtype
610 if( iselastoplastic(mtype) .OR. mtype == norton )
then
616 if(
present(ttc) .AND.
present(ttn) )
then
617 call matlmatrix( gauss, d3, d, time, tincr, coordsys, ttc, isep )
619 call matlmatrix( gauss, d3, d, time, tincr, coordsys, isep=isep)
622 if( flag == infinitesimal )
then
624 gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
625 if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 )
then
626 if(
present(ttc) .AND.
present(ttn) )
then
627 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn )
629 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
631 gauss%stress = real(gauss%stress)
634 else if( flag == totallag )
then
636 if( ishyperelastic(mtype) .OR. mtype == userelastic .OR. mtype == usermaterial )
then
637 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys )
638 else if( ( isviscoelastic(mtype) .OR. mtype == norton ) .AND. tincr /= 0.0d0 )
then
639 gauss%pMaterial%mtype=mtype
640 if(
present(ttc) .AND.
present(ttn) )
then
641 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn )
643 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
646 gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
649 else if( flag == updatelag )
then
653 if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 )
then
654 if(
present(ttc) .AND.
present(ttn) )
then
655 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, tt0 )
657 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
660 dstress = real( matmul( d(1:6,1:6), dstrain(1:6) ) )
661 dumstress(1,1) = gauss%stress_bak(1)
662 dumstress(2,2) = gauss%stress_bak(2)
663 dumstress(3,3) = gauss%stress_bak(3)
664 dumstress(1,2) = gauss%stress_bak(4); dumstress(2,1)=dumstress(1,2)
665 dumstress(2,3) = gauss%stress_bak(5); dumstress(3,2)=dumstress(2,3)
666 dumstress(3,1) = gauss%stress_bak(6); dumstress(1,3)=dumstress(3,1)
669 trd = dstrain(1)+dstrain(2)+dstrain(3)
670 dum(:,:) = dumstress + matmul( rot,dumstress ) - matmul( dumstress, rot ) + dumstress*trd
673 gauss%stress(1) = dum(1,1) + dstress(1)
674 gauss%stress(2) = dum(2,2) + dstress(2)
675 gauss%stress(3) = dum(3,3) + dstress(3)
676 gauss%stress(4) = dum(1,2) + dstress(4)
677 gauss%stress(5) = dum(2,3) + dstress(5)
678 gauss%stress(6) = dum(3,1) + dstress(6)
680 if( mtype == usermaterial )
then
681 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys )
682 else if( mtype == norton )
then
684 if( tincr /= 0.0d0 .AND. any( gauss%stress /= 0.0d0 ) )
then
686 if(
present(ttc) .AND.
present(ttn) )
then
687 call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr, ttc, ttn )
689 call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr )
697 if( iselastoplastic(mtype) )
then
698 if(
present(ttc) )
then
699 call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
700 gauss%istatus(1), gauss%fstatus, ttc )
702 call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
703 gauss%istatus(1), gauss%fstatus )
710 if( flag == infinitesimal )
then
711 gauss%stress_out(1:6) = gauss%stress(1:6)
712 gauss%strain_out(1:6) = gauss%strain(1:6)
715 if( flag == totallag )
then
716 dumstress(1,1) = gauss%stress(1)
717 dumstress(2,2) = gauss%stress(2)
718 dumstress(3,3) = gauss%stress(3)
719 dumstress(1,2) = gauss%stress(4); dumstress(2,1)=dumstress(1,2)
720 dumstress(2,3) = gauss%stress(5); dumstress(3,2)=dumstress(2,3)
721 dumstress(3,1) = gauss%stress(6); dumstress(1,3)=dumstress(3,1)
724 if( det == 0.d0 ) stop
"Fail to convert stress: detF=0"
726 dumstress(1:3,1:3) = matmul(dumstress(1:3,1:3),transpose(f(1:3,1:3)))
727 dumstress(1:3,1:3) = (1.d0/det)*matmul(f(1:3,1:3),dumstress(1:3,1:3))
729 gauss%stress_out(1) = dumstress(1,1)
730 gauss%stress_out(2) = dumstress(2,2)
731 gauss%stress_out(3) = dumstress(3,3)
732 gauss%stress_out(4) = dumstress(1,2)
733 gauss%stress_out(5) = dumstress(2,3)
734 gauss%stress_out(6) = dumstress(3,1)
735 else if( flag == updatelag )
then
736 gauss%stress_out(1:6) = gauss%stress(1:6)
740 dum(1:3,1:3) = matmul(f(1:3,1:3),transpose(f(1:3,1:3)))
750 if( eigval(k) <= 0.d0 ) stop
"Fail to calc log strain: stretch<0"
751 eigval(k) = 0.5d0*dlog(eigval(k))
752 norm = dsqrt(dot_product(princ(1:3,k),princ(1:3,k)))
753 if( norm <= 0.d0 ) stop
"Fail to calc log strain: stretch direction vector=0"
754 princ(1:3,k) = princ(1:3,k)/norm
760 dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
764 gauss%strain_out(1) = dum(1,1)
765 gauss%strain_out(2) = dum(2,2)
766 gauss%strain_out(3) = dum(3,3)
767 gauss%strain_out(4) = 2.d0*dum(1,2)
768 gauss%strain_out(5) = 2.d0*dum(2,3)
769 gauss%strain_out(6) = 2.d0*dum(3,1)
773 gauss%stress_out(1:6) = gauss%stress(1:6)
774 gauss%strain_out(1:6) = gauss%strain(1:6)
782 (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
783 gausses, iter, time, tincr, tt, t0, tn)
793 integer(kind=kint),
intent(in) :: etype
794 integer(kind=kint),
intent(in) :: nn
795 real(kind=kreal),
intent(in) :: ecoord(3, nn)
796 real(kind=kreal),
intent(in) :: u(3, nn)
797 real(kind=kreal),
intent(in) :: ddu(3, nn)
798 integer(kind=kint),
intent(in) :: cdsys_id
799 real(kind=kreal),
intent(inout) :: coords(3, 3)
800 real(kind=kreal),
intent(out) :: qf(nn*3)
802 integer,
intent(in) :: iter
803 real(kind=kreal),
intent(in) :: time
804 real(kind=kreal),
intent(in) :: tincr
805 real(kind=kreal),
intent(in),
optional :: tt(nn)
806 real(kind=kreal),
intent(in),
optional :: t0(nn)
807 real(kind=kreal),
intent(in),
optional :: tn(nn)
810 integer(kind=kint) :: flag
811 integer(kind=kint),
parameter :: ndof = 3
812 real(kind=kreal) :: b(6,ndof*nn), b1(6,ndof*nn), spfunc(nn), ina(1)
813 real(kind=kreal) :: gderiv(nn,3), gderiv1(nn,3), gdispderiv(3,3), f(3,3), det, det1, wg, ttc,tt0, ttn
814 integer(kind=kint) :: j, lx, serr
815 real(kind=kreal) :: naturalcoord(3), rot(3,3), epsth(6)
816 real(kind=kreal) :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3)
817 real(kind=kreal) :: dstrain(6)
818 real(kind=kreal) :: alpo(3)
819 logical :: ierr, matlaniso
823 flag = gausses(1)%pMaterial%nlgeom_flag
824 elem(:,:) = ecoord(:,:)
825 totaldisp(:,:) = u(:,:)+ ddu(:,:)
827 elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
828 elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
830 totaldisp(:,:) = ddu(:,:)
834 if(
present(tt) .and. cdsys_id > 0 )
then
836 call fetch_tabledata(
mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
837 if( .not. ierr ) matlaniso = .true.
845 if( cdsys_id > 0 )
then
846 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
847 if( serr == -1 ) stop
"Fail to setup local coordinate"
848 if( serr == -2 )
then
849 write(*, *)
"WARNING! Cannot setup local coordinate, it is modified automatically"
859 if(
present(tt) .AND.
present(t0) )
then
861 ttc = dot_product(tt, spfunc)
862 tt0 = dot_product(t0, spfunc)
863 ttn = dot_product(tn, spfunc)
869 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
870 dstrain(1) = gdispderiv(1, 1)
871 dstrain(2) = gdispderiv(2, 2)
872 dstrain(3) = gdispderiv(3, 3)
873 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
874 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
875 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
876 dstrain(:) = dstrain(:)-epsth(:)
878 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0;
880 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
884 dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
885 dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
886 dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
887 dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
888 +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
889 dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
890 +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
891 dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
892 +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
894 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
895 f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
899 rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
900 rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
901 rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
903 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
905 call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
906 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) )
911 if(
present(tt) .AND.
present(t0) )
then
912 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
914 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
921 b(1:6, 1:nn*ndof) = 0.0d0
923 b(1,3*j-2) = gderiv(j, 1)
924 b(2,3*j-1) = gderiv(j, 2)
925 b(3,3*j ) = gderiv(j, 3)
926 b(4,3*j-2) = gderiv(j, 2)
927 b(4,3*j-1) = gderiv(j, 1)
928 b(5,3*j-1) = gderiv(j, 3)
929 b(5,3*j ) = gderiv(j, 2)
930 b(6,3*j-2) = gderiv(j, 3)
931 b(6,3*j ) = gderiv(j, 1)
939 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
940 b1(1:6, 1:nn*ndof)=0.0d0
942 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
943 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
944 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
945 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
946 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
947 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
948 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
949 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
950 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
951 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
952 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
953 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
954 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
955 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
956 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
957 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
958 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
959 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
963 b(:,j) = b(:,j)+b1(:,j)
969 b(1:6, 1:nn*ndof) = 0.0d0
971 b(1, 3*j-2) = gderiv(j, 1)
972 b(2, 3*j-1) = gderiv(j, 2)
973 b(3, 3*j ) = gderiv(j, 3)
974 b(4, 3*j-2) = gderiv(j, 2)
975 b(4, 3*j-1) = gderiv(j, 1)
976 b(5, 3*j-1) = gderiv(j, 3)
977 b(5, 3*j ) = gderiv(j, 2)
978 b(6, 3*j-2) = gderiv(j, 3)
979 b(6, 3*j ) = gderiv(j, 1)
987 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
1003 integer(kind=kint),
intent(in) :: etype, nn
1005 real(kind=kreal),
intent(out) :: ndstrain(nn,6)
1006 real(kind=kreal),
intent(out) :: ndstress(nn,6)
1011 real(kind=kreal) :: temp(12)
1020 temp(1:6) = temp(1:6) +gausses(i)%strain_out(1:6)
1021 temp(7:12) = temp(7:12)+gausses(i)%stress_out(1:6)
1024 temp(1:12) = temp(1:12)/ic
1027 ndstrain(i, 1:6) = temp(1:6)
1028 ndstress(i, 1:6) = temp(7:12)
1044 integer(kind=kint),
intent(in) :: etype
1046 real(kind=kreal),
intent(out) :: strain(6)
1047 real(kind=kreal),
intent(out) :: stress(6)
1055 strain(:) = 0.0d0; stress(:) = 0.0d0
1060 strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1061 stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1064 strain(:) = strain(:)/ic
1065 stress(:) = stress(:)/ic
1075 integer(kind=kint),
intent(in) :: etype, nn
1076 real(kind=kreal),
intent(in) :: xx(:), yy(:), zz(:)
1080 real(kind=kreal) :: xj(3, 3), det
1081 integer(kind=kint) :: lx
1082 real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1095 xj(1, 1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
1096 xj(2, 1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
1097 xj(3, 1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
1100 det = xj(1, 1)*xj(2, 2)*xj(3, 3) &
1101 +xj(2, 1)*xj(3, 2)*xj(1, 3) &
1102 +xj(3, 1)*xj(1, 2)*xj(2, 3) &
1103 -xj(3, 1)*xj(2, 2)*xj(1, 3) &
1104 -xj(2, 1)*xj(1, 2)*xj(3, 3) &
1105 -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 coodinate 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 deivatives of shape fucntion in natural coordiante 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 defined coomon 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 stressupdate(gauss, secttype, strain, stress, cdsys, time, dtime, temp, tempn)
Update strain and stress for elastic and hyperelastic materials.
subroutine matlmatrix(gauss, secttype, matrix, time, dtime, cdsys, temperature, isep)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine elementstress_c3(etype, gausses, strain, stress)
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.
real(kind=kreal) function volume_c3(etype, nn, xx, yy, zz)
Volume of element.
subroutine hughes_winget_rotation_3d(rot, stress_in, stress_out)
subroutine update_stress3d(flag, gauss, rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn)
subroutine geomat_c3(stress, mat)
subroutine stf_c3(etype, nn, ecoord, gausses, stiff, cdsys_id, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix of general solid elements.
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)
Distrubuted external load.
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, epsth)
subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
This module provides aux functions.
subroutine calinverse(nn, a)
calculate inverse of matrix a
real(kind=kreal) function determinant33(xj)
Compute determinant for 3*3 matrix.
subroutine get_principal(tensor, eigval, princmatrix)
subroutine transformation(jacob, tm)
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.