FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_2d.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 use hecmw, only : kint, kreal
9 implicit none
10
11contains
12 !***********************************************************************
13 ! 2D Element:
14 ! STF_C2 :Calculate stiff matrix of 2d elements
15 ! DL_C2 :Deal with DLOAD conditions
16 ! TLOAD_C2 :Deal with thermal expansion force
17 ! UpdateST_C2 : Update Strain and stress
18 !***********************************************************************
19 !----------------------------------------------------------------------*
20 subroutine stf_c2( ETYPE,NN,ecoord,gausses,PARAM1,stiff &
21 ,ISET, u)
22 !----------------------------------------------------------------------*
23 use mmechgauss
24 use m_matmatrix
25 integer(kind=kint), intent(in) :: ETYPE
26 integer(kind=kint), intent(in) :: NN
27 type(tgaussstatus), intent(in) :: gausses(:)
28 real(kind=kreal), intent(in) :: ecoord(2,nn),param1
29 real(kind=kreal), intent(out) :: stiff(:,:)
30 integer(kind=kint),intent(in) :: ISET
31 real(kind=kreal), intent(in), optional :: u(:,:)
32
33 ! LOCAL VARIABLES
34 integer(kind=kint) :: flag
35 integer(kind=kint) NDOF
36 parameter(ndof=2)
37 real(kind=kreal) d(4,4),b(4,ndof*nn),db(4,ndof*nn)
38 real(kind=kreal) h(nn),stress(4)
39 real(kind=kreal) thick,pai
40 real(kind=kreal) det,rr,wg
41 real(kind=kreal) localcoord(2)
42 real(kind=kreal) gderiv(nn,2), cdsys(3,3)
43 integer(kind=kint) J,LX
44 real(kind=kreal) gdispderiv(2,2)
45 real(kind=kreal) b1(4,nn*ndof)
46 real(kind=kreal) smat(4,4)
47 real(kind=kreal) bn(4,nn*ndof), sbn(4,nn*ndof)
48
49 stiff =0.d0
50 flag = gausses(1)%pMaterial%nlgeom_flag
51 ! THICKNESS
52 thick=param1
53 ! FOR AX-SYM. ANALYSIS
54 if(iset==2) then
55 thick=1.d0
56 pai=4.d0*atan(1.d0)
57 endif
58 !* LOOP OVER ALL INTEGRATION POINTS
59 do lx=1,numofquadpoints( etype )
60 call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0,cdsys )
61 if( .not. present(u) ) flag=infinitesimal ! enforce to infinitesimal deformation analysis
62
63 if( flag==1 .and. iset == 2 ) then
64 write(*,'(a)') ' PROGRAM STOP : non-TL element for axixsymmetric element'
65 stop
66 endif
67
68 call getquadpoint( etype, lx, localcoord )
69 call getglobalderiv( etype, nn, localcoord, ecoord, det, gderiv )
70 !
71 if(iset==2) then
72 call getshapefunc( etype, localcoord, h(:) )
73 rr=dot_product( h(1:nn), ecoord(1,1:nn) )
74 wg=getweight( etype, lx )*det*rr*2.d0*pai
75 else
76 rr=thick
77 h(:)=0.d0
78 wg=getweight( etype, lx )*det*rr
79 end if
80 do j=1,nn
81 b(1,2*j-1)=gderiv(j,1)
82 b(2,2*j-1)=0.d0
83 b(3,2*j-1)=gderiv(j,2)
84 b(1,2*j )=0.d0
85 b(2,2*j )=gderiv(j,2)
86 b(3,2*j )=gderiv(j,1)
87 b(4,2*j-1)=h(j)/rr
88 b(4,2*j )=0.d0
89 enddo
90 !
91 ! ----- calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
92 if( flag==1 ) then
93 ! ----- dudx(i,j) ==> gdispderiv(i,j)
94 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof,1:nn), gderiv(1:nn,1:ndof) )
95 b1(1:4,1:nn*ndof) = 0.d0
96 do j=1,nn
97 b1(1,2*j-1) = gdispderiv(1,1)*gderiv(j,1)
98 b1(2,2*j-1) = gdispderiv(1,2)*gderiv(j,2)
99 b1(3,2*j-1) = gdispderiv(1,2)*gderiv(j,1)+gdispderiv(1,1)*gderiv(j,2)
100 b1(1,2*j ) = gdispderiv(2,1)*gderiv(j,1)
101 b1(2,2*j ) = gdispderiv(2,2)*gderiv(j,2)
102 b1(3,2*j ) = gdispderiv(2,2)*gderiv(j,1)+gdispderiv(2,1)*gderiv(j,2)
103 b1(4,2*j-1) = 0.d0
104 b1(4,2*j ) = 0.d0
105 enddo
106 do j=1,nn*ndof
107 b(:,j) = b(:,j)+b1(:,j)
108 enddo
109 endif
110 !
111 db(1:4,1:nn*ndof) = 0.d0
112 db(1:4,1:nn*ndof) = matmul( d(1:4,1:4), b(1:4,1:nn*ndof) )
113 stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
114 matmul( transpose( b(1:4,1:nn*ndof) ), db(1:4,1:nn*ndof) )*wg
115 !
116 ! ----- calculate the stress matrix ( TOTAL LAGRANGE METHOD )
117 if( flag==1 ) then
118 stress(1:4)=gausses(lx)%stress(1:4)
119 bn(1:4,1:nn*ndof) = 0.d0
120 do j=1,nn
121 bn(1,2*j-1) = gderiv(j,1)
122 bn(2,2*j ) = gderiv(j,1)
123 bn(3,2*j-1) = gderiv(j,2)
124 bn(4,2*j ) = gderiv(j,2)
125 enddo
126 smat(:,:) = 0.d0
127 do j=1,2
128 smat(j ,j ) = stress(1)
129 smat(j ,j+2) = stress(3)
130 smat(j+2,j ) = stress(3)
131 smat(j+2,j+2) = stress(2)
132 enddo
133 sbn(1:4,1:nn*ndof) = matmul( smat(1:4,1:4), bn(1:4,1:nn*ndof) )
134 stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
135 matmul( transpose( bn(1:4,1:nn*ndof) ), sbn(1:4,1:nn*ndof) )*wg
136 endif
137 !
138 enddo
139
140 end subroutine stf_c2
141
142
143 !----------------------------------------------------------------------*
144 subroutine dl_c2(ETYPE,NN,XX,YY,RHO,PARAM1,LTYPE,PARAMS,VECT,NSIZE,ISET)
145 !----------------------------------------------------------------------*
146 !**
147 !** Deal with DLOAD conditions
148 !**
149 ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
150 ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
151 ! GRAV LTYPE=4 :GRAVITY FORCE
152 ! CENT LTYPE=5 :CENTRIFUGAL LOAD
153 ! P1 LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
154 ! P2 LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
155 ! P3 LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
156 ! P4 LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
157 ! I/F VARIABLES
158 integer(kind=kint), intent(in) :: ETYPE,NN
159 real(kind=kreal), intent(in) :: xx(nn),yy(nn),params(0:6)
160 real(kind=kreal), intent(out) :: vect(nn*2)
161 real(kind=kreal) rho,param1
162 integer(kind=kint) LTYPE,NSIZE,ISET
163 ! LOCAL VARIABLES
164 integer(kind=kint), parameter :: NDOF=2
165 real(kind=kreal) h(nn)
166 real(kind=kreal) plx(nn),ply(nn)
167 real(kind=kreal) xj(2,2),det,rr,wg
168 real(kind=kreal) pai,val,thick
169 integer(kind=kint) IVOL,ISUF
170 real(kind=kreal) ax,ay,rx,ry
171 real(kind=kreal) normal(2)
172 real(kind=kreal) coefx,coefy,xcod,ycod,hx,hy,phx,phy
173 integer(kind=kint) NOD(NN)
174 integer(kind=kint) LX,I,SURTYPE,NSUR
175 real(kind=kreal) vx,vy,localcoord(2),deriv(nn,2),elecoord(2,nn)
176
177 rx = 0.0d0; ry = 0.0d0
178 ay = 0.0d0; ax = 0.0d0
179
180 ! CONSTANT
181 pai=4.d0*atan(1.d0)
182 ! SET VALUE
183 val=params(0)
184 ! THICKNESS
185 thick=param1
186 ! CLEAR VECT
187 nsize=nn*ndof
188 vect(1:nsize)=0.d0
189 !
190 ! SELCTION OF LOAD TYPE
191 !
192 ivol=0
193 isuf=0
194 if( ltype.LT.10 ) then
195 ivol=1
196 if( ltype.EQ.5 ) then
197 ax=params(1)
198 ay=params(2)
199 rx=params(4)
200 ry=params(5)
201 endif
202 else if( ltype.GE.10 ) then
203 isuf=1
204 call getsubface( etype, ltype/10, surtype, nod )
205 nsur = getnumberofnodes( surtype )
206 endif
207 !** SURFACE LOAD
208 if( isuf==1 ) then
209 do i=1,nsur
210 elecoord(1,i)=xx(nod(i))
211 elecoord(2,i)=yy(nod(i))
212 enddo
213 !** INTEGRATION OVER SURFACE
214 do lx=1,numofquadpoints( surtype )
215 call getquadpoint( surtype, lx, localcoord(1:1) )
216 call getshapefunc( surtype, localcoord(1:1), h(1:nsur) )
217 normal=edgenormal( surtype, nsur, localcoord(1:1), elecoord(:,1:nsur) )
218 wg = getweight( surtype, lx )
219 if( iset==2 ) then
220 rr=0.d0
221 do i=1,nsur
222 rr=rr+h(i)*xx(nod(i))
223 enddo
224 wg=wg*rr*2.d0*pai
225 else
226 wg=wg*thick
227 endif
228 do i=1,nsur
229 vect(2*nod(i)-1)=vect(2*nod(i)-1)+val*wg*h(i)*normal(1)
230 vect(2*nod(i) )=vect(2*nod(i) )+val*wg*h(i)*normal(2)
231 enddo
232 enddo
233 endif
234 !** VOLUME LOAD
235 if( ivol==1 ) then
236 plx(:)=0.d0
237 ply(:)=0.d0
238 do lx=1,numofquadpoints( etype )
239 call getquadpoint( etype, lx, localcoord(1:2) )
240 call getshapederiv( etype, localcoord(1:2), deriv )
241 call getshapefunc( etype, localcoord(1:2), h(1:nn) )
242 xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
243 xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
244
245 det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
246
247 wg = getweight( etype, lx )
248 if(iset==2) then
249 rr=dot_product( h(1:nn),xx(1:nn) )
250 wg=wg*det*rr*2.d0*pai
251 else
252 rr=thick
253 wg=wg*det*rr
254 end if
255 coefx=1.d0
256 coefy=1.d0
257 ! CENTRIFUGAL LOAD
258 if( ltype==5 ) then
259 xcod=dot_product( h(1:nn),xx(1:nn) )
260 ycod=dot_product( h(1:nn),yy(1:nn) )
261 hx=ax + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*rx
262 hy=ay + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*ry
263 phx=xcod-hx
264 phy=ycod-hy
265 coefx=rho*val*val*phx
266 coefy=rho*val*val*phy
267 end if
268 do i=1,nn
269 plx(i)=plx(i)+h(i)*wg*coefx
270 ply(i)=ply(i)+h(i)*wg*coefy
271 enddo
272 enddo
273 if( ltype.EQ.1) then
274 do i=1,nn
275 vect(2*i-1)=val*plx(i)
276 vect(2*i )=0.d0
277 enddo
278 else if( ltype.EQ.2 ) then
279 do i=1,nn
280 vect(2*i-1)=0.d0
281 vect(2*i )=val*ply(i)
282 enddo
283 else if( ltype.EQ.4 ) then
284 vx=params(1)
285 vy=params(2)
286 vx=vx/sqrt( params(1)**2 + params(2)**2 )
287 vy=vy/sqrt( params(1)**2 + params(2)**2 )
288 do i=1,nn
289 vect(2*i-1)=val*plx(i)*rho*vx
290 vect(2*i )=val*ply(i)*rho*vy
291 enddo
292 else if( ltype==5 ) then
293 do i=1,nn
294 vect(2*i-1)=plx(i)
295 vect(2*i )=ply(i)
296 enddo
297 end if
298 endif
299 end subroutine dl_c2
300
301 !----------------------------------------------------------------------*
302 subroutine tload_c2(ETYPE,NN,XX,YY,TT,T0,gausses,PARAM1,ISET,VECT)
303 !----------------------------------------------------------------------*
304 !
305 ! THERMAL LOAD OF PLANE ELEMENT
306 !
307 use mmaterial
308 use m_matmatrix
310 ! I/F VARIABLES
311 integer(kind=kint), intent(in) :: ETYPE,NN
312 real(kind=kreal),intent(in) :: xx(nn),yy(nn),tt(nn),t0(nn),param1
313 type(tgaussstatus), intent(in) :: gausses(:)
314 real(kind=kreal),intent(out) :: vect(nn*2)
315 integer(kind=kint) iset
316 ! LOCAL VARIABLES
317 integer(kind=kint) ndof
318 parameter(ndof=2)
319 real(kind=kreal) alp,pp,d(4,4),b(4,ndof*nn)
320 real(kind=kreal) h(nn)
321 real(kind=kreal) eps(4),sgm(4),localcoord(2)
322 real(kind=kreal) deriv(nn,2),dnde(2,nn)
323 real(kind=kreal) xj(2,2),det,rr,dum
324 real(kind=kreal) xji(2,2)
325 real(kind=kreal) pai,thick,wg
326 integer(kind=kint) j,lx
327 real(kind=kreal) tempc,temp0,thermal_eps
328 !*************************
329 ! CONSTANT
330 !*************************
331 pai=4.d0*atan(1.d0)
332 ! CLEAR VECT
333 vect(1:nn*ndof)=0.d0
334 ! THICKNESS
335 thick=param1
336 ! FOR AX-SYM. ANALYSIS
337 if(iset==2) thick=1.d0
338 ! We suppose material properties doesn't varies inside element
339 call calelasticmatrix( gausses(1)%pMaterial, iset, d )
340 alp = gausses(1)%pMaterial%variables(m_exapnsion)
341 pp = gausses(1)%pMaterial%variables(m_poisson)
342 !* LOOP OVER ALL INTEGRATION POINTS
343 do lx=1,numofquadpoints( etype )
344 call getquadpoint( etype, lx, localcoord )
345 call getshapefunc( etype, localcoord, h(1:nn) )
346 call getshapederiv( etype, localcoord, deriv(:,:) )
347 xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
348 xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
349
350 det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
351
352 tempc=dot_product(h(1:nn),tt(1:nn))
353 temp0=dot_product(h(1:nn),t0(1:nn))
354 if(iset==2) then
355 rr=dot_product(h(1:nn),xx(1:nn))
356 wg=getweight( etype, lx )*det*rr*2.d0*pai
357 else
358 rr=thick
359 wg=getweight( etype, lx )*det*rr
360 end if
361 dum=1.d0/det
362 xji(1,1)= xj(2,2)*dum
363 xji(1,2)=-xj(2,1)*dum
364 xji(2,1)=-xj(1,2)*dum
365 xji(2,2)= xj(1,1)*dum
366
367 dnde=matmul( xji, transpose(deriv) )
368 do j=1,nn
369 b(1,2*j-1)=dnde(1,j)
370 b(2,2*j-1)=0.d0
371 b(3,2*j-1)=dnde(2,j)
372 b(1,2*j )=0.d0
373 b(2,2*j )=dnde(2,j)
374 b(3,2*j )=dnde(1,j)
375 b(4,2*j-1)=0.d0
376 b(4,2*j )=0.d0
377 if( iset==2 ) then
378 b(4,2*j-1)=h(j)/rr
379 endif
380 enddo
381 !**
382 !** THERMAL EPS
383 !**
384 thermal_eps=alp*(tempc-temp0)
385 if( iset .EQ. 2 ) then
386 eps(1)=thermal_eps
387 eps(2)=thermal_eps
388 eps(3)=0.d0
389 eps(4)=thermal_eps
390 else if( iset.EQ.0 ) then
391 eps(1)=thermal_eps*(1.d0+pp)
392 eps(2)=thermal_eps*(1.d0+pp)
393 eps(3)=0.d0
394 eps(4)=0.d0
395 else if( iset.EQ.1 ) then
396 eps(1)=thermal_eps
397 eps(2)=thermal_eps
398 eps(3)=0.d0
399 eps(4)=0.d0
400 endif
401 !**
402 !** SET SGM {s}=[D]{e}
403 !**
404 sgm = matmul( eps(1:4), d(1:4,1:4) )
405 !**
406 !** CALCULATE LOAD {F}=[B]T{e}
407 !**
408 vect(1:nn*ndof)=vect(1:nn*ndof)+matmul( sgm(1:4), b(1:4,1:nn*ndof) )*wg
409 enddo
410
411 end subroutine tload_c2
412 !
413 !
415 !---------------------------------------------------------------------*
416 subroutine update_c2( etype,nn,ecoord,gausses,PARAM1,iset, &
417 u, ddu, qf, TT, T0, TN )
418 !---------------------------------------------------------------------*
419 use m_fstr
420 use mmechgauss
421 use m_matmatrix
422 ! I/F VARIAVLES
423 integer(kind=kint), intent(in) :: etype, nn
424 real(kind=kreal), intent(in) :: ecoord(3,nn),param1
425 integer(kind=kint), intent(in) :: iset
426 type(tgaussstatus), intent(inout) :: gausses(:)
427 real(kind=kreal), intent(in) :: u(2,nn)
428 real(kind=kreal), intent(in) :: ddu(2,nn)
429 real(kind=kreal), intent(out) :: qf(:)
430 real(kind=kreal), intent(in), optional :: tt(nn)
431 real(kind=kreal), intent(in), optional :: t0(nn)
432 real(kind=kreal), intent(in), optional :: tn(nn)
433
434
435 integer(kind=kint), parameter :: ndof=2
436 real(kind=kreal) :: d(4,4), b(4,ndof*nn)
437 real(kind=kreal) :: h(nn)
438 real(kind=kreal) :: thick,pai,rr
439 real(kind=kreal) :: det, wg
440 real(kind=kreal) :: localcoord(2)
441 real(kind=kreal) :: gderiv(nn,2), ttc, tt0, ttn
442 real(kind=kreal) :: cdsys(3,3)
443 integer(kind=kint) :: j, LX
444 real(kind=kreal) :: totaldisp(2*nn)
445 real(kind=kreal) :: epsth(4), alp, alp0, ina(1), outa(1)
446 logical :: ierr
447 !
448 qf(:) = 0.d0
449 do j=1,nn
450 totaldisp((j-1)*2+1) = u(1,j) + ddu(1,j)
451 totaldisp(j*2) = u(2,j) + ddu(2,j)
452 enddo
453 !
454 ! THICKNESS
455 thick=param1
456 ! FOR AX-SYM. ANALYSIS
457 if(iset==2) then
458 thick=1.d0
459 pai=4.d0*atan(1.d0)
460 endif
461 !* LOOP OVER ALL INTEGRATION POINTS
462 do lx=1, numofquadpoints(etype)
463 call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0, cdsys )
464
465 call getquadpoint( etype, lx, localcoord(:) )
466 call getglobalderiv( etype, nn, localcoord, ecoord, det, gderiv )
467 !
468 epsth = 0.d0
469 if( present(tt) .AND. present(t0) ) then
470 call getshapefunc( etype, localcoord, h(:) )
471 ttc = dot_product(tt(:), h(:))
472 tt0 = dot_product(t0(:), h(:))
473 ttn = dot_product(tn(:), h(:))
474 ina(1) = ttc
475 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
476 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
477 alp = outa(1)
478 ina(1) = tt0
479 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
480 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
481 alp0 = outa(1)
482 epsth(1:2)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
483 end if
484 !
485 wg=getweight( etype, lx )*det
486 if(iset==2) then
487 call getshapefunc( etype, localcoord, h(:) )
488 rr=dot_product( h(1:nn), ecoord(1,1:nn) )
489 else
490 rr=thick
491 h(:)=0.d0
492 end if
493 do j=1,nn
494 b(1,2*j-1)=gderiv(j,1)
495 b(2,2*j-1)=0.d0
496 b(3,2*j-1)=gderiv(j,2)
497 b(1,2*j )=0.d0
498 b(2,2*j )=gderiv(j,2)
499 b(3,2*j )=gderiv(j,1)
500 b(4,2*j-1)=h(j)/rr
501 b(4,2*j )=0.d0
502 enddo
503
504 gausses(lx)%strain(1:4) = matmul( b(:,:), totaldisp )
505 gausses(lx)%stress(1:4) = matmul( d(1:4, 1:4), gausses(lx)%strain(1:4)-epsth(1:4) )
506
507 !set stress and strain for output
508 gausses(lx)%strain_out(1:4) = gausses(lx)%strain(1:4)
509 gausses(lx)%stress_out(1:4) = gausses(lx)%stress(1:4)
510
511 !
512 ! ----- calculate the Internal Force
513 qf(1:nn*ndof) = qf(1:nn*ndof) + &
514 matmul( gausses(lx)%stress(1:4), b(1:4,1:nn*ndof) )*wg
515 !
516 enddo
517 !
518 end subroutine update_c2
519 !
520 !----------------------------------------------------------------------*
521 subroutine nodalstress_c2(ETYPE,NN,gausses,ndstrain,ndstress)
522 !----------------------------------------------------------------------*
523 !
524 ! Calculate Strain and Stress increment of solid elements
525 !
526 use mmechgauss
527
528 integer(kind=kint), intent(in) :: ETYPE,NN
529 type(tgaussstatus), intent(in) :: gausses(:)
530 real(kind=kreal), intent(out) :: ndstrain(nn,4)
531 real(kind=kreal), intent(out) :: ndstress(nn,4)
532
533 integer :: i,ic
534 real(kind=kreal) :: temp(8)
535
536 temp(:)=0.d0
537 ic = numofquadpoints(etype)
538 do i=1,ic
539 temp(1:4) = temp(1:4) + gausses(i)%strain_out(1:4)
540 temp(5:8) = temp(5:8) + gausses(i)%stress_out(1:4)
541 enddo
542 temp(1:8) = temp(1:8)/ic
543 forall( i=1:nn )
544 ndstrain(i,1:4) = temp(1:4)
545 ndstress(i,1:4) = temp(5:8)
546 end forall
547
548 end subroutine
549
550 !----------------------------------------------------------------------*
551 subroutine elementstress_c2(ETYPE,gausses,strain,stress)
552 !----------------------------------------------------------------------*
553 !
554 ! Calculate Strain and Stress increment of solid elements
555 !
556 use mmechgauss
557 integer(kind=kint), intent(in) :: ETYPE
558 type(tgaussstatus), intent(in) :: gausses(:)
559 real(kind=kreal), intent(out) :: strain(4)
560 real(kind=kreal), intent(out) :: stress(4)
561
562 integer :: i,ic
563
564 strain(:)=0.d0; stress(:)=0.d0
565 ic = numofquadpoints(etype)
566 do i=1,ic
567 strain(:) = strain(:) + gausses(i)%strain_out(1:4)
568 stress(:) = stress(:) + gausses(i)%stress_out(1:4)
569 enddo
570 strain(:) = strain(:)/ic
571 stress(:) = stress(:)/ic
572
573 end subroutine
574
575end module m_static_lib_2d
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
real(kind=kreal) function, dimension(2) edgenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 2d-edge.
Definition: element.f90:898
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
Definition: hecmw.f90:6
This module provides functions for elastic material.
subroutine calelasticmatrix(matl, secttype, d, temp)
Calculate isotropic elastic matrix.
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 Plane deformation elements.
subroutine nodalstress_c2(etype, nn, gausses, ndstrain, ndstress)
subroutine update_c2(etype, nn, ecoord, gausses, param1, iset, u, ddu, qf, tt, t0, tn)
Update strain and stress inside element.
subroutine elementstress_c2(etype, gausses, strain, stress)
subroutine tload_c2(etype, nn, xx, yy, tt, t0, gausses, param1, iset, vect)
subroutine stf_c2(etype, nn, ecoord, gausses, param1, stiff, iset, u)
subroutine dl_c2(etype, nn, xx, yy, rho, param1, ltype, params, vect, nsize, iset)
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter m_exapnsion
Definition: material.f90:97
integer(kind=kint), parameter m_poisson
Definition: material.f90:85
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