FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_Fbar.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!-------------------------------------------------------------------------------
8
9 use hecmw, only : kint, kreal
10 use elementinfo
11
12 implicit none
13
14 real(kind=kreal), private, parameter :: i33(3,3) = &
15 & reshape( (/1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/) )
16
17contains
18
19
21 !----------------------------------------------------------------------*
22 subroutine stf_c3d8fbar &
23 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
24 time, tincr, u, temperature)
25 !----------------------------------------------------------------------*
26
27 use mmechgauss
28 use m_matmatrix
30 use m_static_lib_3d, only: geomat_c3
31 use m_utilities
32
33 !---------------------------------------------------------------------
34
35 integer(kind=kint), intent(in) :: etype
36 integer(kind=kint), intent(in) :: nn
37 real(kind=kreal), intent(in) :: ecoord(3, nn)
38 type(tgaussstatus), intent(in) :: gausses(:)
39 real(kind=kreal), intent(out) :: stiff(:,:)
40 integer(kind=kint), intent(in) :: cdsys_id
41 real(kind=kreal), intent(inout) :: coords(3, 3)
42 real(kind=kreal), intent(in) :: time
43 real(kind=kreal), intent(in) :: tincr
44 real(kind=kreal), intent(in), optional :: u(:, :)
45 real(kind=kreal), intent(in), optional :: temperature(nn)
46
47 !---------------------------------------------------------------------
48
49 integer(kind=kint) :: flag
50 integer(kind=kint), parameter :: ndof = 3
51 real(kind=kreal) :: d(6, 6),b(6, ndof*nn),db(6, ndof*nn)
52 real(kind=kreal) :: gderiv(nn, 3),stress(6),mat(6, 6)
53 real(kind=kreal) :: det, wg, temp, spfunc(nn)
54 integer(kind=kint) :: i, j, p, q, lx, serr
55 real(kind=kreal) :: naturalcoord(3)
56 real(kind=kreal) :: gdispderiv(3, 3)
57 real(kind=kreal) :: b1(6, ndof*nn)
58 real(kind=kreal) :: smat(9, 9), elem(3, nn)
59 real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
60 real(kind=kreal) :: coordsys(3, 3)
61
62 real(kind=kreal) :: elem0(3,nn), elem1(3, nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3,2)
63 real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
64 real(kind=kreal) :: gderiv2_ave(ndof*nn,ndof*nn)
65 real(kind=kreal) :: fbar(3,3), jratio(8), coeff, sff, dstrain(6), ddfs, fs(3,3), gfs(3,2)
66
67 !---------------------------------------------------------------------
68
69 stiff(:, :) = 0.0d0
70 ! we suppose the same material type in the element
71 flag = gausses(1)%pMaterial%nlgeom_flag
72 if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
73 elem(:, :) = ecoord(:, :)
74 elem0(:, :) = ecoord(:, :)
75 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
76 elem1(:, :) = ecoord(:, :)+u(:, :)
77
78 !cal volumetric average of J=detF and dN/dx
79 v0 = 0.d0
80 jacob_ave = 0.d0
81 gderiv1_ave(1:nn,1:ndof) = 0.d0
82 gderiv2_ave(1:ndof*nn,1:ndof*nn) = 0.d0
83 do lx = 1, numofquadpoints(etype)
84 call getquadpoint( etype, lx, naturalcoord(:) )
85 call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
86 wg = getweight(etype, lx)*det
87 if( flag == infinitesimal ) then
88 jacob = 1.d0
89 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv(1:nn, 1:ndof)
90 else
91 gdispderiv(1:3, 1:3) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
92 jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
93 jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob)
94 call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
95 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
96 do p=1,nn
97 do q=1,nn
98 do i=1,ndof
99 do j=1,ndof
100 gderiv2_ave(3*p-3+i,3*q-3+j) = gderiv2_ave(3*p-3+i,3*q-3+j) &
101 & + jacob*wg*(gderiv1(p,i)*gderiv1(q,j)-gderiv1(q,i)*gderiv1(p,j))
102 end do
103 end do
104 end do
105 end do
106 endif
107 v0 = v0 + wg
108 jacob_ave = jacob_ave + jacob*wg
109 enddo
110 if( dabs(v0) > 1.d-12 ) then
111 if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
112 jacob_ave = jacob_ave/v0
113 jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
114 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
115 gderiv2_ave(1:ndof*nn,1:ndof*nn) = gderiv2_ave(1:ndof*nn,1:ndof*nn)/(v0*jacob_ave)
116 else
117 stop 'Error in Update_C3D8Fbar: too small element volume'
118 end if
119
120 do lx = 1, numofquadpoints(etype)
121
122 call getquadpoint( etype, lx, naturalcoord(:) )
123 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
124
125 if( cdsys_id > 0 ) then
126 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
127 if( serr == -1 ) stop "Fail to setup local coordinate"
128 if( serr == -2 ) then
129 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
130 end if
131 end if
132
133 if( present(temperature) ) then
134 call getshapefunc( etype, naturalcoord, spfunc )
135 temp = dot_product( temperature, spfunc )
136 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
137 else
138 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
139 end if
140
141 if( flag == updatelag ) then
142 call geomat_c3( gausses(lx)%stress, mat )
143 d(:, :) = d(:, :)-mat
144 endif
145
146 wg = getweight(etype, lx)*det
147 b(1:6, 1:nn*ndof) = 0.0d0
148 do j = 1, nn
149 b(1, 3*j-2) = gderiv(j, 1)
150 b(2, 3*j-1) = gderiv(j, 2)
151 b(3, 3*j ) = gderiv(j, 3)
152 b(4, 3*j-2) = gderiv(j, 2)
153 b(4, 3*j-1) = gderiv(j, 1)
154 b(5, 3*j-1) = gderiv(j, 3)
155 b(5, 3*j ) = gderiv(j, 2)
156 b(6, 3*j-2) = gderiv(j, 3)
157 b(6, 3*j ) = gderiv(j, 1)
158 end do
159
160 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
161 if( flag == infinitesimal ) then
162 b2(1:6, 1:nn*ndof) = 0.0d0
163 do j = 1, nn
164 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
165 b2(1,3*j-2:3*j) = z1(1:3,1)
166 b2(2,3*j-2:3*j) = z1(1:3,1)
167 b2(3,3*j-2:3*j) = z1(1:3,1)
168 end do
169
170 ! BL = Jratio*(BL0 + BL1)+BL2
171 do j = 1, nn*ndof
172 b(1:3, j) = b(1:3,j)+b2(1:3,j)
173 end do
174
175 else if( flag == totallag ) then
176 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
177 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
178 call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
179
180 ! ---dudx(i,j) ==> gdispderiv(i,j)
181 b1(1:6, 1:nn*ndof) = 0.0d0
182 do j = 1, nn
183 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
184 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
185 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
186 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
187 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
188 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
189 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
190 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
191 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
192 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
193 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
194 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
195 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
196 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
197 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
198 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
199 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
200 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
201 end do
202
203 dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
204 dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
205 dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
206 dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
207 dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
208 dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
209
210 b2(1:6, 1:nn*ndof) = 0.0d0
211 do j = 1, nn
212 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
213 b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3,1)
214 b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3,1)
215 b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3,1)
216 b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3,1)
217 b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3,1)
218 b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3,1)
219 end do
220
221 ! BL = Jratio*(BL0 + BL1)+BL2
222 do j = 1, nn*ndof
223 b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
224 end do
225
226 else if( flag == updatelag ) then
227 wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
228
229 b2(1:3, 1:nn*ndof) = 0.0d0
230 do j = 1, nn
231 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
232 b2(1, 3*j-2:3*j) = z1(1:3,1)
233 b2(2, 3*j-2:3*j) = z1(1:3,1)
234 b2(3, 3*j-2:3*j) = z1(1:3,1)
235 end do
236
237 do j = 1, nn*ndof
238 b(1:3, j) = b(1:3,j)+b2(1:3,j)
239 end do
240
241 end if
242
243 db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
244 forall( i=1:nn*ndof, j=1:nn*ndof )
245 stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
246 end forall
247
248 ! calculate the initial stress matrix(1): dFbar*dFbar*Stress
249 if( flag == totallag .or. flag == updatelag ) then
250 stress(1:6) = gausses(lx)%stress
251 if( flag == totallag ) then
252 coeff = jratio(lx)
253 sff = dot_product(stress(1:6),dstrain(1:6))
254 else if( flag == updatelag ) then
255 coeff = 1.d0
256 sff = stress(1)+stress(2)+stress(3)
257 gderiv1 = gderiv
258 fbar(1:3,1:3) = i33(1:3,1:3)
259 end if
260
261 bn(1:9, 1:nn*ndof) = 0.0d0
262 do j = 1, nn
263 bn(1, 3*j-2) = coeff*gderiv(j, 1)
264 bn(2, 3*j-1) = coeff*gderiv(j, 1)
265 bn(3, 3*j ) = coeff*gderiv(j, 1)
266 bn(4, 3*j-2) = coeff*gderiv(j, 2)
267 bn(5, 3*j-1) = coeff*gderiv(j, 2)
268 bn(6, 3*j ) = coeff*gderiv(j, 2)
269 bn(7, 3*j-2) = coeff*gderiv(j, 3)
270 bn(8, 3*j-1) = coeff*gderiv(j, 3)
271 bn(9, 3*j ) = coeff*gderiv(j, 3)
272 z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
273 bn(1, 3*j-2:3*j) = bn(1, 3*j-2:3*j) + fbar(1,1)*z1(1:3,1)
274 bn(2, 3*j-2:3*j) = bn(2, 3*j-2:3*j) + fbar(2,1)*z1(1:3,1)
275 bn(3, 3*j-2:3*j) = bn(3, 3*j-2:3*j) + fbar(3,1)*z1(1:3,1)
276 bn(4, 3*j-2:3*j) = bn(4, 3*j-2:3*j) + fbar(1,2)*z1(1:3,1)
277 bn(5, 3*j-2:3*j) = bn(5, 3*j-2:3*j) + fbar(2,2)*z1(1:3,1)
278 bn(6, 3*j-2:3*j) = bn(6, 3*j-2:3*j) + fbar(3,2)*z1(1:3,1)
279 bn(7, 3*j-2:3*j) = bn(7, 3*j-2:3*j) + fbar(1,3)*z1(1:3,1)
280 bn(8, 3*j-2:3*j) = bn(8, 3*j-2:3*j) + fbar(2,3)*z1(1:3,1)
281 bn(9, 3*j-2:3*j) = bn(9, 3*j-2:3*j) + fbar(3,3)*z1(1:3,1)
282 end do
283 smat(:, :) = 0.0d0
284 do j= 1, 3
285 smat(j , j ) = stress(1)
286 smat(j , j+3) = stress(4)
287 smat(j , j+6) = stress(6)
288 smat(j+3, j ) = stress(4)
289 smat(j+3, j+3) = stress(2)
290 smat(j+3, j+6) = stress(5)
291 smat(j+6, j ) = stress(6)
292 smat(j+6, j+3) = stress(5)
293 smat(j+6, j+6) = stress(3)
294 end do
295 sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
296 forall( i=1:nn*ndof, j=1:nn*ndof )
297 stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
298 end forall
299
300 ! calculate the initial stress matrix(2): d(dFbar)*Stress
301 fs(1,1) = fbar(1,1)*stress(1)+fbar(1,2)*stress(4)+fbar(1,3)*stress(6)
302 fs(1,2) = fbar(1,1)*stress(4)+fbar(1,2)*stress(2)+fbar(1,3)*stress(5)
303 fs(1,3) = fbar(1,1)*stress(6)+fbar(1,2)*stress(5)+fbar(1,3)*stress(3)
304 fs(2,1) = fbar(2,1)*stress(1)+fbar(2,2)*stress(4)+fbar(2,3)*stress(6)
305 fs(2,2) = fbar(2,1)*stress(4)+fbar(2,2)*stress(2)+fbar(2,3)*stress(5)
306 fs(2,3) = fbar(2,1)*stress(6)+fbar(2,2)*stress(5)+fbar(2,3)*stress(3)
307 fs(3,1) = fbar(3,1)*stress(1)+fbar(3,2)*stress(4)+fbar(3,3)*stress(6)
308 fs(3,2) = fbar(3,1)*stress(4)+fbar(3,2)*stress(2)+fbar(3,3)*stress(5)
309 fs(3,3) = fbar(3,1)*stress(6)+fbar(3,2)*stress(5)+fbar(3,3)*stress(3)
310 do i=1,nn
311 z1(1:3,1) = (gderiv1_ave(i,1:3)-gderiv1(i,1:3))/3.d0
312 gfs(1:3,1) = coeff*matmul(fs,gderiv(i,1:3))
313 do j=1,nn
314 z1(1:3,2) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
315 gfs(1:3,2) = coeff*matmul(fs,gderiv(j,1:3))
316 do p=1,ndof
317 do q=1,ndof
318 ddfs = z1(p,1)*z1(q,2)
319 ddfs = ddfs + (gderiv2_ave(3*i-3+p,3*j-3+q)-gderiv1_ave(i,p)*gderiv1_ave(j,q))/3.d0
320 ddfs = ddfs + gderiv1(i,q)*gderiv1(j,p)/3.d0
321 ddfs = sff*ddfs + z1(p,1)*gfs(q,2)+z1(q,2)*gfs(p,1)
322 stiff(3*i-3+p, 3*j-3+q) = stiff(3*i-3+p, 3*j-3+q) + ddfs*wg
323 end do
324 end do
325 end do
326 end do
327 end if
328
329
330 end do ! gauss roop
331
332 end subroutine stf_c3d8fbar
333
334
336 !----------------------------------------------------------------------*
337 subroutine update_c3d8fbar &
338 (etype, nn, ecoord, u, du, cdsys_id, coords, &
339 qf ,gausses, iter, time, tincr, tt,t0, tn )
340 !----------------------------------------------------------------------*
341
342 use m_fstr
343 use mmaterial
344 use mmechgauss
345 use m_matmatrix
347 use mhyperelastic
348 use m_utilities
350
351 !---------------------------------------------------------------------
352
353 integer(kind=kint), intent(in) :: etype
354 integer(kind=kint), intent(in) :: nn
355 real(kind=kreal), intent(in) :: ecoord(3, nn)
356 real(kind=kreal), intent(in) :: u(3, nn)
357 real(kind=kreal), intent(in) :: du(3, nn)
358 integer(kind=kint), intent(in) :: cdsys_id
359 real(kind=kreal), intent(inout) :: coords(3, 3)
360 real(kind=kreal), intent(out) :: qf(nn*3)
361 type(tgaussstatus), intent(inout) :: gausses(:)
362 integer, intent(in) :: iter
363 real(kind=kreal), intent(in) :: time
364 real(kind=kreal), intent(in) :: tincr
365 real(kind=kreal), intent(in), optional :: tt(nn)
366 real(kind=kreal), intent(in), optional :: t0(nn)
367 real(kind=kreal), intent(in), optional :: tn(nn)
368
369 !---------------------------------------------------------------------
370
371 integer(kind=kint) :: flag
372 integer(kind=kint), parameter :: ndof = 3
373 real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
374 real(kind=kreal) :: gderiv(nn, 3), gdispderiv(3, 3), det, wg
375 integer(kind=kint) :: j, lx, serr
376 real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
377 real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
378 real(kind=kreal) :: dstrain(6)
379 real(kind=kreal) :: dvol
380 real(kind=kreal) :: ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
381 logical :: ierr, matlaniso
382
383 real(kind=kreal) :: elem0(3,nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3)
384 real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
385 real(kind=kreal) :: fbar(3,3), jratio(8)
386 real(kind=kreal) :: jacob_ave05, gderiv05_ave(nn,ndof)
387
388 !---------------------------------------------------------------------
389
390 qf(:) = 0.0d0
391 ! we suppose the same material type in the element
392 flag = gausses(1)%pMaterial%nlgeom_flag
393 elem0(1:ndof,1:nn) = ecoord(1:ndof,1:nn)
394 totaldisp(:, :) = u(:, :)+du(:, :)
395 if( flag == infinitesimal ) then
396 elem(:, :) = ecoord(:, :)
397 elem1(:, :) = ecoord(:, :)
398 else if( flag == totallag ) then
399 elem(:, :) = ecoord(:, :)
400 elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
401 else if( flag == updatelag ) then
402 elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
403 elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
404 totaldisp(:, :) = du(:, :)
405 end if
406
407 matlaniso = .false.
408 if( cdsys_id > 0 .AND. present(tt) ) then
409 ina = tt(1)
410 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
411 if( .not. ierr ) matlaniso = .true.
412 end if
413
414 !cal volumetric average of J=detF and dN/dx
415 v0 = 0.d0
416 jacob_ave = 0.d0
417 gderiv1_ave(1:nn,1:ndof) = 0.d0
418 if( flag == updatelag ) then
419 jacob_ave05 = 0.d0
420 gderiv05_ave(1:nn,1:ndof) = 0.d0
421 endif
422 do lx = 1, numofquadpoints(etype)
423 call getquadpoint( etype, lx, naturalcoord(:) )
424 call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
425 wg = getweight(etype, lx)*det
426 if( flag == infinitesimal ) then
427 jacob = 1.d0
428 gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
429 else
430 gdispderiv(1:3, 1:3) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
431 jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
432 jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob) ! Jratio(LX) = jacob**(-1.d0/3.d0)
433
434 call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
435 endif
436 v0 = v0 + wg
437 jacob_ave = jacob_ave + jacob*wg
438 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
439 if( flag == updatelag ) then
440 call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv1)
441 wg = getweight(etype, lx)*det
442 jacob_ave05 = jacob_ave05 + wg
443 gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof) + wg*gderiv1(1:nn, 1:ndof)
444 endif
445 enddo
446 if( dabs(v0) > 1.d-12 ) then
447 if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
448 jacob_ave = jacob_ave/v0
449 jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
450 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
451 if( flag == updatelag ) gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof)/jacob_ave05
452 else
453 stop 'Error in Update_C3D8Fbar: too small element volume'
454 end if
455
456 do lx = 1, numofquadpoints(etype)
457
458 call getquadpoint( etype, lx, naturalcoord(:) )
459 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
460
461 if( cdsys_id > 0 ) then
462 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
463 if( serr == -1 ) stop "Fail to setup local coordinate"
464 if( serr == -2 ) then
465 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
466 end if
467 end if
468
469 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
470
471 ! ========================================================
472 ! UPDATE STRAIN and STRESS
473 ! ========================================================
474
475 ! Thermal Strain
476 epsth = 0.0d0
477 if( present(tt) .AND. present(t0) ) then
478 call getshapefunc(etype, naturalcoord, spfunc)
479 ttc = dot_product(tt, spfunc)
480 tt0 = dot_product(t0, spfunc)
481 ttn = dot_product(tn, spfunc)
482 call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
483 end if
484
485 ! Update strain
486 if( flag == infinitesimal ) then
487 dvol = dot_product( totaldisp(1,1:nn), gderiv1_ave(1:nn,1) ) !du1/dx1
488 dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv1_ave(1:nn,2) ) !du2/dx2
489 dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv1_ave(1:nn,3) ) !du3/dx3
490 dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
491 dstrain(1) = gdispderiv(1, 1) + dvol
492 dstrain(2) = gdispderiv(2, 2) + dvol
493 dstrain(3) = gdispderiv(3, 3) + dvol
494 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
495 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
496 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
497 dstrain(:) = dstrain(:)-epsth(:)
498
499 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
500
501 else if( flag == totallag ) then
502 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
503
504 ! Green-Lagrange strain
505 dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
506 dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
507 dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
508 dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
509 dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
510 dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
511 dstrain(:) = dstrain(:)-epsth(:)
512
513 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
514
515 else if( flag == updatelag ) then
516 dvol = dot_product( totaldisp(1,1:nn), gderiv05_ave(1:nn,1) ) !du1/dx1
517 dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv05_ave(1:nn,2) ) !du2/dx2
518 dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv05_ave(1:nn,3) ) !du3/dx3
519 dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
520 dstrain(1) = gdispderiv(1, 1) + dvol
521 dstrain(2) = gdispderiv(2, 2) + dvol
522 dstrain(3) = gdispderiv(3, 3) + dvol
523 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
524 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
525 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
526 dstrain(:) = dstrain(:)-epsth(:)
527
528 rot = 0.0d0
529 rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
530 rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
531 rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
532
533 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
534
535 call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv1)
536 gdispderiv(1:ndof, 1:ndof) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
537 fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
538
539 end if
540
541 ! Update stress
542 if( present(tt) .AND. present(t0) ) then
543 call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr, ttc, tt0, ttn )
544 else
545 call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr )
546 end if
547
548 ! ========================================================
549 ! calculate the internal force ( equivalent nodal force )
550 ! ========================================================
551 ! Small strain
552 b(1:6, 1:nn*ndof) = 0.0d0
553 do j = 1, nn
554 b(1,3*j-2) = gderiv(j, 1)
555 b(2,3*j-1) = gderiv(j, 2)
556 b(3,3*j ) = gderiv(j, 3)
557 b(4,3*j-2) = gderiv(j, 2)
558 b(4,3*j-1) = gderiv(j, 1)
559 b(5,3*j-1) = gderiv(j, 3)
560 b(5,3*j ) = gderiv(j, 2)
561 b(6,3*j-2) = gderiv(j, 3)
562 b(6,3*j ) = gderiv(j, 1)
563 end do
564
565 wg=getweight( etype, lx )*det
566 if( flag == infinitesimal ) then
567 gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
568
569 b2(1:6, 1:nn*ndof) = 0.0d0
570 do j = 1, nn
571 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
572 b2(1,3*j-2:3*j) = z1(1:3)
573 b2(2,3*j-2:3*j) = z1(1:3)
574 b2(3,3*j-2:3*j) = z1(1:3)
575 end do
576
577 ! BL = BL0 + BL2
578 do j = 1, nn*ndof
579 b(:, j) = b(:,j)+b2(:,j)
580 end do
581
582 else if( flag == totallag ) then
583
584 b1(1:6, 1:nn*ndof) = 0.0d0
585 do j = 1, nn
586 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
587 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
588 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
589 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
590 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
591 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
592 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
593 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
594 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
595 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
596 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
597 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
598 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
599 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
600 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
601 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
602 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
603 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
604 end do
605
606 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv1)
607
608 b2(1:6, 1:nn*ndof) = 0.0d0
609 do j = 1, nn
610 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
611 b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3)
612 b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3)
613 b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3)
614 b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3)
615 b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3)
616 b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3)
617 end do
618
619 ! BL = R^2*(BL0 + BL1)+BL2
620 do j = 1, nn*ndof
621 b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
622 end do
623
624 else if( flag == updatelag ) then
625
626 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
627 wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
628 b(1:6, 1:nn*ndof) = 0.0d0
629 do j = 1, nn
630 b(1, 3*j-2) = gderiv(j, 1)
631 b(2, 3*j-1) = gderiv(j, 2)
632 b(3, 3*j ) = gderiv(j, 3)
633 b(4, 3*j-2) = gderiv(j, 2)
634 b(4, 3*j-1) = gderiv(j, 1)
635 b(5, 3*j-1) = gderiv(j, 3)
636 b(5, 3*j ) = gderiv(j, 2)
637 b(6, 3*j-2) = gderiv(j, 3)
638 b(6, 3*j ) = gderiv(j, 1)
639 end do
640
641 b2(1:3, 1:nn*ndof) = 0.0d0
642 do j = 1, nn
643 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
644 b2(1, 3*j-2:3*j) = z1(1:3)
645 b2(2, 3*j-2:3*j) = z1(1:3)
646 b2(3, 3*j-2:3*j) = z1(1:3)
647 end do
648
649 do j = 1, nn*ndof
650 b(1:3, j) = b(1:3,j)+b2(1:3,j)
651 end do
652
653 end if
654
655 qf(1:nn*ndof) &
656 = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
657
658 end do
659
660 end subroutine update_c3d8fbar
661
663 !----------------------------------------------------------------------*
664 subroutine tload_c3d8fbar &
665 (etype, nn, xx, yy, zz, tt, t0, &
666 gausses, vect, cdsys_id, coords)
667 !----------------------------------------------------------------------*
668
669 use m_fstr
670 use mmechgauss
671 use m_matmatrix
672 use m_utilities
673
674 !---------------------------------------------------------------------
675
676 integer(kind=kint), parameter :: ndof = 3
677 integer(kind=kint), intent(in) :: etype, nn
678 type(tgaussstatus), intent(in) :: gausses(:)
679 real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
680 real(kind=kreal), intent(in) :: tt(nn), t0(nn)
681 real(kind=kreal), intent(out) :: vect(nn*ndof)
682 integer(kind=kint), intent(in) :: cdsys_ID
683 real(kind=kreal), intent(inout) :: coords(3, 3)
684
685 !---------------------------------------------------------------------
686
687 real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
688 real(kind=kreal) :: z1(3), det, ecoord(3, nn)
689 integer(kind=kint) :: j, LX, serr
690 real(kind=kreal) :: estrain(6), sgm(6), h(nn)
691 real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
692 real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
693 real(kind=kreal) :: tempc, temp0, v0, jacob_ave, thermal_eps, tm(6,6)
694 logical :: ierr, matlaniso
695
696 !---------------------------------------------------------------------
697
698 matlaniso = .false.
699
700 if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
701 ina = tt(1)
702 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
703 if( .not. ierr ) matlaniso = .true.
704 end if
705
706 vect(:) = 0.0d0
707
708 ecoord(1, :) = xx(:)
709 ecoord(2, :) = yy(:)
710 ecoord(3, :) = zz(:)
711
712 !cal volumetric average of J=detF and dN/dx
713 v0 = 0.d0
714 jacob_ave = 0.d0
715 gderiv1_ave(1:nn,1:ndof) = 0.d0
716 do lx = 1, numofquadpoints(etype)
717 call getquadpoint( etype, lx, naturalcoord(:) )
718 call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv)
719 wg = getweight(etype, lx)*det
720 v0 = v0 + wg
721 jacob_ave = jacob_ave + wg
722 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + wg*gderiv(1:nn, 1:ndof)
723 enddo
724 if( dabs(v0) > 1.d-12 ) then
725 if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in TLOAD_C3D8Fbar: too small average jacob'
726 jacob_ave = jacob_ave/v0
727 gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
728 else
729 stop 'Error in TLOAD_C3D8Fbar: too small element volume'
730 end if
731
732 ! LOOP FOR INTEGRATION POINTS
733 do lx = 1, numofquadpoints(etype)
734
735 call getquadpoint( etype, lx, naturalcoord(:) )
736 call getshapefunc( etype, naturalcoord, h(1:nn) )
737 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
738
739 if( matlaniso ) then
740 call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
741 if( serr == -1 ) stop "Fail to setup local coordinate"
742 if( serr == -2 ) then
743 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
744 end if
745 end if
746
747 ! WEIGHT VALUE AT GAUSSIAN POINT
748 wg = getweight(etype, lx)*det
749 b(1:6, 1:nn*ndof) = 0.0d0
750 do j = 1, nn
751 b(1,3*j-2) = gderiv(j, 1)
752 b(2,3*j-1) = gderiv(j, 2)
753 b(3,3*j ) = gderiv(j, 3)
754 b(4,3*j-2) = gderiv(j, 2)
755 b(4,3*j-1) = gderiv(j, 1)
756 b(5,3*j-1) = gderiv(j, 3)
757 b(5,3*j ) = gderiv(j, 2)
758 b(6,3*j-2) = gderiv(j, 3)
759 b(6,3*j ) = gderiv(j, 1)
760 end do
761
762 do j = 1, nn
763 z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
764 b(1,3*j-2:3*j) = b(1,3*j-2:3*j)+z1(1:3)
765 b(2,3*j-2:3*j) = b(2,3*j-2:3*j)+z1(1:3)
766 b(3,3*j-2:3*j) = b(3,3*j-2:3*j)+z1(1:3)
767 end do
768
769 tempc = dot_product( h(1:nn),tt(1:nn) )
770 temp0 = dot_product( h(1:nn),t0(1:nn) )
771
772 call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
773
774 ina(1) = tempc
775 if( matlaniso ) then
776 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
777 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
778 else
779 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
780 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
781 alp = outa(1)
782 end if
783 ina(1) = temp0
784 if( matlaniso ) then
785 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
786 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
787 else
788 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
789 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
790 alp0 = outa(1)
791 end if
792
793 !**
794 !** THERMAL strain
795 !**
796 if( matlaniso ) then
797 do j = 1, 3
798 estrain(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
799 end do
800 estrain(4:6) = 0.0d0
801 call transformation(coordsys, tm)
802 estrain(:) = matmul( estrain(:), tm ) ! to global coord
803 estrain(4:6) = estrain(4:6)*2.0d0
804 else
805 thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
806 estrain(1:3) = thermal_eps
807 estrain(4:6) = 0.0d0
808 end if
809
810 !**
811 !** SET SGM {s}=[D]{e}
812 !**
813 sgm(:) = matmul( d(:, :), estrain(:) )
814
815 !**
816 !** CALCULATE LOAD {F}=[B]T{e}
817 !**
818 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
819
820 end do
821
822 end subroutine tload_c3d8fbar
823
824
825end module m_static_lib_fbar
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
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
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
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
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, secttype, matrix, time, dtime, cdsys, temperature, isep)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine update_stress3d(flag, gauss, rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn)
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, epsth)
This module contains several strategy to free locking problem in Eight-node hexagonal element.
subroutine stf_c3d8fbar(etype, nn, ecoord, gausses, stiff, cdsys_id, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using F-bar method.
subroutine update_c3d8fbar(etype, nn, ecoord, u, du, cdsys_id, coords, qf, gausses, iter, time, tincr, tt, t0, tn)
Update Strain stress of this element.
subroutine tload_c3d8fbar(etype, nn, xx, yy, zz, tt, t0, gausses, vect, cdsys_id, coords)
This subroutien calculate thermal loading.
This module provides aux functions.
Definition: utilities.f90:6
real(kind=kreal) function determinant33(xj)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
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