FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_3d.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
7
8 use hecmw, only : kint, kreal
10
11 implicit none
12
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/))
15
16contains
17 !----------------------------------------------------------------------*
18 subroutine geomat_c3(stress, mat)
19 !----------------------------------------------------------------------*
20
21 real(kind=kreal), intent(in) :: stress(6)
22 real(kind=kreal), intent(out) :: mat(6, 6)
23
24 !---------------------------------------------------------------------
25
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)
32
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) );
39
40 end subroutine
41
42
43 !=====================================================================*
45 !
49 !----------------------------------------------------------------------*
50 subroutine stf_c3 &
51 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
52 time, tincr, u ,temperature)
53 !----------------------------------------------------------------------*
54
55 use mmechgauss
56 use m_matmatrix
58
59 !---------------------------------------------------------------------
60
61 integer(kind=kint), intent(in) :: etype
62 integer(kind=kint), intent(in) :: nn
63 real(kind=kreal), intent(in) :: ecoord(3,nn)
64 type(tgaussstatus), intent(in) :: gausses(:)
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(:,:)
72
73 !---------------------------------------------------------------------
74
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)
86
87 !---------------------------------------------------------------------
88
89 stiff(:, :) = 0.0d0
90 ! we suppose the same material type in the element
91 flag = gausses(1)%pMaterial%nlgeom_flag
92 if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
93 elem(:, :) = ecoord(:, :)
94 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
95
96 do lx = 1, numofquadpoints(etype)
97
98 call getquadpoint( etype, lx, naturalcoord(:) )
99 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
100
101 if( cdsys_id > 0 ) then
102 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
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"
106 end if
107 end if
108
109 if( present(temperature) ) then
110 call getshapefunc(etype, naturalcoord, spfunc)
111 temp = dot_product(temperature, spfunc)
112 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
113 else
114 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
115 end if
116
117 if( flag == updatelag ) then
118 call geomat_c3( gausses(lx)%stress, mat )
119 d(:, :) = d(:, :)-mat
120 endif
121
122 wg = getweight(etype, lx)*det
123 b(1:6, 1:nn*ndof) = 0.0d0
124 do j = 1, nn
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)
134 enddo
135
136 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
137 if( flag == totallag ) then
138 ! ---dudx(i, j) ==> gdispderiv(i, j)
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
141 do j=1, nn
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)
160 end do
161 ! ---BL = BL0 + BL1
162 do j=1, nn*ndof
163 b(:, j) = b(:, j)+b1(:, j)
164 end do
165 end if
166
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
170 end forall
171
172 ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
173 if( flag == totallag .OR. flag==updatelag ) then
174 stress(1:6) = gausses(lx)%stress
175 bn(1:9, 1:nn*ndof) = 0.0d0
176 do j = 1, nn
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)
186 end do
187 smat(:, :) = 0.0d0
188 do j = 1, 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)
198 end do
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
202 end forall
203
204 end if
205
206 enddo ! gauss roop
207
208 end subroutine stf_c3
209
210
212 !----------------------------------------------------------------------*
213 subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
214 !----------------------------------------------------------------------*
215 !**
216 !** SET DLOAD
217 !**
218 ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
219 ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
220 ! BZ LTYPE=3 :BODY FORCE IN Z-DIRECTION
221 ! GRAV LTYPE=4 :GRAVITY FORCE
222 ! CENT LTYPE=5 :CENTRIFUGAL LOAD
223 ! P1 LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
224 ! P2 LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
225 ! P3 LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
226 ! P4 LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
227 ! P5 LTYPE=50 :TRACTION IN NORMAL-DIRECTION FOR FACE-5
228 ! P6 LTYPE=60 :TRACTION IN NORMAL-DIRECTION FOR FACE-6
229 ! I/F VARIABLES
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(:)
234 real(kind=kreal) rho
235 integer(kind=kint) ltype,nsize
236
237 ! LOCAL VARIABLES
238 integer(kind=kint) ndof
239 parameter(ndof=3)
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)
251
252 ax = 0.0d0; ay = 0.0d0; az = 0.0d0; rx = 0.0d0; ry = 0.0d0; rz = 0.0d0;
253 !
254 ! SET VALUE
255 !
256 val = params(0)
257 !
258 ! SELCTION OF LOAD TYPE
259 !
260 ivol=0
261 isuf=0
262 if( ltype.LT.10 ) then
263 ivol=1
264 if( ltype.EQ.5 ) then
265 ax=params(1)
266 ay=params(2)
267 az=params(3)
268 rx=params(4)
269 ry=params(5)
270 rz=params(6)
271 endif
272 else if( ltype.GE.10 ) then
273 isuf=1
274 call getsubface( etype, ltype/10, surtype, nod )
275 nsur = getnumberofnodes( surtype )
276 endif
277 ! CLEAR VECT
278 nsize=nn*ndof
279 vect(1:nsize)=0.0d0
280 !** SURFACE LOAD
281 if( isuf==1 ) then
282 ! INTEGRATION OVER SURFACE
283 do i=1,nsur
284 elecoord(1,i)=xx(nod(i))
285 elecoord(2,i)=yy(nod(i))
286 elecoord(3,i)=zz(nod(i))
287 enddo
288 do ig2=1,numofquadpoints( surtype )
289 call getquadpoint( surtype, ig2, localcoord(1:2) )
290 call getshapefunc( surtype, localcoord(1:2), h(1:nsur) )
291
292 wg=getweight( surtype, ig2 )
293 normal=surfacenormal( surtype, nsur, localcoord(1:2), elecoord(:,1:nsur) )
294 do i=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)
298 enddo
299 enddo
300 endif
301 !** VOLUME LOAD
302 if( ivol==1 ) then
303 plx(:)=0.0d0
304 ply(:)=0.0d0
305 plz(:)=0.0d0
306 ! LOOP FOR INTEGRATION POINTS
307 do lx=1,numofquadpoints( etype )
308 call getquadpoint( etype, lx, localcoord )
309 call getshapefunc( etype, localcoord, h(1:nn) )
310 call getshapederiv( etype, localcoord, deriv )
311 ! JACOBI MATRIX
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) )
315 !DETERMINANT OF JACOBIAN
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)
322
323 coefx=1.0
324 coefy=1.0
325 coefz=1.0
326 ! CENTRIFUGAL LOAD
327 if( ltype==5 ) then
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
334 phx=xcod-hx
335 phy=ycod-hy
336 phz=zcod-hz
337 coefx=rho*val*val*phx
338 coefy=rho*val*val*phy
339 coefz=rho*val*val*phz
340 end if
341
342 wg=getweight( etype, lx )*det
343 do i=1,nn
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
347 enddo
348 enddo
349 if( ltype.EQ.1) then
350 do i=1,nn
351 vect(3*i-2)=val*plx(i)
352 enddo
353 else if( ltype.EQ.2 ) then
354 do i=1,nn
355 vect(3*i-1)=val*ply(i)
356 enddo
357 else if( ltype.EQ.3 ) then
358 do i=1,nn
359 vect(3*i )=val*plz(i)
360 enddo
361 else if( ltype.EQ.4 ) then
362 vx=params(1)
363 vy=params(2)
364 vz=params(3)
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)
368 do i=1,nn
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
372 enddo
373 else if( ltype.EQ.5 ) then
374 do i=1,nn
375 vect(3*i-2)=plx(i)
376 vect(3*i-1)=ply(i)
377 vect(3*i )=plz(i)
378 enddo
379 end if
380 endif
381
382 end subroutine dl_c3
383
385 !----------------------------------------------------------------------*
386 subroutine tload_c3 &
387 (etype, nn, xx, yy, zz, tt, t0, gausses, &
388 vect, cdsys_id, coords)
389 !----------------------------------------------------------------------*
390
391 use m_fstr
392 use mmechgauss
393 use m_matmatrix
394 use m_utilities
395
396 !---------------------------------------------------------------------
397
398 integer(kind=kint), parameter :: ndof = 3
399 integer(kind=kint), intent(in) :: etype, nn
400 type(tgaussstatus), intent(in) :: gausses(:)
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)
406
407 !---------------------------------------------------------------------
408
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
417
418 !---------------------------------------------------------------------
419
420 matlaniso = .false.
421
422 if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
423 ina = tt(1)
424 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
425 if( .not. ierr ) matlaniso = .true.
426 end if
427
428 vect(:) = 0.0d0
429
430 ecoord(1, :) = xx(:)
431 ecoord(2, :) = yy(:)
432 ecoord(3, :) = zz(:)
433
434 ! LOOP FOR INTEGRATION POINTS
435 do lx = 1, numofquadpoints(etype)
436
437 call getquadpoint( etype, lx, naturalcoord(:) )
438 call getshapefunc( etype, naturalcoord, h(1:nn) )
439 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv )
440
441 if( matlaniso ) 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"
446 end if
447 end if
448
449 ! WEIGHT VALUE AT GAUSSIAN POINT
450 wg = getweight(etype, lx)*det
451 b(1:6,1:nn*ndof)=0.0d0
452 do j=1,nn
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)
462 enddo
463
464 tempc = dot_product( h(1:nn), tt(1:nn) )
465 temp0 = dot_product( h(1:nn), t0(1:nn) )
466
467 call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
468
469 ina(1) = tempc
470 if( matlaniso ) then
471 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
472 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
473 else
474 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
475 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
476 alp = outa(1)
477 end if
478 ina(1) = temp0
479 if( matlaniso ) then
480 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
481 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
482 else
483 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
484 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
485 alp0 = outa(1)
486 end if
487
488 !**
489 !** THERMAL strain
490 !**
491 if( matlaniso ) then
492 do j=1,3
493 epsth(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
494 end do
495 epsth(4:6) = 0.0d0
496 call transformation(coordsys, tm)
497 epsth(:) = matmul( epsth(:), tm ) ! to global coord
498 epsth(4:6) = epsth(4:6)*2.0d0
499 else
500 thermal_eps=alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
501 epsth(1:3) = thermal_eps
502 epsth(4:6) = 0.0d0
503 end if
504
505 !**
506 !** SET SGM {s}=[D]{e}
507 !**
508 sgm(:) = matmul( d(:, :), epsth(:) )
509 !**
510 !** CALCULATE LOAD {F}=[B]T{e}
511 !**
512 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:),b(:, :) )*wg
513
514 end do
515
516 end subroutine tload_c3
517
518 subroutine cal_thermal_expansion_c3( tt0, ttc, material, coordsys, matlaniso, EPSTH )
519 use m_fstr
520 use m_utilities
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)
527
528 integer(kind=kint) :: j
529 real(kind=kreal) :: ina(1), outa(1)
530 logical :: ierr
531 real(kind=kreal) :: alp, alp0, alpo(3), alpo0(3), tm(6,6)
532
533 ina(1) = ttc
534 if( matlaniso ) then
535 call fetch_tabledata( mc_orthoexp, material%dict, alpo(:), ierr, ina )
536 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
537 else
538 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
539 if( ierr ) outa(1) = material%variables(m_exapnsion)
540 alp = outa(1)
541 end if
542 ina(1) = tt0
543 if( matlaniso ) then
544 call fetch_tabledata( mc_orthoexp, material%dict, alpo0(:), ierr, ina )
545 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
546 else
547 call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
548 if( ierr ) outa(1) = material%variables(m_exapnsion)
549 alp0 = outa(1)
550 end if
551 if( matlaniso ) then
552 do j=1,3
553 epsth(j) = alpo(j)*(ttc-ref_temp)-alpo0(j)*(tt0-ref_temp)
554 end do
555 call transformation( coordsys(:, :), tm)
556 epsth(:) = matmul( epsth(:), tm ) ! to global coord
557 epsth(4:6) = epsth(4:6)*2.0d0
558 else
559 epsth(1:3)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
560 end if
561 end subroutine
562
563 !Hughes, T. J. R., and J. Winget, ?gFinite Rotation Effects in Numerical Integration of
564 ! Rate Constitutive Equations Arising in Large Deformation Analysis,?h
565 ! International Journal for Numerical Methods in Engineering, vol. 15, pp. 1862-1867, 1980.
566 subroutine hughes_winget_rotation_3d( rot, stress_in, stress_out )
567 use m_utilities
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)
571
572 real(kind=kreal) :: dr(3,3)
573
574 !calc dR=(I-0.5*dW)^-1*(I+0.5*dW)
575 dr = i3-0.5d0*rot
576 call calinverse(3, dr)
577 dr = matmul(dr,i3+0.5d0*rot)
578
579 stress_out = matmul(dr,stress_in)
580 stress_out = matmul(stress_out,transpose(dr))
581 end subroutine
582
583 subroutine update_stress3d( flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn )
584 use m_fstr
585 use m_matmatrix
586 use mmechgauss
587 use m_utilities
588
589 type(tgaussstatus), intent(inout) :: gauss
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) !deformation gradient (used for ss_out)
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
600
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
607
608 mtype = gauss%pMaterial%mtype
609
610 if( iselastoplastic(mtype) .OR. mtype == norton )then
611 isep = 1
612 else
613 isep = 0
614 endif
615
616 if( present(ttc) .AND. present(ttn) ) then
617 call matlmatrix( gauss, d3, d, time, tincr, coordsys, ttc, isep )
618 else
619 call matlmatrix( gauss, d3, d, time, tincr, coordsys, isep=isep)
620 end if
621
622 if( flag == infinitesimal ) then
623
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 )
628 else
629 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
630 end if
631 gauss%stress = real(gauss%stress)
632 end if
633
634 else if( flag == totallag ) then
635
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 )
642 else
643 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
644 end if
645 else
646 gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
647 end if
648
649 else if( flag == updatelag ) then
650 ! CALL GEOMAT_C3( gauss%stress, mat )
651 ! D(:, :) = D(:, :)+mat(:, :)
652
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 )
656 else
657 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
658 end if
659 else
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)
667
668 !stress integration
669 trd = dstrain(1)+dstrain(2)+dstrain(3)
670 dum(:,:) = dumstress + matmul( rot,dumstress ) - matmul( dumstress, rot ) + dumstress*trd
671 !call Hughes_Winget_rotation_3D( rot, dumstress, dum )
672
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)
679
680 if( mtype == usermaterial ) then
681 call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys )
682 else if( mtype == norton ) then
683 !gauss%pMaterial%mtype = mtype
684 if( tincr /= 0.0d0 .AND. any( gauss%stress /= 0.0d0 ) ) then
685 !gauss%pMaterial%mtype = mtype
686 if( present(ttc) .AND. present(ttn) ) then
687 call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr, ttc, ttn )
688 else
689 call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr )
690 end if
691 end if
692 end if
693 end if
694
695 end if
696
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 )
701 else
702 call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
703 gauss%istatus(1), gauss%fstatus )
704 end if
705 end if
706
707 !convert stress/strain measure for output
708 if( opsstype == kopss_solution ) then
709
710 if( flag == infinitesimal ) then !linear
711 gauss%stress_out(1:6) = gauss%stress(1:6)
712 gauss%strain_out(1:6) = gauss%strain(1:6)
713 else !nonlinear
714 !convert stress
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)
722
723 det = determinant33(f)
724 if( det == 0.d0 ) stop "Fail to convert stress: detF=0"
725 ! cauchy stress = (1/detF)*F*(2ndPK stress)*F^T
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))
728
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)
737 endif
738
739 !calc logarithmic strain
740 dum(1:3,1:3) = matmul(f(1:3,1:3),transpose(f(1:3,1:3)))
741 tensor(1) = dum(1,1)
742 tensor(2) = dum(2,2)
743 tensor(3) = dum(3,3)
744 tensor(4) = dum(1,2)
745 tensor(5) = dum(2,3)
746 tensor(6) = dum(3,1)
747 call get_principal(tensor, eigval, princ)
748
749 do k=1,3
750 if( eigval(k) <= 0.d0 ) stop "Fail to calc log strain: stretch<0"
751 eigval(k) = 0.5d0*dlog(eigval(k)) !log(sqrt(lambda))
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
755 end do
756 do i=1,3
757 do j=1,3
758 dum(i,j) = 0.d0
759 do k=1,3
760 dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
761 end do
762 end do
763 end do
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)
770 endif
771
772 else
773 gauss%stress_out(1:6) = gauss%stress(1:6)
774 gauss%strain_out(1:6) = gauss%strain(1:6)
775 end if
776
777 end subroutine
778
780 !---------------------------------------------------------------------*
781 subroutine update_c3 &
782 (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
783 gausses, iter, time, tincr, tt, t0, tn)
784 !---------------------------------------------------------------------*
785
786 use m_fstr
787 use mmaterial
788 use mmechgauss
789 use m_matmatrix
791 use m_utilities
792
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)
801 type(tgaussstatus), intent(inout) :: gausses(:)
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)
808
809 ! LCOAL VARIAVLES
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
820
821 qf(:) = 0.0d0
822 ! we suppose the same material type in the element
823 flag = gausses(1)%pMaterial%nlgeom_flag
824 elem(:,:) = ecoord(:,:)
825 totaldisp(:,:) = u(:,:)+ ddu(:,:)
826 if( flag == updatelag ) then
827 elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
828 elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
829 ! elem = elem1
830 totaldisp(:,:) = ddu(:,:)
831 end if
832
833 matlaniso = .false.
834 if( present(tt) .and. cdsys_id > 0 ) then
835 ina = tt(1)
836 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
837 if( .not. ierr ) matlaniso = .true.
838 end if
839
840 do lx = 1, numofquadpoints(etype)
841
842 call getquadpoint( etype, lx, naturalcoord(:) )
843 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
844
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"
850 end if
851 end if
852
853 ! ========================================================
854 ! UPDATE STRAIN and STRESS
855 ! ========================================================
856
857 ! Thermal Strain
858 epsth = 0.0d0
859 if( present(tt) .AND. present(t0) ) then
860 call getshapefunc(etype, naturalcoord, spfunc)
861 ttc = dot_product(tt, spfunc)
862 tt0 = dot_product(t0, spfunc)
863 ttn = dot_product(tn, spfunc)
864 call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
865 end if
866
867 ! Update strain
868 ! Small strain
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(:) ! allright?
877
878 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
879 if( flag == infinitesimal ) then
880 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
881
882 else if( flag == totallag ) then
883 ! Green-Lagrange strain
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) )
893
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)
896
897 else if( flag == updatelag ) then
898 rot = 0.0d0
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)
902
903 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
904
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) )
907
908 end if
909
910 ! Update stress
911 if( present(tt) .AND. present(t0) ) then
912 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
913 else
914 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
915 end if
916
917 ! ========================================================
918 ! calculate the internal force ( equivalent nodal force )
919 ! ========================================================
920 ! Small strain
921 b(1:6, 1:nn*ndof) = 0.0d0
922 do j=1,nn
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)
932 end do
933
934 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
935 if( flag == infinitesimal ) then
936
937 else if( flag == totallag ) then
938
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
941 do j = 1,nn
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)
960 end do
961 ! BL = BL0 + BL1
962 do j=1,nn*ndof
963 b(:,j) = b(:,j)+b1(:,j)
964 end do
965
966 else if( flag == updatelag ) then
967
968 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
969 b(1:6, 1:nn*ndof) = 0.0d0
970 do j = 1, nn
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)
980 end do
981
982 end if
983
984 ! calculate the Internal Force
985 wg=getweight( etype, lx )*det
986 qf(1:nn*ndof) &
987 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
988
989 end do
990
991 end subroutine update_c3
992 !
993 !----------------------------------------------------------------------*
994 subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
995 !----------------------------------------------------------------------*
996 !
997 ! Calculate Strain and Stress increment of solid elements
998 !
999 use mmechgauss
1000
1001 !---------------------------------------------------------------------
1002
1003 integer(kind=kint), intent(in) :: etype, nn
1004 type(tgaussstatus), intent(in) :: gausses(:)
1005 real(kind=kreal), intent(out) :: ndstrain(nn,6)
1006 real(kind=kreal), intent(out) :: ndstress(nn,6)
1007
1008 !---------------------------------------------------------------------
1009
1010 integer :: i, ic
1011 real(kind=kreal) :: temp(12)
1012
1013 !---------------------------------------------------------------------
1014
1015 temp(:) = 0.0d0
1016
1017 ic = numofquadpoints(etype)
1018
1019 do i = 1, ic
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)
1022 end do
1023
1024 temp(1:12) = temp(1:12)/ic
1025
1026 forall( i=1:nn )
1027 ndstrain(i, 1:6) = temp(1:6)
1028 ndstress(i, 1:6) = temp(7:12)
1029 end forall
1030
1031 end subroutine nodalstress_c3
1032
1033
1034 !----------------------------------------------------------------------*
1035 subroutine elementstress_c3(etype, gausses, strain, stress)
1036 !----------------------------------------------------------------------*
1037 !
1038 ! Calculate Strain and Stress increment of solid elements
1039 !
1040 use mmechgauss
1041
1042 !---------------------------------------------------------------------
1043
1044 integer(kind=kint), intent(in) :: etype
1045 type(tgaussstatus), intent(in) :: gausses(:)
1046 real(kind=kreal), intent(out) :: strain(6)
1047 real(kind=kreal), intent(out) :: stress(6)
1048
1049 !---------------------------------------------------------------------
1050
1051 integer :: i, ic
1052
1053 !---------------------------------------------------------------------
1054
1055 strain(:) = 0.0d0; stress(:) = 0.0d0
1056
1057 ic = numofquadpoints(etype)
1058
1059 do i = 1, ic
1060 strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1061 stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1062 enddo
1063
1064 strain(:) = strain(:)/ic
1065 stress(:) = stress(:)/ic
1066
1067 end subroutine elementstress_c3
1068
1069
1071 !----------------------------------------------------------------------*
1072 real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
1073 !----------------------------------------------------------------------*
1074
1075 integer(kind=kint), intent(in) :: etype, nn
1076 real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
1077
1078 !---------------------------------------------------------------------
1079
1080 real(kind=kreal) :: xj(3, 3), det
1081 integer(kind=kint) :: lx
1082 real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1083
1084 !---------------------------------------------------------------------
1085
1086 volume_c3 = 0.0d0
1087
1088 ! LOOP FOR INTEGRATION POINTS
1089 do lx = 1, numofquadpoints(etype)
1090
1091 call getquadpoint(etype, lx, localcoord)
1092 call getshapederiv(etype, localcoord, deriv)
1093
1094 ! JACOBI MATRIX
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) )
1098
1099 ! DETERMINANT OF JACOBIAN
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)
1106
1107 volume_c3 = volume_c3+getweight(etype, lx)*det
1108
1109 end do
1110
1111 end function volume_c3
1112
1113
1114end module m_static_lib_3d
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:126
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:188
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate deivatives of shape fucntion in natural coordiante system.
Definition: element.f90:571
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:863
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kopss_solution
Definition: m_fstr.f90:114
integer(kind=kint) opsstype
Definition: m_fstr.f90:116
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
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.
Definition: utilities.f90:6
subroutine calinverse(nn, a)
calculate inverse of matrix a
Definition: utilities.f90:258
real(kind=kreal) function determinant33(xj)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:374
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13