1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> \brief Eight-node hexagonal element with imcompatible mode 6!> \see R.L.Taylor,P.J.Bereford, and E.L.Wilson, "A Nonconforming element 7!> for Stress Analysis", Intl. J. Numer. Methods Engng, 10(6), pp1211-1219 8!> ,1976 9module m_static_LIB_3dIC 10 11 use hecmw, only : kint, kreal 12 use m_utilities 13 use elementInfo 14 15 implicit none 16 17contains 18 19 !> CALCULATION STIFF Matrix for C3D8IC ELEMENT 20 !----------------------------------------------------------------------* 21 subroutine STF_C3D8IC & 22 (etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, & 23 time, tincr, u, aux, temperature) 24 !----------------------------------------------------------------------* 25 26 use mMechGauss 27 use m_MatMatrix 28 use m_common_struct 29 use m_static_LIB_3d, only: GEOMAT_C3 30 31 !--------------------------------------------------------------------- 32 33 integer(kind=kint), intent(in) :: etype !< element type, not used here 34 integer(kind=kint), intent(in) :: nn !< number of elements nodes 35 real(kind=kreal), intent(in) :: ecoord(3, nn) !< nodal coord of curr element 36 type(tGaussStatus), intent(in) :: gausses(:) !< Info of qudrature points 37 real(kind=kreal), intent(out) :: stiff(:, :) !< stiffness matrix 38 integer(kind=kint), intent(in) :: cdsys_ID 39 real(kind=kreal), intent(inout) :: coords(3, 3) !< variables to define matreial coordinate system 40 real(kind=kreal), intent(in) :: time !< current time 41 real(kind=kreal), intent(in) :: tincr !< time increment 42 real(kind=kreal), intent(in), optional :: u(3, nn) !< nodal displacemwent 43 real(kind=kreal), intent(in), optional :: aux(3, 3) !< enhanced disp of bending mode 44 real(kind=kreal), intent(in), optional :: temperature(nn) !< temperature 45 46 !--------------------------------------------------------------------- 47 48 integer(kind=kint) :: flag 49 integer(kind=kint), parameter :: ndof = 3 50 real(kind=kreal) :: D(6, 6), B(6, ndof*(nn+3)), DB(6, ndof*(nn+3)) 51 real(kind=kreal) :: gderiv(nn+3, 3), stress(6) 52 real(kind=kreal) :: xj(9, 9), jacobian(3, 3), inverse(3, 3) 53 real(kind=kreal) :: tmpstiff((nn+3)*3, (nn+3)*3), tmpk(nn*3, 9) 54 real(kind=kreal) :: det, wg, elem(3, nn), mat(6, 6) 55 integer(kind=kint) :: i, j, LX 56 integer(kind=kint) :: serr 57 real(kind=kreal) :: naturalCoord(3), unode(3, nn+3) 58 real(kind=kreal) :: gdispderiv(3, 3), coordsys(3, 3) 59 real(kind=kreal) :: B1(6, ndof*(nn+3)) 60 real(kind=kreal) :: Smat(9, 9) 61 real(kind=kreal) :: BN(9, ndof*(nn+3)), SBN(9, ndof*(nn+3)) 62 real(kind=kreal) :: spfunc(nn), temp 63 64 !--------------------------------------------------------------------- 65 66 if( present(u) .AND. present(aux) ) then 67 unode(:, 1:nn) = u(:, :) 68 unode(:, nn+1:nn+3) = aux(:, :) 69 end if 70 71 ! we suppose the same material type in the element 72 flag = gausses(1)%pMaterial%nlgeom_flag 73 if( .not. present(u) ) flag = INFINITESIMAL ! enforce to infinitesimal deformation analysis 74 elem(:, :) = ecoord(:, :) 75 if( flag == UPDATELAG ) elem(:, :) = ecoord(:, :)+unode(:, 1:nn) 76 77 ! --- Inverse of Jacobian at elemental center 78 naturalcoord(:) = 0.0D0 79 call getJacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse) 80 inverse(:, :)= inverse(:, :)*det 81 ! ---- We now calculate stiff matrix include imcompatible mode 82 ! [ Kdd Kda ] 83 ! [ Kad Kaa ] 84 tmpstiff(:, :) = 0.0D0 85 B1(1:6, 1:(nn+3)*ndof) = 0.0D0 86 BN(1:9, 1:(nn+3)*ndof) = 0.0D0 87 88 do LX = 1, NumOfQuadPoints(etype) 89 90 call getQuadPoint(etype, LX, naturalCoord) 91 call getGlobalDeriv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) ) 92 93 if( cdsys_ID > 0 ) then 94 call set_localcoordsys( coords, g_LocalCoordSys(cdsys_ID), coordsys(:, :), serr ) 95 if( serr == -1 ) stop "Fail to setup local coordinate" 96 if( serr == -2 ) then 97 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically" 98 end if 99 end if 100 101 if( present(temperature) ) then 102 call getShapeFunc( etype, naturalcoord, spfunc ) 103 temp = dot_product( temperature, spfunc ) 104 call MatlMatrix( gausses(LX), D3, D, time, tincr, coordsys, temp ) 105 else 106 call MatlMatrix( gausses(LX), D3, D, time, tincr, coordsys ) 107 end if 108 109 if( flag == UPDATELAG ) then 110 call GEOMAT_C3( gausses(LX)%stress, mat ) 111 D(:, :) = D(:, :)-mat 112 endif 113 114 ! -- Derivative of shape function of imcompatible mode -- 115 ! [ -2*a 0, 0 ] 116 ! [ 0, -2*b, 0 ] 117 ! [ 0, 0, -2*c ] 118 ! we don't call getShapeDeriv but use above value directly for computation efficiency 119 gderiv(nn+1, :) = -2.0D0*naturalcoord(1)*inverse(1, :)/det 120 gderiv(nn+2, :) = -2.0D0*naturalcoord(2)*inverse(2, :)/det 121 gderiv(nn+3, :) = -2.0D0*naturalcoord(3)*inverse(3, :)/det 122 123 wg = getWeight(etype, LX)*det 124 B(1:6, 1:(nn+3)*ndof)=0.0D0 125 do j = 1, nn+3 126 B(1, 3*j-2) = gderiv(j, 1) 127 B(2, 3*j-1) = gderiv(j, 2) 128 B(3, 3*j ) = gderiv(j, 3) 129 B(4, 3*j-2) = gderiv(j, 2) 130 B(4, 3*j-1) = gderiv(j, 1) 131 B(5, 3*j-1) = gderiv(j, 3) 132 B(5, 3*j ) = gderiv(j, 2) 133 B(6, 3*j-2) = gderiv(j, 3) 134 B(6, 3*j ) = gderiv(j, 1) 135 end do 136 137 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD ) 138 if( flag == TOTALLAG ) then 139 ! ---dudx(i,j) ==> gdispderiv(i,j) 140 gdispderiv(1:ndof,1:ndof) = matmul( unode(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) ) 141 do j = 1, nn+3 142 143 B1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1) 144 B1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1) 145 B1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1) 146 B1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2) 147 B1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2) 148 B1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2) 149 B1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3) 150 B1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3) 151 B1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3) 152 B1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2) 153 B1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2) 154 B1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2) 155 B1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2) 156 B1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2) 157 B1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2) 158 B1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3) 159 B1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3) 160 B1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3) 161 162 end do 163 ! ---BL = BL0 + BL1 164 do j = 1, (nn+3)*ndof 165 B(:, j) = B(:, j)+B1(:, j) 166 end do 167 168 end if 169 170 DB(1:6, 1:(nn+3)*ndof) = matmul( D, B(1:6, 1:(nn+3)*ndof) ) 171 forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof ) 172 tmpstiff(i, j) = tmpstiff(i, j)+dot_product( B(:, i), DB(:, j) )*wg 173 end forall 174 175 ! calculate the stress matrix ( TOTAL LAGRANGE METHOD ) 176 if( flag == TOTALLAG .OR. flag == UPDATELAG ) then 177 stress(1:6) = gausses(LX)%stress 178 do j = 1, nn+3 179 BN(1, 3*j-2) = gderiv(j, 1) 180 BN(2, 3*j-1) = gderiv(j, 1) 181 BN(3, 3*j ) = gderiv(j, 1) 182 BN(4, 3*j-2) = gderiv(j, 2) 183 BN(5, 3*j-1) = gderiv(j, 2) 184 BN(6, 3*j ) = gderiv(j, 2) 185 BN(7, 3*j-2) = gderiv(j, 3) 186 BN(8, 3*j-1) = gderiv(j, 3) 187 BN(9, 3*j ) = gderiv(j, 3) 188 end do 189 Smat(:, :) = 0.0D0 190 do j = 1, 3 191 Smat(j , j ) = stress(1) 192 Smat(j , j+3) = stress(4) 193 Smat(j , j+6) = stress(6) 194 Smat(j+3, j ) = stress(4) 195 Smat(j+3, j+3) = stress(2) 196 Smat(j+3, j+6) = stress(5) 197 Smat(j+6, j ) = stress(6) 198 Smat(j+6, j+3) = stress(5) 199 Smat(j+6, j+6) = stress(3) 200 end do 201 SBN(1:9, 1:(nn+3)*ndof) = matmul( Smat(1:9, 1:9), BN(1:9, 1:(nn+3)*ndof) ) 202 forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof ) 203 tmpstiff(i, j) = tmpstiff(i, j)+dot_product( BN(:, i), SBN(:, j) )*wg 204 end forall 205 end if 206 207 end do 208 209 ! -----Condense tmpstiff to stiff 210 xj(1:9, 1:9)= tmpstiff(nn*ndof+1:(nn+3)*ndof, nn*ndof+1:(nn+3)*ndof) 211 call calInverse(9, xj) 212 tmpk = matmul( tmpstiff( 1:nn*ndof, nn*ndof+1:(nn+3)*ndof ), xj ) 213 stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)-matmul( tmpk, tmpstiff(nn*ndof+1:(nn+3)*ndof, 1:nn*ndof) ) 214 215 end subroutine STF_C3D8IC 216 217 218 !> Update strain and stress inside element 219 !---------------------------------------------------------------------* 220 subroutine UPDATE_C3D8IC & 221 (etype,nn,ecoord, u, du, ddu, cdsys_ID, coords, & 222 qf, gausses, iter, time, tincr, aux, ddaux, TT, T0, TN ) 223 !---------------------------------------------------------------------* 224 225 use m_fstr 226 use mMaterial 227 use mMechGauss 228 use m_MatMatrix 229 use m_ElastoPlastic 230 use m_utilities 231 use m_static_LIB_3d 232 233 integer(kind=kint), intent(in) :: etype !< \param [in] element type 234 integer(kind=kint), intent(in) :: nn !< \param [in] number of elemental nodes 235 real(kind=kreal), intent(in) :: ecoord(3, nn) !< \param [in] coordinates of elemental nodes 236 real(kind=kreal), intent(in) :: u(3, nn) !< \param [in] nodal displacements 237 real(kind=kreal), intent(in) :: du(3, nn) !< \param [in] increment of nodal displacements 238 real(kind=kreal), intent(in) :: ddu(3, nn) !< \param [in] correction of nodal displacements 239 integer(kind=kint), intent(in) :: cdsys_ID 240 real(kind=kreal), intent(inout) :: coords(3, 3) !< variables to define matreial coordinate system 241 real(kind=kreal), intent(out) :: qf(nn*3) !< \param [out] Internal Force 242 type(tGaussStatus), intent(inout) :: gausses(:) !< \param [out] status of qudrature points 243 integer, intent(in) :: iter 244 real(kind=kreal), intent(in) :: time !< current time 245 real(kind=kreal), intent(in) :: tincr !< time increment 246 real(kind=kreal), intent(inout) :: aux(3, 3) !< \param [in] incompatible dof 247 real(kind=kreal), intent(out) :: ddaux(3, 3) !< \param [in] increment of incompatible dof 248 real(kind=kreal), intent(in), optional :: TT(nn) !< current temperature 249 real(kind=kreal), intent(in), optional :: T0(nn) !< reference temperature 250 real(kind=kreal), intent(in), optional :: TN(nn) !< reference temperature 251 252 ! LCOAL VARIAVLES 253 integer(kind=kint) :: flag 254 integer(kind=kint), parameter :: ndof = 3 255 real(kind=kreal) :: D(6,6), B(6,ndof*(nn+3)), B1(6,ndof*(nn+3)), spfunc(nn), ina(1) 256 real(kind=kreal) :: BN(9,ndof*(nn+3)), DB(6, ndof*(nn+3)), stress(6), Smat(9, 9), SBN(9, ndof*(nn+3)) 257 real(kind=kreal) :: gderiv(nn+3,3), gderiv0(nn+3,3), gdispderiv(3,3), F(3,3), det, det0, WG, ttc, tt0, ttn 258 integer(kind=kint) :: i, j, LX, mtype, serr 259 real(kind=kreal) :: naturalCoord(3), rot(3,3), mat(6,6), EPSTH(6) 260 real(kind=kreal) :: totaldisp(3,nn+3), elem(3,nn), elem1(3,nn), coordsys(3,3) 261 real(kind=kreal) :: dstrain(6) 262 real(kind=kreal) :: alpo(3) 263 logical :: ierr, matlaniso 264 real(kind=kreal) :: stiff_ad(9, nn*3), stiff_aa(9, 9), qf_a(9) 265 real(kind=kreal) :: xj(9, 9) 266 real(kind=kreal) :: tmpforce(9), tmpdisp(9), tmp_d(nn*3), tmp_a(9) 267 real(kind=kreal) :: jacobian(3, 3), inverse(3, 3), inverse1(3, 3), inverse0(3, 3) 268 269 !--------------------------------------------------------------------- 270 271 totaldisp(:, 1:nn) = u(:, :) + du(:, :) - ddu(:, :) 272 totaldisp(:, nn+1:nn+3) = aux(:, :) 273 274 ! we suppose the same material type in the element 275 flag = gausses(1)%pMaterial%nlgeom_flag 276 elem(:, :) = ecoord(:, :) 277 if( flag == UPDATELAG ) elem(:, :) = ecoord(:, :)+totaldisp(:, 1:nn) 278 279 ! --- Inverse of Jacobian at elemental center 280 naturalcoord(:) = 0.0D0 281 call getJacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse) 282 inverse(:, :)= inverse(:, :)*det 283 ! ---- We now calculate stiff matrix include imcompatible mode 284 ! [ Kdd Kda ] 285 ! [ Kad Kaa ] 286 stiff_ad(:, :) = 0.0D0 287 stiff_aa(:, :) = 0.0D0 288 B1(1:6, 1:(nn+3)*ndof) = 0.0D0 289 BN(1:9, 1:(nn+3)*ndof) = 0.0D0 290 291 do LX = 1, NumOfQuadPoints(etype) 292 293 call getQuadPoint(etype, LX, naturalCoord) 294 call getGlobalDeriv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) ) 295 296 if( cdsys_ID > 0 ) then 297 call set_localcoordsys( coords, g_LocalCoordSys(cdsys_ID), coordsys(:, :), serr ) 298 if( serr == -1 ) stop "Fail to setup local coordinate" 299 if( serr == -2 ) then 300 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically" 301 end if 302 end if 303 304 if( present(tt) ) then 305 call getShapeFunc( etype, naturalcoord, spfunc ) 306 ttc = dot_product( tt, spfunc ) 307 call MatlMatrix( gausses(LX), D3, D, time, tincr, coordsys, ttc ) 308 else 309 call MatlMatrix( gausses(LX), D3, D, time, tincr, coordsys ) 310 end if 311 312 if( flag == UPDATELAG ) then 313 call GEOMAT_C3( gausses(LX)%stress, mat ) 314 D(:, :) = D(:, :)-mat 315 endif 316 317 ! -- Derivative of shape function of imcompatible mode -- 318 ! [ -2*a 0, 0 ] 319 ! [ 0, -2*b, 0 ] 320 ! [ 0, 0, -2*c ] 321 ! we don't call getShapeDeriv but use above value directly for computation efficiency 322 gderiv(nn+1, :) = -2.0D0*naturalcoord(1)*inverse(1, :)/det 323 gderiv(nn+2, :) = -2.0D0*naturalcoord(2)*inverse(2, :)/det 324 gderiv(nn+3, :) = -2.0D0*naturalcoord(3)*inverse(3, :)/det 325 326 wg = getWeight(etype, LX)*det 327 B(1:6, 1:(nn+3)*ndof)=0.0D0 328 do j = 1, nn+3 329 B(1, 3*j-2) = gderiv(j, 1) 330 B(2, 3*j-1) = gderiv(j, 2) 331 B(3, 3*j ) = gderiv(j, 3) 332 B(4, 3*j-2) = gderiv(j, 2) 333 B(4, 3*j-1) = gderiv(j, 1) 334 B(5, 3*j-1) = gderiv(j, 3) 335 B(5, 3*j ) = gderiv(j, 2) 336 B(6, 3*j-2) = gderiv(j, 3) 337 B(6, 3*j ) = gderiv(j, 1) 338 end do 339 340 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD ) 341 if( flag == TOTALLAG ) then 342 ! ---dudx(i,j) ==> gdispderiv(i,j) 343 gdispderiv(1:ndof,1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) ) 344 do j = 1, nn+3 345 346 B1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1) 347 B1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1) 348 B1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1) 349 B1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2) 350 B1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2) 351 B1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2) 352 B1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3) 353 B1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3) 354 B1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3) 355 B1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2) 356 B1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2) 357 B1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2) 358 B1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2) 359 B1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2) 360 B1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2) 361 B1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3) 362 B1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3) 363 B1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3) 364 365 end do 366 ! ---BL = BL0 + BL1 367 do j = 1, (nn+3)*ndof 368 B(:, j) = B(:, j)+B1(:, j) 369 end do 370 371 end if 372 373 DB(1:6, 1:(nn+3)*ndof) = matmul( D, B(1:6, 1:(nn+3)*ndof) ) 374 forall( i=1:3*ndof, j=1:nn*ndof ) 375 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( B(:, nn*ndof+i), DB(:, j) )*wg 376 end forall 377 forall( i=1:3*ndof, j=1:3*ndof ) 378 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( B(:, nn*ndof+i), DB(:, nn*ndof+j) )*wg 379 end forall 380 381 ! calculate the stress matrix ( TOTAL LAGRANGE METHOD ) 382 if( flag == TOTALLAG .OR. flag == UPDATELAG ) then 383 stress(1:6) = gausses(LX)%stress 384 do j = 1, nn+3 385 BN(1, 3*j-2) = gderiv(j, 1) 386 BN(2, 3*j-1) = gderiv(j, 1) 387 BN(3, 3*j ) = gderiv(j, 1) 388 BN(4, 3*j-2) = gderiv(j, 2) 389 BN(5, 3*j-1) = gderiv(j, 2) 390 BN(6, 3*j ) = gderiv(j, 2) 391 BN(7, 3*j-2) = gderiv(j, 3) 392 BN(8, 3*j-1) = gderiv(j, 3) 393 BN(9, 3*j ) = gderiv(j, 3) 394 end do 395 Smat(:, :) = 0.0D0 396 do j = 1, 3 397 Smat(j , j ) = stress(1) 398 Smat(j , j+3) = stress(4) 399 Smat(j , j+6) = stress(6) 400 Smat(j+3, j ) = stress(4) 401 Smat(j+3, j+3) = stress(2) 402 Smat(j+3, j+6) = stress(5) 403 Smat(j+6, j ) = stress(6) 404 Smat(j+6, j+3) = stress(5) 405 Smat(j+6, j+6) = stress(3) 406 end do 407 SBN(1:9, 1:(nn+3)*ndof) = matmul( Smat(1:9, 1:9), BN(1:9, 1:(nn+3)*ndof) ) 408 forall( i=1:3*ndof, j=1:nn*ndof ) 409 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( BN(:, nn*ndof+i), SBN(:, j) )*wg 410 end forall 411 forall( i=1:3*ndof, j=1:3*ndof ) 412 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( BN(:, nn*ndof+i), SBN(:, nn*ndof+j) )*wg 413 end forall 414 end if 415 416 end do 417 418 ! -----recover incompatible dof of nodal displacement 419 xj(1:3*ndof, 1:3*ndof)= stiff_aa(1:3*ndof, 1:3*ndof) 420 call calInverse(3*ndof, xj) 421 ! --- [Kad] * ddu 422 do i=1,nn 423 tmp_d((i-1)*ndof+1:i*ndof) = ddu(1:ndof, i) 424 enddo 425 tmpforce(1:3*ndof) = matmul( stiff_ad(1:3*ndof, 1:nn*ndof), tmp_d(1:nn*ndof) ) 426 ! --- ddaux = -[Kaa]-1 * ([Kad] * ddu) 427 tmp_a(1:3*ndof) = -matmul( xj(1:3*ndof, 1:3*ndof), tmpforce(1:3*ndof) ) 428 do i=1,3 429 ddaux(1:ndof, i) = tmp_a((i-1)*ndof+1:i*ndof) 430 enddo 431 432 433 !--------------------------------------------------------------------- 434 435 totaldisp(:,1:nn) = u(:,:)+ du(:,:) 436 totaldisp(:,nn+1:nn+3) = aux(:,:) + ddaux(:,:) 437 438 !--------------------------------------------------------------------- 439 440 qf(:) = 0.0D0 441 qf_a(:) = 0.0D0 442 443 stiff_ad(:, :) = 0.0D0 444 stiff_aa(:, :) = 0.0D0 445 446 !--------------------------------------------------------------------- 447 448 if( flag == UPDATELAG ) then 449 elem(:,:) = (0.5D0*du(:,:)+u(:,:) ) +ecoord(:,:) 450 elem1(:,:) = (du(:,:)+u(:,:) ) +ecoord(:,:) 451 ! elem = elem1 452 totaldisp(:,1:nn) = du(:,:) 453 if( iter == 1 ) aux(:,:) = 0.0D0 !--- is this correct??? 454 totaldisp(:,nn+1:nn+3) = aux(:,:) 455 end if 456 457 matlaniso = .FALSE. 458 if( present(TT) .and. cdsys_ID > 0 ) then 459 ina = TT(1) 460 call fetch_TableData( MC_ORTHOEXP, gausses(1)%pMaterial%dict, alpo(:), ierr, ina ) 461 if( .not. ierr ) matlaniso = .true. 462 end if 463 464 ! --- Inverse of Jacobian at elemental center 465 naturalcoord(:) = 0.0D0 466 call getJacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse) 467 inverse(:, :)= inverse(:, :)*det 468 if( flag == UPDATELAG ) then 469 call getJacobian(etype, nn, naturalcoord, elem1, det, jacobian, inverse1) 470 inverse1(:, :)= inverse1(:, :)*det 471 call getJacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse0) 472 inverse0(:, :)= inverse0(:, :)*det !for deformation gradient F 473 endif 474 475 do LX = 1, NumOfQuadPoints(etype) 476 477 mtype = gausses(LX)%pMaterial%mtype 478 479 call getQuadPoint( etype, LX, naturalCoord(:) ) 480 call getGlobalDeriv(etype, nn, naturalcoord, elem, det, gderiv) 481 482 if( cdsys_ID > 0 ) then 483 call set_localcoordsys( coords, g_LocalCoordSys(cdsys_ID), coordsys(:,:), serr ) 484 if( serr == -1 ) stop "Fail to setup local coordinate" 485 if( serr == -2 ) then 486 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically" 487 end if 488 end if 489 490 ! ======================================================== 491 ! UPDATE STRAIN and STRESS 492 ! ======================================================== 493 494 ! Thermal Strain 495 EPSTH = 0.0D0 496 if( present(tt) .AND. present(t0) ) then 497 call getShapeFunc(etype, naturalcoord, spfunc) 498 ttc = dot_product(TT, spfunc) 499 tt0 = dot_product(T0, spfunc) 500 ttn = dot_product(TN, spfunc) 501 call Cal_Thermal_expansion_C3( tt0, ttc, gausses(LX)%pMaterial, coordsys, matlaniso, EPSTH ) 502 end if 503 504 ! -- Derivative of shape function of imcompatible mode -- 505 ! [ -2*a 0, 0 ] 506 ! [ 0, -2*b, 0 ] 507 ! [ 0, 0, -2*c ] 508 ! we don't call getShapeDeriv but use above value directly for computation efficiency 509 gderiv(nn+1, :) = -2.0D0*naturalcoord(1)*inverse(1, :)/det 510 gderiv(nn+2, :) = -2.0D0*naturalcoord(2)*inverse(2, :)/det 511 gderiv(nn+3, :) = -2.0D0*naturalcoord(3)*inverse(3, :)/det 512 513 ! Small strain 514 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) ) 515 dstrain(1) = gdispderiv(1, 1) 516 dstrain(2) = gdispderiv(2, 2) 517 dstrain(3) = gdispderiv(3, 3) 518 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) ) 519 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) ) 520 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) ) 521 dstrain(:) = dstrain(:)-EPSTH(:) ! allright? 522 523 F(1:3,1:3) = 0.d0; F(1,1)=1.d0; F(2,2)=1.d0; F(3,3)=1.d0; !deformation gradient 524 if( flag == INFINITESIMAL ) then 525 gausses(LX)%strain(1:6) = dstrain(1:6)+EPSTH(:) 526 F(1:3,1:3) = F(1:3,1:3) + gdispderiv(1:3,1:3) 527 528 else if( flag == TOTALLAG ) then 529 ! Green-Lagrange strain 530 dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) ) 531 dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) ) 532 dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) ) 533 dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) & 534 +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) ) 535 dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) & 536 +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) ) 537 dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) & 538 +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) ) 539 540 gausses(LX)%strain(1:6) = dstrain(1:6)+EPSTH(:) 541 542 else if( flag == UPDATELAG ) then 543 ! CALL GEOMAT_C3( gausses(LX)%stress, mat ) 544 ! D(:, :) = D(:, :)+mat(:, :) 545 rot = 0.0D0 546 rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2) 547 rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3) 548 rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3) 549 550 gausses(LX)%strain(1:6) = gausses(LX)%strain_bak(1:6)+ dstrain(1:6)+EPSTH(:) 551 call getGlobalDeriv(etype, nn, naturalcoord, ecoord, det0, gderiv0) 552 gderiv0(nn+1, :) = -2.0D0*naturalcoord(1)*inverse0(1, :)/det0 553 gderiv0(nn+2, :) = -2.0D0*naturalcoord(2)*inverse0(2, :)/det0 554 gderiv0(nn+3, :) = -2.0D0*naturalcoord(3)*inverse0(3, :)/det0 555 F(1:3,1:3) = F(1:3,1:3) + matmul( u(1:ndof, 1:nn+3)+du(1:ndof, 1:nn+3), gderiv0(1:nn+3, 1:ndof) ) 556 557 end if 558 559 ! Update stress 560 if( present(tt) .AND. present(t0) ) then 561 call Update_Stress3D( flag, gausses(LX), rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn ) 562 else 563 call Update_Stress3D( flag, gausses(LX), rot, dstrain, F, coordsys, time, tincr ) 564 end if 565 566 ! ======================================================== 567 ! calculate the internal force ( equivalent nodal force ) 568 ! ======================================================== 569 ! Small strain 570 B(1:6, 1:(nn+3)*ndof) = 0.0D0 571 do J=1,nn+3 572 B(1,3*j-2) = gderiv(j, 1) 573 B(2,3*j-1) = gderiv(j, 2) 574 B(3,3*j ) = gderiv(j, 3) 575 B(4,3*j-2) = gderiv(j, 2) 576 B(4,3*j-1) = gderiv(j, 1) 577 B(5,3*j-1) = gderiv(j, 3) 578 B(5,3*j ) = gderiv(j, 2) 579 B(6,3*j-2) = gderiv(j, 3) 580 B(6,3*j ) = gderiv(j, 1) 581 end do 582 583 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD ) 584 if( flag == INFINITESIMAL ) then 585 586 else if( flag == TOTALLAG ) then 587 588 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) ) 589 B1(1:6, 1:(nn+3)*ndof)=0.0D0 590 do j = 1,nn+3 591 B1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1) 592 B1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1) 593 B1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1) 594 B1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2) 595 B1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2) 596 B1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2) 597 B1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3) 598 B1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3) 599 B1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3) 600 B1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2) 601 B1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2) 602 B1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2) 603 B1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2) 604 B1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2) 605 B1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2) 606 B1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3) 607 B1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3) 608 B1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3) 609 end do 610 ! BL = BL0 + BL1 611 do j=1,(nn+3)*ndof 612 B(:,j) = B(:,j)+B1(:,j) 613 end do 614 615 else if( flag == UPDATELAG ) then 616 617 call getGlobalDeriv(etype, nn, naturalcoord, elem1, det, gderiv(1:nn,1:3)) 618 619 ! -- Derivative of shape function of imcompatible mode -- 620 ! [ -2*a 0, 0 ] 621 ! [ 0, -2*b, 0 ] 622 ! [ 0, 0, -2*c ] 623 ! we don't call getShapeDeriv but use above value directly for computation efficiency 624 gderiv(nn+1, :) = -2.0D0*naturalcoord(1)*inverse1(1, :)/det 625 gderiv(nn+2, :) = -2.0D0*naturalcoord(2)*inverse1(2, :)/det 626 gderiv(nn+3, :) = -2.0D0*naturalcoord(3)*inverse1(3, :)/det 627 628 B(1:6, 1:(nn+3)*ndof) = 0.0D0 629 do j = 1, nn+3 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 end if 642 643 WG=getWeight( etype, LX )*DET 644 645 DB(1:6, 1:(nn+3)*ndof) = matmul( D, B(1:6, 1:(nn+3)*ndof) ) 646 forall( i=1:3*ndof, j=1:nn*ndof ) 647 stiff_ad(i, j) = stiff_ad(i, j)+dot_product( B(:, nn*ndof+i), DB(:, j) )*wg 648 end forall 649 forall( i=1:3*ndof, j=1:3*ndof ) 650 stiff_aa(i, j) = stiff_aa(i, j)+dot_product( B(:, nn*ndof+i), DB(:, nn*ndof+j) )*wg 651 end forall 652 653 ! calculate the Internal Force 654 qf(1:nn*ndof) & 655 = qf(1:nn*ndof)+matmul( gausses(LX)%stress(1:6), B(1:6,1:nn*ndof) )*WG 656 qf_a(1:3*ndof) & 657 = qf_a(1:3*ndof)+matmul( gausses(LX)%stress(1:6), B(1:6,nn*ndof+1:(nn+3)*ndof) )*WG 658 end do 659 660 ! condence ( qf - [Kda] * [Kaa]-1 * qf_a ) 661 xj(1:9, 1:9)= stiff_aa(1:3*ndof, 1:3*ndof) 662 call calInverse(9, xj) 663 tmpdisp(:) = matmul( xj(:,:), qf_a(:) ) 664 do i=1,nn*ndof 665 qf(i) = qf(i) - dot_product( stiff_ad(:,i), tmpdisp(:) ) 666 enddo 667 end subroutine UPDATE_C3D8IC 668 669 670 !> This subroutine calculatess thermal loading 671 !----------------------------------------------------------------------* 672 subroutine TLOAD_C3D8IC & 673 (etype, nn, xx, yy, zz, tt, t0, & 674 gausses, VECT, cdsys_ID, coords) 675 !----------------------------------------------------------------------* 676 677 use m_fstr 678 use mMechGauss 679 use m_MatMatrix 680 use m_common_struct 681 682 !--------------------------------------------------------------------- 683 684 integer(kind=kint), parameter :: ndof=3 685 integer(kind=kint), intent(in) :: etype !< element type, not used here 686 integer(kind=kint), intent(in) :: nn !< number of element nodes 687 real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn) !< nodes coordinate of element 688 real(kind=kreal), intent(in) :: tt(nn), t0(nn) !< current and ref temprature 689 type(tGaussStatus), intent(inout) :: gausses(:) !< info about qudrature points 690 real(kind=kreal), intent(out) :: VECT(nn*ndof) !< load vector 691 integer(kind=kint), intent(in) :: cdsys_ID 692 real(kind=kreal), intent(inout) :: coords(3, 3) !< variables to define matreial coordinate system 693 694 !--------------------------------------------------------------------- 695 696 real(kind=kreal) :: ALP, ALP0 697 real(kind=kreal) :: D(6, 6), B(6, ndof*(nn+3)), DB_a(6, ndof*3) 698 real(kind=kreal) :: det, wg, ecoord(3, nn) 699 integer(kind=kint) :: i, j, IC, serr 700 real(kind=kreal) :: EPSTH(6), SGM(6), H(nn), alpo(3), alpo0(3), coordsys(3, 3), xj(9, 9) 701 real(kind=kreal) :: naturalcoord(3), gderiv(nn+3, 3) 702 real(kind=kreal) :: jacobian(3, 3),inverse(3, 3) 703 real(kind=kreal) :: stiff_da(nn*3, 3*3), stiff_aa(3*3, 3*3) 704 real(kind=kreal) :: VECT_a(3*3) 705 real(kind=kreal) :: icdisp(9) 706 real(kind=kreal) :: TEMPC, TEMP0, THERMAL_EPS, tm(6, 6), outa(1), ina(1) 707 logical :: ierr, matlaniso 708 709 !--------------------------------------------------------------------- 710 711 matlaniso = .FALSE. 712 if( cdsys_ID > 0 ) then ! cannot define aniso exapansion when no local coord defined 713 ina = TT(1) 714 call fetch_TableData( MC_ORTHOEXP, gausses(1)%pMaterial%dict, alpo(:), ierr, ina ) 715 if( .not. ierr ) matlaniso = .true. 716 end if 717 718 ecoord(1, :) = XX(:) 719 ecoord(2, :) = YY(:) 720 ecoord(3, :) = ZZ(:) 721 ! ---- Calculate enhanced displacement at first 722 ! --- Inverse of Jacobian at elemental center 723 naturalcoord(:) = 0.0D0 724 call getJacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse) 725 inverse(:, :) = inverse(:, :)*det 726 ! ---- We now calculate stiff matrix include imcompatible mode 727 ! [ Kdd Kda ] 728 ! [ Kad Kaa ] 729 stiff_da(:, :) = 0.0D0 730 stiff_aa(:, :) = 0.0D0 731 B(1:6, 1:(nn+3)*ndof) = 0.0D0 732 VECT(1:nn*ndof) = 0.0D0 733 VECT_a(1:3*ndof) = 0.0D0 734 735 do IC = 1, NumOfQuadPoints(etype) 736 737 call getQuadPoint(etype, IC, naturalCoord) 738 call getGlobalDeriv( etype, nn, naturalcoord, ecoord, det, gderiv(1:nn, 1:3) ) 739 740 if( cdsys_ID > 0 ) then 741 call set_localcoordsys( coords, g_LocalCoordSys(cdsys_ID), coordsys(:, :), serr ) 742 if( serr == -1 ) stop "Fail to setup local coordinate" 743 if( serr == -2 ) then 744 write(*,*) "WARNING! Cannot setup local coordinate, it is modified automatically" 745 end if 746 end if 747 748 call getShapeFunc( etype, naturalcoord, H(1:nn) ) 749 TEMPC = dot_product( H(1:nn), TT(1:nn) ) 750 TEMP0 = dot_product( H(1:nn), T0(1:nn) ) 751 call MatlMatrix( gausses(IC), D3, D, 1.d0, 1.0D0, coordsys, TEMPC ) 752 753 ! -- Derivative of shape function of imcompatible mode -- 754 ! [ -2*a 0, 0 ] 755 ! [ 0, -2*b, 0 ] 756 ! [ 0, 0, -2*c ] 757 ! we don't call getShapeDeriv but use above value directly for computation efficiency 758 gderiv(nn+1, :) = -2.0D0*naturalcoord(1)*inverse(1, :)/det 759 gderiv(nn+2, :) = -2.0D0*naturalcoord(2)*inverse(2, :)/det 760 gderiv(nn+3, :) = -2.0D0*naturalcoord(3)*inverse(3, :)/det 761 762 wg = getWeight(etype, IC)*det 763 do j = 1, nn+3 764 B(1, 3*j-2) = gderiv(j, 1) 765 B(2, 3*j-1) = gderiv(j, 2) 766 B(3, 3*j ) = gderiv(j, 3) 767 B(4, 3*j-2) = gderiv(j, 2) 768 B(4, 3*j-1) = gderiv(j, 1) 769 B(5, 3*j-1) = gderiv(j, 3) 770 B(5, 3*j ) = gderiv(j, 2) 771 B(6, 3*j-2) = gderiv(j, 3) 772 B(6, 3*j ) = gderiv(j, 1) 773 end do 774 775 DB_a(1:6, 1:3*ndof) = matmul( D, B(1:6, nn*ndof+1:(nn+3)*ndof) ) 776 forall( i=1:nn*ndof, j=1:3*ndof ) 777 stiff_da(i, j) = stiff_da(i, j) + dot_product( B(1:6, i), DB_a(1:6, j) ) * wg 778 end forall 779 forall( i=1:3*ndof, j=1:3*ndof ) 780 stiff_aa(i, j) = stiff_aa(i, j) + dot_product( B(1:6, nn*ndof+i), DB_a(1:6, j) ) * wg 781 end forall 782 783 ina(1) = TEMPC 784 if( matlaniso ) then 785 call fetch_TableData( MC_ORTHOEXP, gausses(IC)%pMaterial%dict, alpo(:), ierr, ina ) 786 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!" 787 else 788 call fetch_TableData( MC_THEMOEXP, gausses(IC)%pMaterial%dict, outa(:), ierr, ina ) 789 if( ierr ) outa(1) = gausses(IC)%pMaterial%variables(M_EXAPNSION) 790 alp = outa(1) 791 end if 792 ina(1) = TEMP0 793 if( matlaniso ) then 794 call fetch_TableData( MC_ORTHOEXP, gausses(IC)%pMaterial%dict, alpo0(:), ierr, ina ) 795 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!" 796 else 797 call fetch_TableData( MC_THEMOEXP, gausses(IC)%pMaterial%dict, outa(:), ierr, ina ) 798 if( ierr ) outa(1) = gausses(IC)%pMaterial%variables(M_EXAPNSION) 799 alp0 = outa(1) 800 end if 801 802 !** 803 !** THERMAL strain 804 !** 805 if( matlaniso ) then 806 do j = 1, 3 807 EPSTH(j) = ALPO(j)*(TEMPC-ref_temp)-alpo0(j)*(TEMP0-ref_temp) 808 end do 809 EPSTH(4:6) = 0.0D0 810 call transformation(coordsys, tm) 811 EPSTH(:) = matmul( EPSTH(:), tm ) ! to global coord 812 EPSTH(4:6) = EPSTH(4:6)*2.0D0 813 else 814 THERMAL_EPS = ALP*(TEMPC-ref_temp)-alp0*(TEMP0-ref_temp) 815 EPSTH(1:3) = THERMAL_EPS 816 EPSTH(4:6) = 0.0D0 817 end if 818 819 !** 820 !** SET SGM {s}=[D]{e} 821 !** 822 SGM(:) = matmul( D(:, :), EPSTH(:) ) 823 824 !** 825 !** CALCULATE LOAD {F}=[B]T{e} 826 !** 827 VECT(1:nn*ndof) = VECT(1:nn*ndof)+matmul( SGM(1:6), B(1:6, 1:nn*ndof) )*wg 828 VECT_a(1:3*ndof) = VECT_a(1:3*ndof)+matmul( SGM(1:6), B(1:6, nn*ndof+1:(nn+3)*ndof) )*wg 829 830 end do 831 832 ! --- condense vect 833 xj(1:9,1:9)= stiff_aa(1:9, 1:9) 834 call calInverse(9, xj) 835 ! --- [Kaa]-1 * fa 836 icdisp(1:9) = matmul( xj(:, :), VECT_a(1:3*ndof) ) 837 ! --- fd - [Kda] * [Kaa]-1 * fa 838 VECT(1:nn*ndof) = VECT(1:nn*ndof) - matmul( stiff_da(1:nn*ndof, 1:3*ndof), icdisp(1:3*ndof) ) 839 840 end subroutine TLOAD_C3D8IC 841 842 843end module m_static_LIB_3dIC 844