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