1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> \brief This module contains several strategy to free locking problem 6!> in Eight-node hexagonal element 7module m_static_LIB_Fbar 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 20 !> This subroutine calculate stiff matrix using F-bar method 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 29 use m_common_struct 30 use m_static_LIB_3d, only: GEOMAT_C3 31 use m_utilities 32 33 !--------------------------------------------------------------------- 34 35 integer(kind=kint), intent(in) :: etype !< element type 36 integer(kind=kint), intent(in) :: nn !< number of elemental nodes 37 real(kind=kreal), intent(in) :: ecoord(3, nn) !< coordinates of elemental nodes 38 type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points 39 real(kind=kreal), intent(out) :: stiff(:,:) !< stiff matrix 40 integer(kind=kint), intent(in) :: cdsys_ID 41 real(kind=kreal), intent(inout) :: coords(3, 3) !< variables to define matreial coordinate system 42 real(kind=kreal), intent(in) :: time !< current time 43 real(kind=kreal), intent(in) :: tincr !< time increment 44 real(kind=kreal), intent(in), optional :: u(:, :) !< nodal displacemwent 45 real(kind=kreal), intent(in), optional :: temperature(nn) !< temperature 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 335 !> Update Strain stress of this element 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 346 use m_ElastoPlastic 347 use mHyperElastic 348 use m_utilities 349 use m_static_LIB_3d 350 351 !--------------------------------------------------------------------- 352 353 integer(kind=kint), intent(in) :: etype !< \param [in] element type 354 integer(kind=kint), intent(in) :: nn !< \param [in] number of elemental nodes 355 real(kind=kreal), intent(in) :: ecoord(3, nn) !< \param [in] coordinates of elemental nodes 356 real(kind=kreal), intent(in) :: u(3, nn) !< \param [in] nodal dislplacements 357 real(kind=kreal), intent(in) :: du(3, nn) !< \param [in] nodal displacement increment 358 integer(kind=kint), intent(in) :: cdsys_ID 359 real(kind=kreal), intent(inout) :: coords(3, 3) !< variables to define matreial coordinate system 360 real(kind=kreal), intent(out) :: qf(nn*3) !< \param [out] Internal Force 361 type(tGaussStatus), intent(inout) :: gausses(:) !< \param [out] status of qudrature points 362 integer, intent(in) :: iter 363 real(kind=kreal), intent(in) :: time !< current time 364 real(kind=kreal), intent(in) :: tincr !< time increment 365 real(kind=kreal), intent(in), optional :: TT(nn) !< current temperature 366 real(kind=kreal), intent(in), optional :: T0(nn) !< reference temperature 367 real(kind=kreal), intent(in), optional :: TN(nn) !< reference temperature 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 662 !> This subroutien calculate thermal loading 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(:) !< status of qudrature points 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) !< variables to define matreial coordinate system 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 826