1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!>   This module provide common functions of Plane deformation elements
6module m_static_LIB_2d
7  use hecmw, only : kint, kreal
8  use elementInfo
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
309    use m_ElasticLinear
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(:)          !< status of qudrature points
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  !
414  !> Update strain and stress inside element
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)   !< current temperature
431    real(kind=kreal), intent(in), optional :: T0(nn)   !< reference temperature
432    real(kind=kreal), intent(in), optional :: TN(nn)   !< reference temperature
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
576