1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5 6module m_static_LIB_3d_vp 7 8 use hecmw, only : kint, kreal 9 use elementInfo 10 11 implicit none 12 13contains 14 !-------------------------------------------------------------------- 15 subroutine STF_C3_vp & 16 (etype, nn, ecoord, gausses, stiff, tincr, & 17 v, temperature) 18 !-------------------------------------------------------------------- 19 20 use mMechGauss 21 22 !-------------------------------------------------------------------- 23 24 integer(kind=kint), intent(in) :: etype !< element type 25 integer(kind=kint), intent(in) :: nn !< the number of elemental nodes 26 real(kind=kreal), intent(in) :: ecoord(3, nn) !< coordinates of elemental nodes 27 type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points 28 real(kind=kreal), intent(out) :: stiff(:,:) !< stiff matrix 29 real(kind=kreal), intent(in) :: tincr !< time increment 30 real(kind=kreal), intent(in), optional :: v(:, :) !< nodal velocity 31 real(kind=kreal), intent(in), optional :: temperature(nn) !< temperature 32 33 !-------------------------------------------------------------------- 34 35 integer(kind=kint) :: i, j 36 integer(kind=kint) :: na, nb 37 integer(kind=kint) :: isize, jsize 38 integer(kind=kint) :: LX 39 40 real(kind=kreal) :: MM(nn, nn), AA(nn, nn), DD(nn, nn, 3, 3), & 41 trD(nn, nn), BB(nn, nn), CC(nn, nn, 3), & 42 MS(nn, nn), AS(nn, nn), CS(nn, nn, 3), & 43 MP(nn, nn, 3), AP(nn, nn, 3), CP(nn, nn) 44 real(kind=kreal) :: spfunc(nn), gderiv(nn, 3) 45 real(kind=kreal) :: elem(3, nn) 46 real(kind=kreal) :: naturalCoord(3) 47 real(kind=kreal) :: dndx(nn, 3) 48 real(kind=kreal) :: tincr_inv 49 real(kind=kreal) :: volume, volume_inv 50 real(kind=kreal) :: mu 51 real(kind=kreal) :: rho, rho_inv 52 real(kind=kreal) :: vx, vy, vz 53 real(kind=kreal) :: t1, t2, t3 54 real(kind=kreal) :: v_dot_v 55 real(kind=kreal) :: d 56 real(kind=kreal) :: det, wg 57 real(kind=kreal) :: tau 58 real(kind=kreal), parameter :: gamma = 0.5D0 59 60 !-------------------------------------------------------------------- 61 62 tincr_inv = 1.0D0/tincr 63 64 !-------------------------------------------------------------------- 65 66 elem(:, :) = ecoord(:, :) 67 68 !-------------------------------------------------------------------- 69 70 t1 = 2.0D0*tincr_inv 71 72 !--------------------------------------------------------------- 73 74 volume = 0.0D0 75 76 loopVolume: do LX = 1, NumOfQuadPoints(etype) 77 78 !---------------------------------------------------------- 79 80 call getQuadPoint( etype, LX, naturalCoord(:) ) 81 call getShapeFunc(etype, naturalCoord, spfunc) 82 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 83 84 !---------------------------------------------------------- 85 86 wg = getWeight(etype, LX)*det 87 88 !---------------------------------------------------------- 89 90 volume = volume+wg 91 92 !---------------------------------------------------------- 93 94 end do loopVolume 95 96 volume_inv = 1.0D0/volume 97 98 !--------------------------------------------------------------- 99 100 naturalCoord(1) = 0.25D0 101 naturalCoord(2) = 0.25D0 102 naturalCoord(3) = 0.25D0 103 104 call getShapeFunc(etype, naturalCoord, spfunc) 105 106 vx = 0.0D0 107 vy = 0.0D0 108 vz = 0.0D0 109 110 do na = 1, nn 111 112 vx = vx+spfunc(na)*v(1, na) 113 vy = vy+spfunc(na)*v(2, na) 114 vz = vz+spfunc(na)*v(3, na) 115 116 end do 117 118 v_dot_v = vx*vx+vy*vy+vz*vz 119 120 !--------------------------------------------------------------- 121 122 mu = 0.0D0 123 rho = 0.0D0 124 125 do na = 1, nn 126 127 dndx(na, 1) = 0.0D0 128 dndx(na, 2) = 0.0D0 129 dndx(na, 3) = 0.0D0 130 131 end do 132 133 loopGlobalDeriv: do LX = 1, NumOfQuadPoints(etype) 134 135 !---------------------------------------------------------- 136 137 call getQuadPoint( etype, LX, naturalCoord(1:3) ) 138 call getShapeFunc(etype, naturalCoord, spfunc) 139 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 140 141 !---------------------------------------------------------- 142 143 wg = getWeight(etype, LX)*det 144 145 !---------------------------------------------------------- 146 147 mu = mu+wg*gausses(LX)%pMaterial%variables(M_VISCOCITY) 148 149 rho = rho+wg*gausses(LX)%pMaterial%variables(M_DENSITY) 150 151 !---------------------------------------------------------- 152 153 do na = 1, nn 154 155 dndx(na, 1) = dndx(na, 1)+wg*gderiv(na, 1) 156 dndx(na, 2) = dndx(na, 2)+wg*gderiv(na, 2) 157 dndx(na, 3) = dndx(na, 3)+wg*gderiv(na, 3) 158 159 end do 160 161 !---------------------------------------------------------- 162 163 end do loopGlobalDeriv 164 165 mu = volume_inv*mu 166 rho = volume_inv*rho 167 168 do na = 1, nn 169 170 dndx(na, 1) = volume_inv*dndx(na, 1) 171 dndx(na, 2) = volume_inv*dndx(na, 2) 172 dndx(na, 3) = volume_inv*dndx(na, 3) 173 174 end do 175 176 !--------------------------------------------------------------- 177 178 d = 0.0D0 179 180 do na = 1, nn 181 182 d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) ) 183 184 end do 185 186 ! h_es3d = 2.0D0/( d/DSQRT( v_dot_v ) ) 187 188 !--------------------------------------------------------------- 189 190 ! t2 = 2.0D0*DSQRT( v_dot_v )/h_es3d 191 t2 = d 192 193 !---------------------------------------------------------------- 194 195 if( v_dot_v .LT. 1.0D-15 ) then 196 197 t3 = 4.0D0*mu/( rho*volume**(2.0D0/3.0D0) ) 198 199 else 200 201 t3 = mu*d*d/( rho*v_dot_v ) 202 203 end if 204 205 !---------------------------------------------------------------- 206 207 tau = 1.0D0/dsqrt( t1*t1+t2*t2+t3*t3 ) 208 209 !-------------------------------------------------------------------- 210 211 stiff(:, :) = 0.0D0 212 213 loopGauss: do LX = 1, NumOfQuadPoints(etype) 214 215 !---------------------------------------------------------- 216 217 mu = gausses(LX)%pMaterial%variables(M_VISCOCITY) 218 219 rho = gausses(LX)%pMaterial%variables(M_DENSITY) 220 rho_inv = 1.0D0/rho 221 222 !---------------------------------------------------------- 223 224 call getQuadPoint( etype, LX, naturalCoord(1:3) ) 225 call getShapeFunc( etype, naturalCoord(:), spfunc(:) ) 226 call getGlobalDeriv( etype, nn, naturalCoord(:), elem(:,:), & 227 det, gderiv(:,:) ) 228 229 !---------------------------------------------------------- 230 231 wg = getWeight(etype, LX)*det 232 233 !---------------------------------------------------------- 234 235 vx = 0.0D0 236 vy = 0.0D0 237 vz = 0.0D0 238 239 do na = 1, nn 240 241 vx = vx+spfunc(na)*v(1, na) 242 vy = vy+spfunc(na)*v(2, na) 243 vz = vz+spfunc(na)*v(3, na) 244 245 end do 246 247 !---------------------------------------------------------- 248 249 forall(na = 1:nn, nb = 1:nn) 250 251 MM(na, nb) = spfunc(na)*spfunc(nb) 252 AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) & 253 +vy*( spfunc(na)*gderiv(nb, 2) ) & 254 +vz*( spfunc(na)*gderiv(nb, 3) ) 255 DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1) 256 DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2) 257 DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3) 258 DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1) 259 DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2) 260 DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3) 261 DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1) 262 DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2) 263 DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3) 264 trD(na, nb) = DD(na, nb, 1, 1) & 265 +DD(na, nb, 2, 2) & 266 +DD(na, nb, 3, 3) 267 BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) & 268 +( vx*vy )*DD(na, nb, 1, 2) & 269 +( vx*vz )*DD(na, nb, 1, 3) & 270 +( vy*vx )*DD(na, nb, 2, 1) & 271 +( vy*vy )*DD(na, nb, 2, 2) & 272 +( vy*vz )*DD(na, nb, 2, 3) & 273 +( vz*vx )*DD(na, nb, 3, 1) & 274 +( vz*vy )*DD(na, nb, 3, 2) & 275 +( vz*vz )*DD(na, nb, 3, 3) 276 CC(na, nb, 1) = gderiv(na, 1)*spfunc(nb) 277 CC(na, nb, 2) = gderiv(na, 2)*spfunc(nb) 278 CC(na, nb, 3) = gderiv(na, 3)*spfunc(nb) 279 280 MS(nb, na) = AA(na, nb) 281 AS(na, nb) = BB(na, nb) 282 CS(na, nb, 1) = vx*DD(na, nb, 1, 1) & 283 +vy*DD(na, nb, 2, 1) & 284 +vz*DD(na, nb, 3, 1) 285 CS(na, nb, 2) = vx*DD(na, nb, 1, 2) & 286 +vy*DD(na, nb, 2, 2) & 287 +vz*DD(na, nb, 3, 2) 288 CS(na, nb, 3) = vx*DD(na, nb, 1, 3) & 289 +vy*DD(na, nb, 2, 3) & 290 +vz*DD(na, nb, 3, 3) 291 MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1) 292 MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2) 293 MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3) 294 AP(nb, na, 1) = CS(na, nb, 1) 295 AP(nb, na, 2) = CS(na, nb, 2) 296 AP(nb, na, 3) = CS(na, nb, 3) 297 CP(na, nb) = trD(na, nb) 298 299 end forall 300 301 !---------------------------------------------------------- 302 303 do nb = 1, nn 304 305 do na = 1, nn 306 307 i = 1 308 j = 1 309 isize = 4*(na-1)+i 310 jsize = 4*(nb-1)+j 311 312 stiff(isize, jsize) & 313 = stiff(isize, jsize) & 314 +wg & 315 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 316 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 317 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 318 319 i = 1 320 j = 2 321 isize = 4*(na-1)+i 322 jsize = 4*(nb-1)+j 323 324 stiff(isize, jsize) & 325 = stiff(isize, jsize) & 326 +wg & 327 *( gamma*mu*DD(na, nb, 2, 1) ) 328 329 i = 1 330 j = 3 331 isize = 4*(na-1)+i 332 jsize = 4*(nb-1)+j 333 334 stiff(isize, jsize) & 335 = stiff(isize, jsize) & 336 +wg & 337 *( gamma*mu*DD(na, nb, 3, 1) ) 338 339 i = 1 340 j = 4 341 isize = 4*(na-1)+i 342 jsize = 4*(nb-1)+j 343 344 stiff(isize, jsize) & 345 = stiff(isize, jsize) & 346 +wg & 347 *( -CC(na, nb, 1)+tau*CS(na, nb, 1) ) 348 349 i = 2 350 j = 1 351 isize = 4*(na-1)+i 352 jsize = 4*(nb-1)+j 353 354 stiff(isize, jsize) & 355 = stiff(isize, jsize) & 356 +wg & 357 *( gamma*mu*DD(na, nb, 1, 2) ) 358 359 i = 2 360 j = 2 361 isize = 4*(na-1)+i 362 jsize = 4*(nb-1)+j 363 364 stiff(isize, jsize) & 365 = stiff(isize, jsize) & 366 +wg & 367 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 368 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 369 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 370 371 i = 2 372 j = 3 373 isize = 4*(na-1)+i 374 jsize = 4*(nb-1)+j 375 376 stiff(isize, jsize) & 377 = stiff(isize, jsize) & 378 +wg & 379 *( gamma*mu*DD(na, nb, 3, 2) ) 380 381 i = 2 382 j = 4 383 isize = 4*(na-1)+i 384 jsize = 4*(nb-1)+j 385 386 stiff(isize, jsize) & 387 = stiff(isize, jsize) & 388 +wg & 389 *( -CC(na, nb, 2)+tau*CS(na, nb, 2) ) 390 391 i = 3 392 j = 1 393 isize = 4*(na-1)+i 394 jsize = 4*(nb-1)+j 395 396 stiff(isize, jsize) & 397 = stiff(isize, jsize) & 398 +wg & 399 *( gamma*mu*DD(na, nb, 1, 3) ) 400 401 i = 3 402 j = 2 403 isize = 4*(na-1)+i 404 jsize = 4*(nb-1)+j 405 406 stiff(isize, jsize) & 407 = stiff(isize, jsize) & 408 +wg & 409 *( gamma*mu*DD(na, nb, 2, 3) ) 410 411 i = 3 412 j = 3 413 isize = 4*(na-1)+i 414 jsize = 4*(nb-1)+j 415 416 stiff(isize, jsize) & 417 = stiff(isize, jsize) & 418 +wg & 419 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 420 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 421 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 422 423 i = 3 424 j = 4 425 isize = 4*(na-1)+i 426 jsize = 4*(nb-1)+j 427 428 stiff(isize, jsize) & 429 = stiff(isize, jsize) & 430 +wg & 431 *( -CC(na, nb, 3)+tau*CS(na, nb, 3) ) 432 433 i = 4 434 j = 1 435 isize = 4*(na-1)+i 436 jsize = 4*(nb-1)+j 437 438 stiff(isize, jsize) & 439 = stiff(isize, jsize) & 440 +wg & 441 *( CC(nb, na, j) & 442 +tincr_inv*tau*MP(na, nb, j) & 443 +gamma*tau*AP(na, nb, j) ) 444 445 i = 4 446 j = 2 447 isize = 4*(na-1)+i 448 jsize = 4*(nb-1)+j 449 450 stiff(isize, jsize) & 451 = stiff(isize, jsize) & 452 +wg & 453 *( CC(nb, na, j) & 454 +tincr_inv*tau*MP(na, nb, j) & 455 +gamma*tau*AP(na, nb, j) ) 456 457 i = 4 458 j = 3 459 isize = 4*(na-1)+i 460 jsize = 4*(nb-1)+j 461 462 stiff(isize, jsize) & 463 = stiff(isize, jsize) & 464 +wg & 465 *( CC(nb, na, j) & 466 +tincr_inv*tau*MP(na, nb, j) & 467 +gamma*tau*AP(na, nb, j) ) 468 469 i = 4 470 j = 4 471 isize = 4*(na-1)+i 472 jsize = 4*(nb-1)+j 473 474 stiff(isize, jsize) & 475 = stiff(isize, jsize) & 476 +wg & 477 *( rho_inv*tau*trD(na, nb) ) 478 479 end do 480 481 end do 482 483 !---------------------------------------------------------- 484 485 end do loopGauss 486 487 !-------------------------------------------------------------------- 488 end subroutine STF_C3_vp 489 !-------------------------------------------------------------------- 490 491 492 !-------------------------------------------------------------------- 493 subroutine UPDATE_C3_vp & 494 (etype, nn, ecoord, v, dv, gausses) 495 !-------------------------------------------------------------------- 496 497 use mMechGauss 498 499 !-------------------------------------------------------------------- 500 501 integer(kind=kint), intent(in) :: etype !< element type 502 integer(kind=kint), intent(in) :: nn !< the number of elemental nodes 503 real(kind=kreal), intent(in) :: ecoord(3, nn) !< coordinates of elemental nodes 504 real(kind=kreal), intent(in) :: v(4, nn) !< nodal velcoity 505 real(kind=kreal), intent(in) :: dv(4, nn) !< nodal velocity increment 506 type(tGaussStatus), intent(inout) :: gausses(:) !< status of qudrature points 507 508 !-------------------------------------------------------------------- 509 510 integer(kind=kint) :: LX 511 512 real(kind=kreal) :: elem(3, nn) 513 real(kind=kreal) :: totalvelo(4, nn) 514 real(kind=kreal) :: spfunc(nn), gderiv(nn, 3) 515 real(kind=kreal) :: gveloderiv(3, 3) 516 real(kind=kreal) :: naturalCoord(3) 517 real(kind=kreal) :: det 518 real(kind=kreal) :: mu 519 real(kind=kreal) :: p 520 521 !-------------------------------------------------------------------- 522 523 elem(:, :) = ecoord(:, :) 524 525 totalvelo(:, :) = v(:, :)+dv(:, :) 526 527 !-------------------------------------------------------------------- 528 529 loopMatrix: do LX = 1, NumOfQuadPoints(etype) 530 531 !---------------------------------------------------------- 532 533 mu = gausses(LX)%pMaterial%variables(M_VISCOCITY) 534 535 !---------------------------------------------------------- 536 537 call getQuadPoint( etype, LX, naturalCoord(:) ) 538 call getShapeFunc(etype, naturalCoord, spfunc) 539 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 540 541 !---------------------------------------------------------- 542 543 ! Deformation rate tensor 544 gveloderiv(1:3, 1:3) = matmul( totalvelo(1:3, 1:nn), gderiv(1:nn, 1:3) ) 545 gausses(LX)%strain(1) = gveloderiv(1, 1) 546 gausses(LX)%strain(2) = gveloderiv(2, 2) 547 gausses(LX)%strain(3) = gveloderiv(3, 3) 548 gausses(LX)%strain(4) = 0.5D0*( gveloderiv(1, 2)+gveloderiv(2, 1) ) 549 gausses(LX)%strain(5) = 0.5D0*( gveloderiv(2, 3)+gveloderiv(3, 2) ) 550 gausses(LX)%strain(6) = 0.5D0*( gveloderiv(3, 1)+gveloderiv(1, 3) ) 551 552 !---------------------------------------------------------- 553 554 ! Pressure 555 p = dot_product(totalvelo(4, 1:nn), spfunc(1:nn)) 556 557 ! Cauchy stress tensor 558 gausses(LX)%stress(1) = -p+2.0D0*mu*gausses(LX)%strain(1) 559 gausses(LX)%stress(2) = -p+2.0D0*mu*gausses(LX)%strain(2) 560 gausses(LX)%stress(3) = -p+2.0D0*mu*gausses(LX)%strain(3) 561 gausses(LX)%stress(4) = 2.0D0*mu*gausses(LX)%strain(4) 562 gausses(LX)%stress(5) = 2.0D0*mu*gausses(LX)%strain(5) 563 gausses(LX)%stress(6) = 2.0D0*mu*gausses(LX)%strain(6) 564 565 !---------------------------------------------------------- 566 567 !set stress and strain for output 568 gausses(LX)%strain_out(1:6) = gausses(LX)%strain(1:6) 569 gausses(LX)%stress_out(1:6) = gausses(LX)%stress(1:6) 570 571 end do loopMatrix 572 573 !---------------------------------------------------------------- 574 575 !-------------------------------------------------------------------- 576 end subroutine UPDATE_C3_vp 577 !-------------------------------------------------------------------- 578 579 580 !-------------------------------------------------------------------- 581 subroutine LOAD_C3_vp & 582 (etype, nn, ecoord, v, dv, r, gausses, tincr) 583 !-------------------------------------------------------------------- 584 585 use mMechGauss 586 587 !-------------------------------------------------------------------- 588 589 integer(kind=kint), intent(in) :: etype !< element type 590 integer(kind=kint), intent(in) :: nn !< the number of elemental nodes 591 real(kind=kreal), intent(in) :: ecoord(3, nn) !< coordinates of elemental nodes 592 real(kind=kreal), intent(in) :: v(4, nn) !< nodal dislplacements 593 real(kind=kreal), intent(in) :: dv(4, nn) !< nodal velocity increment 594 real(kind=kreal), intent(out) :: r(4*nn) !< elemental residual 595 type(tGaussStatus), intent(inout) :: gausses(:) !< status of qudrature points 596 real(kind=kreal), intent(in) :: tincr !< time increment 597 598 !-------------------------------------------------------------------- 599 600 integer(kind=kint) :: i, j, k 601 integer(kind=kint) :: na, nb 602 integer(kind=kint) :: isize, jsize 603 integer(kind=kint) :: LX 604 605 real(kind=kreal) :: elem(3, nn) 606 real(kind=kreal) :: velo_new(4*nn) 607 real(kind=kreal) :: stiff(4*nn, 4*nn) 608 real(kind=kreal) :: b(4*nn) 609 real(kind=kreal) :: MM(nn, nn), AA(nn, nn), DD(nn, nn, 3, 3), & 610 trD(nn, nn), BB(nn, nn), CC(nn, nn, 3), & 611 MS(nn, nn), AS(nn, nn), CS(nn, nn, 3), & 612 MP(nn, nn, 3), AP(nn, nn, 3), CP(nn, nn) 613 real(kind=kreal) :: spfunc(nn), gderiv(nn, 3) 614 real(kind=kreal) :: naturalCoord(3) 615 real(kind=kreal) :: dndx(nn, 3) 616 real(kind=kreal) :: tincr_inv 617 real(kind=kreal) :: volume, volume_inv 618 real(kind=kreal) :: mu 619 real(kind=kreal) :: rho, rho_inv 620 real(kind=kreal) :: vx, vy, vz 621 real(kind=kreal) :: t1, t2, t3 622 real(kind=kreal) :: v_a_dot_v_a 623 real(kind=kreal) :: d 624 real(kind=kreal) :: det, wg 625 real(kind=kreal) :: tau 626 real(kind=kreal) :: m_v(3), a_v(3), d_v(3, 3, 3), & 627 ms_v(3), as_v(3), mp_dot_v, ap_dot_v 628 real(kind=kreal) :: stiff_velo 629 real(kind=kreal), parameter :: gamma = 0.5D0 630 631 !-------------------------------------------------------------------- 632 633 tincr_inv = 1.0D0/tincr 634 635 !-------------------------------------------------------------------- 636 637 elem(:, :) = ecoord(:, :) 638 639 forall(na = 1:nn, i = 1:4) 640 641 velo_new( 4*(na-1)+i ) = v(i, na)+dv(i, na) 642 643 end forall 644 645 !-------------------------------------------------------------------- 646 647 t1 = 2.0D0*tincr_inv 648 649 !--------------------------------------------------------------- 650 651 volume = 0.0D0 652 653 loopVolume: do LX = 1, NumOfQuadPoints(etype) 654 655 !---------------------------------------------------------- 656 657 call getQuadPoint( etype, LX, naturalCoord(:) ) 658 call getShapeFunc(etype, naturalCoord, spfunc) 659 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 660 661 !---------------------------------------------------------- 662 663 wg = getWeight(etype, LX)*det 664 665 !---------------------------------------------------------- 666 667 volume = volume+wg 668 669 !---------------------------------------------------------- 670 671 end do loopVolume 672 673 volume_inv = 1.0D0/volume 674 675 !--------------------------------------------------------------- 676 677 naturalCoord(1) = 0.25D0 678 naturalCoord(2) = 0.25D0 679 naturalCoord(3) = 0.25D0 680 681 call getShapeFunc(etype, naturalCoord, spfunc) 682 683 vx = 0.0D0 684 vy = 0.0D0 685 vz = 0.0D0 686 687 do na = 1, nn 688 689 vx = vx+spfunc(na)*v(1, na) 690 vy = vy+spfunc(na)*v(2, na) 691 vz = vz+spfunc(na)*v(3, na) 692 693 end do 694 695 v_a_dot_v_a = vx*vx+vy*vy+vz*vz 696 697 !--------------------------------------------------------------- 698 699 do na = 1, nn 700 701 dndx(na, 1) = 0.0D0 702 dndx(na, 2) = 0.0D0 703 dndx(na, 3) = 0.0D0 704 705 end do 706 707 mu = 0.0d0 708 rho = 0.0d0 709 710 loopGlobalDeriv: do LX = 1, NumOfQuadPoints(etype) 711 712 !---------------------------------------------------------- 713 714 call getQuadPoint( etype, LX, naturalCoord(1:3) ) 715 call getShapeFunc(etype, naturalCoord, spfunc) 716 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 717 718 !---------------------------------------------------------- 719 720 wg = getWeight(etype, LX)*det 721 722 !---------------------------------------------------------- 723 724 mu = mu +wg*gausses(LX)%pMaterial%variables(M_VISCOCITY) 725 rho = rho+wg*gausses(LX)%pMaterial%variables(M_DENSITY) 726 727 !---------------------------------------------------------- 728 729 do na = 1, nn 730 731 dndx(na, 1) = dndx(na, 1)+wg*gderiv(na, 1) 732 dndx(na, 2) = dndx(na, 2)+wg*gderiv(na, 2) 733 dndx(na, 3) = dndx(na, 3)+wg*gderiv(na, 3) 734 735 end do 736 737 !---------------------------------------------------------- 738 739 end do loopGlobalDeriv 740 741 mu = volume_inv*mu 742 rho = volume_inv*rho 743 744 do na = 1, nn 745 746 dndx(na, 1) = volume_inv*dndx(na, 1) 747 dndx(na, 2) = volume_inv*dndx(na, 2) 748 dndx(na, 3) = volume_inv*dndx(na, 3) 749 750 end do 751 752 !--------------------------------------------------------------- 753 754 d = 0.0D0 755 756 do na = 1, nn 757 758 d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) ) 759 760 end do 761 762 ! h_es3d = 2.0D0/( d/DSQRT( v_dot_v ) ) 763 764 !--------------------------------------------------------------- 765 766 ! t2 = 2.0D0*DSQRT( v_dot_v )/h_es3d 767 t2 = d 768 769 !---------------------------------------------------------------- 770 771 if( v_a_dot_v_a .LT. 1.0D-15 ) then 772 773 t3 = 4.0D0*mu/( rho*volume**(2.0D0/3.0D0) ) 774 775 else 776 777 t3 = mu*d*d/( rho*v_a_dot_v_a ) 778 779 end if 780 781 !---------------------------------------------------------------- 782 783 tau = 1.0D0/dsqrt( t1*t1+t2*t2+t3*t3 ) 784 785 !-------------------------------------------------------------------- 786 787 stiff(:, :) = 0.0D0 788 789 loopMatrix: do LX = 1, NumOfQuadPoints(etype) 790 791 !---------------------------------------------------------- 792 793 mu = gausses(LX)%pMaterial%variables(M_VISCOCITY) 794 795 rho = gausses(LX)%pMaterial%variables(M_DENSITY) 796 rho_inv = 1.0D0/rho 797 798 !---------------------------------------------------------- 799 800 call getQuadPoint( etype, LX, naturalCoord(:) ) 801 call getShapeFunc(etype, naturalCoord, spfunc) 802 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 803 804 !---------------------------------------------------------- 805 806 wg = getWeight(etype, LX)*det 807 808 !---------------------------------------------------------- 809 810 vx = 0.0D0 811 vy = 0.0D0 812 vz = 0.0D0 813 814 do na = 1, nn 815 816 vx = vx+spfunc(na)*v(1, na) 817 vy = vy+spfunc(na)*v(2, na) 818 vz = vz+spfunc(na)*v(3, na) 819 820 end do 821 822 !---------------------------------------------------------- 823 824 forall(na = 1:nn, nb = 1:nn) 825 826 MM(na, nb) = spfunc(na)*spfunc(nb) 827 AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) & 828 +vy*( spfunc(na)*gderiv(nb, 2) ) & 829 +vz*( spfunc(na)*gderiv(nb, 3) ) 830 DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1) 831 DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2) 832 DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3) 833 DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1) 834 DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2) 835 DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3) 836 DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1) 837 DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2) 838 DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3) 839 trD(na, nb) = DD(na, nb, 1, 1) & 840 +DD(na, nb, 2, 2) & 841 +DD(na, nb, 3, 3) 842 BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) & 843 +( vx*vy )*DD(na, nb, 1, 2) & 844 +( vx*vz )*DD(na, nb, 1, 3) & 845 +( vy*vx )*DD(na, nb, 2, 1) & 846 +( vy*vy )*DD(na, nb, 2, 2) & 847 +( vy*vz )*DD(na, nb, 2, 3) & 848 +( vz*vx )*DD(na, nb, 3, 1) & 849 +( vz*vy )*DD(na, nb, 3, 2) & 850 +( vz*vz )*DD(na, nb, 3, 3) 851 CC(na, nb, 1) = gderiv(na, 1)*spfunc(nb) 852 CC(na, nb, 2) = gderiv(na, 2)*spfunc(nb) 853 CC(na, nb, 3) = gderiv(na, 3)*spfunc(nb) 854 855 MS(nb, na) = AA(na, nb) 856 AS(na, nb) = BB(na, nb) 857 CS(na, nb, 1) = vx*DD(na, nb, 1, 1) & 858 +vy*DD(na, nb, 2, 1) & 859 +vz*DD(na, nb, 3, 1) 860 CS(na, nb, 2) = vx*DD(na, nb, 1, 2) & 861 +vy*DD(na, nb, 2, 2) & 862 +vz*DD(na, nb, 3, 2) 863 CS(na, nb, 3) = vx*DD(na, nb, 1, 3) & 864 +vy*DD(na, nb, 2, 3) & 865 +vz*DD(na, nb, 3, 3) 866 MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1) 867 MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2) 868 MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3) 869 AP(nb, na, 1) = CS(na, nb, 1) 870 AP(nb, na, 2) = CS(na, nb, 2) 871 AP(nb, na, 3) = CS(na, nb, 3) 872 CP(na, nb) = trD(na, nb) 873 874 end forall 875 876 !---------------------------------------------------------- 877 878 do nb = 1, nn 879 880 do na = 1, nn 881 882 i = 1 883 j = 1 884 isize = 4*(na-1)+i 885 jsize = 4*(nb-1)+j 886 887 stiff(isize, jsize) & 888 = stiff(isize, jsize) & 889 +wg & 890 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 891 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 892 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 893 894 i = 1 895 j = 2 896 isize = 4*(na-1)+i 897 jsize = 4*(nb-1)+j 898 899 stiff(isize, jsize) & 900 = stiff(isize, jsize) & 901 +wg & 902 *( gamma*mu*DD(na, nb, 2, 1) ) 903 904 i = 1 905 j = 3 906 isize = 4*(na-1)+i 907 jsize = 4*(nb-1)+j 908 909 stiff(isize, jsize) & 910 = stiff(isize, jsize) & 911 +wg & 912 *( gamma*mu*DD(na, nb, 3, 1) ) 913 914 i = 1 915 j = 4 916 isize = 4*(na-1)+i 917 jsize = 4*(nb-1)+j 918 919 stiff(isize, jsize) & 920 = stiff(isize, jsize) & 921 +wg & 922 *( -CC(na, nb, 1)+tau*CS(na, nb, 1) ) 923 924 i = 2 925 j = 1 926 isize = 4*(na-1)+i 927 jsize = 4*(nb-1)+j 928 929 stiff(isize, jsize) & 930 = stiff(isize, jsize) & 931 +wg & 932 *( gamma*mu*DD(na, nb, 1, 2) ) 933 934 i = 2 935 j = 2 936 isize = 4*(na-1)+i 937 jsize = 4*(nb-1)+j 938 939 stiff(isize, jsize) & 940 = stiff(isize, jsize) & 941 +wg & 942 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 943 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 944 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 945 946 i = 2 947 j = 3 948 isize = 4*(na-1)+i 949 jsize = 4*(nb-1)+j 950 951 stiff(isize, jsize) & 952 = stiff(isize, jsize) & 953 +wg & 954 *( gamma*mu*DD(na, nb, 3, 2) ) 955 956 i = 2 957 j = 4 958 isize = 4*(na-1)+i 959 jsize = 4*(nb-1)+j 960 961 stiff(isize, jsize) & 962 = stiff(isize, jsize) & 963 +wg & 964 *( -CC(na, nb, 2)+tau*CS(na, nb, 2) ) 965 966 i = 3 967 j = 1 968 isize = 4*(na-1)+i 969 jsize = 4*(nb-1)+j 970 971 stiff(isize, jsize) & 972 = stiff(isize, jsize) & 973 +wg & 974 *( gamma*mu*DD(na, nb, 1, 3) ) 975 976 i = 3 977 j = 2 978 isize = 4*(na-1)+i 979 jsize = 4*(nb-1)+j 980 981 stiff(isize, jsize) & 982 = stiff(isize, jsize) & 983 +wg & 984 *( gamma*mu*DD(na, nb, 2, 3) ) 985 986 i = 3 987 j = 3 988 isize = 4*(na-1)+i 989 jsize = 4*(nb-1)+j 990 991 stiff(isize, jsize) & 992 = stiff(isize, jsize) & 993 +wg & 994 *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) & 995 +gamma*rho*( AA(na, nb)+tau*AS(na, nb) ) & 996 +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) ) 997 998 i = 3 999 j = 4 1000 isize = 4*(na-1)+i 1001 jsize = 4*(nb-1)+j 1002 1003 stiff(isize, jsize) & 1004 = stiff(isize, jsize) & 1005 +wg & 1006 *( -CC(na, nb, 3)+tau*CS(na, nb, 3) ) 1007 1008 i = 4 1009 j = 1 1010 isize = 4*(na-1)+i 1011 jsize = 4*(nb-1)+j 1012 1013 stiff(isize, jsize) & 1014 = stiff(isize, jsize) & 1015 +wg & 1016 *( CC(nb, na, j) & 1017 +tincr_inv*tau*MP(na, nb, j) & 1018 +gamma*tau*AP(na, nb, j) ) 1019 1020 i = 4 1021 j = 2 1022 isize = 4*(na-1)+i 1023 jsize = 4*(nb-1)+j 1024 1025 stiff(isize, jsize) & 1026 = stiff(isize, jsize) & 1027 +wg & 1028 *( CC(nb, na, j) & 1029 +tincr_inv*tau*MP(na, nb, j) & 1030 +gamma*tau*AP(na, nb, j) ) 1031 1032 i = 4 1033 j = 3 1034 isize = 4*(na-1)+i 1035 jsize = 4*(nb-1)+j 1036 1037 stiff(isize, jsize) & 1038 = stiff(isize, jsize) & 1039 +wg & 1040 *( CC(nb, na, j) & 1041 +tincr_inv*tau*MP(na, nb, j) & 1042 +gamma*tau*AP(na, nb, j) ) 1043 1044 i = 4 1045 j = 4 1046 isize = 4*(na-1)+i 1047 jsize = 4*(nb-1)+j 1048 1049 stiff(isize, jsize) & 1050 = stiff(isize, jsize) & 1051 +wg & 1052 *( rho_inv*tau*trD(na, nb) ) 1053 1054 end do 1055 1056 end do 1057 1058 !---------------------------------------------------------- 1059 1060 end do loopMatrix 1061 1062 !-------------------------------------------------------------------- 1063 1064 b(:) = 0.0D0 1065 1066 loopVector: do LX = 1, NumOfQuadPoints(etype) 1067 1068 !---------------------------------------------------------- 1069 1070 mu = gausses(LX)%pMaterial%variables(M_VISCOCITY) 1071 1072 rho = gausses(LX)%pMaterial%variables(M_DENSITY) 1073 rho_inv = 1.0D0/rho 1074 1075 !---------------------------------------------------------- 1076 1077 call getQuadPoint( etype, LX, naturalCoord(:) ) 1078 call getShapeFunc(etype, naturalCoord, spfunc) 1079 call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv) 1080 1081 !---------------------------------------------------------- 1082 1083 wg = getWeight(etype, LX)*det 1084 1085 !---------------------------------------------------------- 1086 1087 vx = 0.0D0 1088 vy = 0.0D0 1089 vz = 0.0D0 1090 1091 do na = 1, nn 1092 1093 vx = vx+spfunc(na)*v(1, na) 1094 vy = vy+spfunc(na)*v(2, na) 1095 vz = vz+spfunc(na)*v(3, na) 1096 1097 end do 1098 1099 !---------------------------------------------------------- 1100 1101 forall(na = 1:nn, nb = 1:nn) 1102 1103 MM(na, nb) = spfunc(na)*spfunc(nb) 1104 AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) & 1105 +vy*( spfunc(na)*gderiv(nb, 2) ) & 1106 +vz*( spfunc(na)*gderiv(nb, 3) ) 1107 DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1) 1108 DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2) 1109 DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3) 1110 DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1) 1111 DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2) 1112 DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3) 1113 DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1) 1114 DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2) 1115 DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3) 1116 BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) & 1117 +( vx*vy )*DD(na, nb, 1, 2) & 1118 +( vx*vz )*DD(na, nb, 1, 3) & 1119 +( vy*vx )*DD(na, nb, 2, 1) & 1120 +( vy*vy )*DD(na, nb, 2, 2) & 1121 +( vy*vz )*DD(na, nb, 2, 3) & 1122 +( vz*vx )*DD(na, nb, 3, 1) & 1123 +( vz*vy )*DD(na, nb, 3, 2) & 1124 +( vz*vz )*DD(na, nb, 3, 3) 1125 1126 MS(nb, na) = AA(na, nb) 1127 AS(na, nb) = BB(na, nb) 1128 CS(na, nb, 1) = vx*DD(na, nb, 1, 1) & 1129 +vy*DD(na, nb, 2, 1) & 1130 +vz*DD(na, nb, 3, 1) 1131 CS(na, nb, 2) = vx*DD(na, nb, 1, 2) & 1132 +vy*DD(na, nb, 2, 2) & 1133 +vz*DD(na, nb, 3, 2) 1134 CS(na, nb, 3) = vx*DD(na, nb, 1, 3) & 1135 +vy*DD(na, nb, 2, 3) & 1136 +vz*DD(na, nb, 3, 3) 1137 MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1) 1138 MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2) 1139 MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3) 1140 AP(nb, na, 1) = CS(na, nb, 1) 1141 AP(nb, na, 2) = CS(na, nb, 2) 1142 AP(nb, na, 3) = CS(na, nb, 3) 1143 1144 end forall 1145 1146 !---------------------------------------------------------- 1147 1148 do na = 1, nn 1149 1150 !---------------------------------------------------- 1151 1152 do i = 1, 3 1153 1154 m_v(i) = 0.0D0 1155 a_v(i) = 0.0D0 1156 do j = 1, 3 1157 do k = 1, 3 1158 d_v(j, k, i) = 0.0D0 1159 end do 1160 end do 1161 ms_v(i) = 0.0D0 1162 as_v(i) = 0.0D0 1163 mp_dot_v = 0.0D0 1164 ap_dot_v = 0.0D0 1165 1166 do nb = 1, nn 1167 1168 ! Unsteady term 1169 m_v(i) = m_v(i)+MM(na, nb)*v(i, nb) 1170 ! Advection term 1171 a_v(i) = a_v(i)+AA(na, nb)*v(i, nb) 1172 ! Diffusion term 1173 do j = 1, 3 1174 do k = 1, 3 1175 d_v(j, k, i) = d_v(j, k, i)+DD(na, nb, j, k)*v(i, nb) 1176 end do 1177 end do 1178 ! Unsteady term (SUPG) 1179 ms_v(i) = ms_v(i)+MS(na, nb)*v(i, nb) 1180 ! Advection term (SUPG) 1181 as_v(i) = as_v(i)+AS(na, nb)*v(i, nb) 1182 ! Unsteady term (PSPG) 1183 mp_dot_v = mp_dot_v+( MP(na, nb, 1)*v(1, nb) & 1184 +MP(na, nb, 2)*v(2, nb) & 1185 +MP(na, nb, 3)*v(3, nb) ) 1186 ! Advection term (PSPG) 1187 ap_dot_v = ap_dot_v+( AP(na, nb, 1)*v(1, nb) & 1188 +AP(na, nb, 2)*v(2, nb) & 1189 +AP(na, nb, 3)*v(3, nb) ) 1190 end do 1191 1192 end do 1193 1194 !---------------------------------------------------- 1195 1196 do i = 1, 3 1197 1198 isize = 4*(na-1)+i 1199 1200 b(isize) & 1201 = b(isize) & 1202 +wg & 1203 *( tincr_inv*rho*( m_v(i)+tau*ms_v(i) ) & 1204 -( 1.0D0-gamma )*rho*( a_v(i)+tau*as_v(i) ) & 1205 -( 1.0D0-gamma ) & 1206 *mu*( ( d_v(1, 1, i)+d_v(2, 2, i)+d_v(3, 3, i) ) & 1207 +( d_v(1, i, 1)+d_v(2, i, 2)+d_v(3, i, 3) ) ) ) 1208 1209 end do 1210 1211 i = 4 1212 isize = 4*(na-1)+i 1213 1214 b(isize) & 1215 = b(isize) & 1216 +wg & 1217 *( tincr_inv*tau*( mp_dot_v ) & 1218 -( 1.0D0-gamma )*tau*( ap_dot_v ) ) 1219 1220 !---------------------------------------------------- 1221 1222 end do 1223 1224 !---------------------------------------------------------- 1225 1226 end do loopVector 1227 1228 !---------------------------------------------------------------- 1229 1230 do isize = 1, 4*nn 1231 1232 stiff_velo = 0.0D0 1233 1234 do jsize = 1, 4*nn 1235 1236 stiff_velo = stiff_velo+stiff(isize, jsize)*velo_new(jsize) 1237 1238 end do 1239 1240 r(isize) = b(isize)-stiff_velo 1241 1242 end do 1243 1244 !-------------------------------------------------------------------- 1245 end subroutine LOAD_C3_vp 1246 !-------------------------------------------------------------------- 1247 1248 1249end module m_static_LIB_3d_vp 1250