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