1 subroutine bc3LHS (iBC, BC, iens, xKebe ) 2c 3c---------------------------------------------------------------------- 4c 5c This routine satisfies the BC of LHS mass matrix for all 6c elements in this block. 7c 8c input: 9c iBC (nshg) : boundary condition code 10c BC (nshg,ndofBC) : Dirichlet BC constraint parameters 11c ien (npro,nshape) : ien array for this element 12c xKebe (npro,9,nshl,nshl) : element consistent mass matrix before BC 13c 14c output: 15c xKebe (npro,9,nshl,nshl) : LHS mass matrix after BC is satisfied 16c 17c 18c Farzin Shakib, Winter 1987. 19c Zdenek Johan, Spring 1990. (Modified for general divariant gas) 20c Ken Jansen, Summer 2000. Incompressible (only needed on xKebe) 21c---------------------------------------------------------------------- 22c 23 include "common.h" 24c 25 dimension iBC(nshg), ien(npro,nshl), 26 & BC(nshg,ndofBC), xKebe(npro,9,nshl,nshl) 27 integer iens(npro,nshl) 28c 29c prefer to show explicit absolute value needed for cubic modes and 30c higher rather than inline abs on pointer as in past versions 31c iens is the signed ien array ien is unsigned 32c 33 ien=abs(iens) 34c 35c.... loop over elements 36c 37c return 38 do iel = 1, npro 39c 40c.... loop over number of shape functions for this element 41c 42 do inod = 1, nshl 43c 44c.... set up parameters 45c 46 in = abs(ien(iel,inod)) 47 if (ibits(iBC(in),3,3) .eq. 0) goto 5000 ! NO velocity BC's 48 if (ibits(iBC(in),3,3) .eq. 7) goto 5000 ! 3 components ok 49 50c.... 1 or 2 component velocities 51c 52c 53c.... x1-velocity 54c 55 if ( ibits(iBC(in),3,3) .eq. 1) then 56c 57! we want to project out the x1 component of the velocity from the tangent 58! matix which is, mathematically, M^e = S^T M^e S. We will do the M^e S 59! product first. It has the effect of 60! subtracting the column of the block-9 matrix from each column of the block-9 61! matrix that is going to survive (weighted by the coefficient in the 62! BC array associated with that row) FOR EACH column of the 63! nshl by nshl matrix FOR EACH element. THEN the transpose of the 64! operation is carried out (replace the word "column" by row 65! EVERYWHERE). The following code has been set up so that we only have to 66! give the starting position in each case since we know the block-9 matrix is 67! ordered like this 68! 1 2 3 69! 4 5 6 70! 7 8 9 71 72c 73c adjusting the second column for the eventual removal of the first 74c column of the block-9 submatrix 75c 76 irem1=1 77 irem2=irem1+3 78 irem3=irem2+3 79 80 iadj1=2 81 iadj2=iadj1+3 82 iadj3=iadj2+3 83 do i = 1, nshl 84 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 85 & - BC(in,4) * xKebe(iel,irem1,i,inod) 86 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 87 & - BC(in,4) * xKebe(iel,irem2,i,inod) 88 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 89 & - BC(in,4) * xKebe(iel,irem3,i,inod) 90 91 enddo 92! block status ' denotes colunn 1 projected off. 93! 1 2' 3 94! 4 5' 6 95! 7 8' 9 96c 97c adjusting the third column for the eventual removal of the first 98c column of the block-9 submatrix 99c 100 iadj1=3 101 iadj2=iadj1+3 102 iadj3=iadj2+3 103 do i = 1, nshl 104 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 105 & - BC(in,5) * xKebe(iel,irem1,i,inod) 106 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 107 & - BC(in,5) * xKebe(iel,irem2,i,inod) 108 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 109 & - BC(in,5) * xKebe(iel,irem3,i,inod) 110! block status 111! 1 2' 3' 112! 4 5' 6' 113! 7 8' 9' 114 enddo 115 do i=1,nshl 116c 117c done with the first columns_block-9 for columns AND rows of nshl 118c 119 xKebe(iel,irem1,i,inod) = zero 120 xKebe(iel,irem2,i,inod) = zero 121 xKebe(iel,irem3,i,inod) = zero 122 123 124! block status 125! 0 2' 3' 126! 0 5' 6' 127! 0 8' 9' 128 129 enddo 130c 131c now adjust the second row_block-9 for EACH row nshl for EACH element 132c 133 134 iadj1=4 135 iadj2=iadj1+1 136 iadj3=iadj2+1 137 irem1=1 138 irem2=irem1+1 139 irem3=irem2+1 140 do i = 1, nshl 141 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 142 & - BC(in,4) * xKebe(iel,irem1,inod,i) 143 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 144 & - BC(in,4) * xKebe(iel,irem2,inod,i) 145 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 146 & - BC(in,4) * xKebe(iel,irem3,inod,i) 147 148 enddo 149! block status 150! 0 2' 3' 151! 0 5'' 6'' 152! 0 8' 9' 153 154 155 iadj1=7 156 iadj2=iadj1+1 157 iadj3=iadj2+1 158 do i = 1, nshl 159 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 160 & - BC(in,5) * xKebe(iel,irem1,inod,i) 161 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 162 & - BC(in,5) * xKebe(iel,irem2,inod,i) 163 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 164 & - BC(in,5) * xKebe(iel,irem3,inod,i) 165 166! block status 167! 0 2' 3' 168! 0 5'' 6'' 169! 0 8'' 9'' 170 enddo 171 do i=1,nshl 172 173c 174c eliminate the first row of block-9 for all rows 175c 176 xKebe(iel,irem1,inod,i) = zero 177 xKebe(iel,irem2,inod,i) = zero 178 xKebe(iel,irem3,inod,i) = zero 179 180 enddo 181 182! block status 183! 0 0 0 184! 0 5'' 6'' 185! 0 8'' 9'' 186 187! Be aware that this simple status of the block does not reflect that when 188! we eliminated columns we did it for columns in nshl as well for the given 189! inod. Conversely when we eliminated rows in the block we did so for ALL 190! rows in nshl as can be seen by the transpose of i and inod. 191 192 xKebe(iel,1,inod,inod)=one 193! block status 194! 1 0 0 195! 0 5'' 6'' 196! 0 8'' 9'' 197 endif 198c 199c.... x2-velocity 200c 201 if ( ibits(iBC(in),3,3) .eq. 2) then 202c 203! See comment above. Now we are eliminating the 2nd column then row of 204! the block-9 matrix 205! 1 2 3 206! 4 5 6 207! 7 8 9 208c 209c adjusting the first column for the eventual removal of the second 210c column of the block-9 submatrix 211c 212 irem1=2 213 irem2=irem1+3 214 irem3=irem2+3 215 216 iadj1=1 217 iadj2=iadj1+3 218 iadj3=iadj2+3 219 do i = 1, nshl 220 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 221 & - BC(in,4) * xKebe(iel,irem1,i,inod) 222 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 223 & - BC(in,4) * xKebe(iel,irem2,i,inod) 224 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 225 & - BC(in,4) * xKebe(iel,irem3,i,inod) 226 227 enddo 228c 229c adjusting the third column for the eventual removal of the second 230c column of the block-9 submatrix 231c 232 iadj1=3 233 iadj2=iadj1+3 234 iadj3=iadj2+3 235 do i = 1, nshl 236 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 237 & - BC(in,5) * xKebe(iel,irem1,i,inod) 238 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 239 & - BC(in,5) * xKebe(iel,irem2,i,inod) 240 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 241 & - BC(in,5) * xKebe(iel,irem3,i,inod) 242 243 enddo 244 do i=1,nshl 245c 246c done with the second columns_block-9 for columns 247c 248 249 xKebe(iel,irem1,i,inod) = zero 250 xKebe(iel,irem2,i,inod) = zero 251 xKebe(iel,irem3,i,inod) = zero 252 253 enddo 254c 255c now adjust the 1st row_block-9 for EACH row nshl for EACH element 256c 257 258 iadj1=1 259 iadj2=iadj1+1 260 iadj3=iadj2+1 261 irem1=4 262 irem2=irem1+1 263 irem3=irem2+1 264 do i = 1, nshl 265 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 266 & - BC(in,4) * xKebe(iel,irem1,inod,i) 267 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 268 & - BC(in,4) * xKebe(iel,irem2,inod,i) 269 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 270 & - BC(in,4) * xKebe(iel,irem3,inod,i) 271 272 enddo 273 iadj1=7 274 iadj2=iadj1+1 275 iadj3=iadj2+1 276 do i = 1, nshl 277 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 278 & - BC(in,5) * xKebe(iel,irem1,inod,i) 279 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 280 & - BC(in,5) * xKebe(iel,irem2,inod,i) 281 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 282 & - BC(in,5) * xKebe(iel,irem3,inod,i) 283 enddo 284 do i=1,nshl 285 286c 287c eliminate the second row of block-9 for all rows 288c 289 xKebe(iel,irem1,inod,i) = zero 290 xKebe(iel,irem2,inod,i) = zero 291 xKebe(iel,irem3,inod,i) = zero 292 enddo 293 xKebe(iel,5,inod,inod)=one 294 endif 295c 296c.... x3-velocity 297c 298 if ( ibits(iBC(in),3,3) .eq. 4) then 299c 300! See comment above. Now we are eliminating the 3rd column then row of 301! the block-9 matrix 302! 1 2 3 303! 4 5 6 304! 7 8 9 305c 306c adjusting the 1st column for the eventual removal of the 3rd 307c column of the block-9 submatrix 308c 309 irem1=3 310 irem2=irem1+3 311 irem3=irem2+3 312 313 iadj1=1 314 iadj2=iadj1+3 315 iadj3=iadj2+3 316 do i = 1, nshl 317 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 318 & - BC(in,4) * xKebe(iel,irem1,i,inod) 319 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 320 & - BC(in,4) * xKebe(iel,irem2,i,inod) 321 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 322 & - BC(in,4) * xKebe(iel,irem3,i,inod) 323 324 enddo 325c 326c adjusting the second column for the eventual removal of the 3rd 327c column of the block-9 submatrix 328c 329 iadj1=2 330 iadj2=iadj1+3 331 iadj3=iadj2+3 332 do i = 1, nshl 333 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 334 & - BC(in,5) * xKebe(iel,irem1,i,inod) 335 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 336 & - BC(in,5) * xKebe(iel,irem2,i,inod) 337 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 338 & - BC(in,5) * xKebe(iel,irem3,i,inod) 339 enddo 340 do i=1,nshl 341 342c 343c done with the 3rd columns_block-9 for columns 344c 345 346 xKebe(iel,irem1,i,inod) = zero 347 xKebe(iel,irem2,i,inod) = zero 348 xKebe(iel,irem3,i,inod) = zero 349 350 enddo 351c 352c now adjust the 1st row_block-9 for EACH row nshl for EACH element 353c 354 355 iadj1=1 356 iadj2=iadj1+1 357 iadj3=iadj2+1 358 irem1=7 359 irem2=irem1+1 360 irem3=irem2+1 361 do i = 1, nshl 362 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 363 & - BC(in,4) * xKebe(iel,irem1,inod,i) 364 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 365 & - BC(in,4) * xKebe(iel,irem2,inod,i) 366 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 367 & - BC(in,4) * xKebe(iel,irem3,inod,i) 368 369 enddo 370 iadj1=4 371 iadj2=iadj1+1 372 iadj3=iadj2+1 373 do i = 1, nshl 374 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 375 & - BC(in,5) * xKebe(iel,irem1,inod,i) 376 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 377 & - BC(in,5) * xKebe(iel,irem2,inod,i) 378 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 379 & - BC(in,5) * xKebe(iel,irem3,inod,i) 380 381 enddo 382 do i=1,nshl 383 xKebe(iel,irem1,inod,i) = zero 384 xKebe(iel,irem2,inod,i) = zero 385 xKebe(iel,irem3,inod,i) = zero 386 387 enddo 388 xKebe(iel,9,inod,inod)=one 389 endif 390c 391c.... x1-velocity and x2-velocity 392c 393 if ( ibits(iBC(in),3,3) .eq. 3 ) then 394c 395! See comment above. Now we are eliminating the 2nd and 1st column then 396! same rows of 397! the block-9 matrix 398! 1 2 3 399! 4 5 6 400! 7 8 9 401c 402c adjusting the 3rd column for the eventual removal of the first and second 403c column of the block-9 submatrix 404c 405 irem1=1 406 irem2=irem1+3 407 irem3=irem2+3 408 409 ire21=2 410 ire22=ire21+3 411 ire23=ire22+3 412 413 iadj1=3 414 iadj2=iadj1+3 415 iadj3=iadj2+3 416 do i = 1, nshl 417 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 418 & - BC(in,4) * xKebe(iel,irem1,i,inod) 419 & - BC(in,6) * xKebe(iel,ire21,i,inod) 420 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 421 & - BC(in,4) * xKebe(iel,irem2,i,inod) 422 & - BC(in,6) * xKebe(iel,ire22,i,inod) 423 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 424 & - BC(in,4) * xKebe(iel,irem3,i,inod) 425 & - BC(in,6) * xKebe(iel,ire23,i,inod) 426 427! Status of the block-9 matrix 428! 1 2 3' 429! 4 5 6' 430! 7 8 9' 431 enddo 432 do i=1,nshl 433c 434c done with the first and second columns_block-9 for columns AND rows of nshl 435c 436 437 xKebe(iel,irem1,i,inod) = zero 438 xKebe(iel,irem2,i,inod) = zero 439 xKebe(iel,irem3,i,inod) = zero 440 441 xKebe(iel,ire21,i,inod) = zero 442 xKebe(iel,ire22,i,inod) = zero 443 xKebe(iel,ire23,i,inod) = zero 444 445! Status of the block-9 matrix 446! 0 0 3' 447! 0 0 6' 448! 0 0 9' 449 450 enddo 451c 452c now adjust the 3rd row_block-9 for EACH row nshl for EACH element 453c 454 455 iadj1=7 456 iadj2=iadj1+1 457 iadj3=iadj2+1 458 irem1=1 459 irem2=irem1+1 460 irem3=irem2+1 461 ire21=4 462 ire22=ire21+1 463 ire23=ire22+1 464 do i = 1, nshl 465 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 466 & - BC(in,4) * xKebe(iel,irem1,inod,i) 467 & - BC(in,6) * xKebe(iel,ire21,inod,i) 468 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 469 & - BC(in,4) * xKebe(iel,irem2,inod,i) 470 & - BC(in,6) * xKebe(iel,ire22,inod,i) 471 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 472 & - BC(in,4) * xKebe(iel,irem3,inod,i) 473 & - BC(in,6) * xKebe(iel,ire23,inod,i) 474 475 476! Status of the block-9 matrix 477! 0 0 3' 478! 0 0 6' 479! 0 0 9'' 480 enddo 481 do i=1,nshl 482 xKebe(iel,irem1,inod,i) = zero 483 xKebe(iel,irem2,inod,i) = zero 484 xKebe(iel,irem3,inod,i) = zero 485 486 xKebe(iel,ire21,inod,i) = zero 487 xKebe(iel,ire22,inod,i) = zero 488 xKebe(iel,ire23,inod,i) = zero 489 490! Status of the block-9 matrix 491! 0 0 0 492! 0 0 0 493! 0 0 9'' 494 495 enddo 496 xKebe(iel,1,inod,inod)=one 497 xKebe(iel,5,inod,inod)=one 498 endif 499c 500c.... x1-velocity and x3-velocity 501c 502 if ( ibits(iBC(in),3,3) .eq. 5 ) then 503c 504! See comment above. Now we are eliminating the 1 and 3 column then 505! same rows of 506! the block-9 matrix 507! 1 2 3 508! 4 5 6 509! 7 8 9 510c 511c adjusting the 3rd column for the eventual removal of the first and second 512c column of the block-9 submatrix 513c 514 irem1=1 515 irem2=irem1+3 516 irem3=irem2+3 517 518 ire21=3 519 ire22=ire21+3 520 ire23=ire22+3 521 522 iadj1=2 523 iadj2=iadj1+3 524 iadj3=iadj2+3 525 do i = 1, nshl 526 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 527 & - BC(in,4) * xKebe(iel,irem1,i,inod) 528 & - BC(in,6) * xKebe(iel,ire21,i,inod) 529 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 530 & - BC(in,4) * xKebe(iel,irem2,i,inod) 531 & - BC(in,6) * xKebe(iel,ire22,i,inod) 532 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 533 & - BC(in,4) * xKebe(iel,irem3,i,inod) 534 & - BC(in,6) * xKebe(iel,ire23,i,inod) 535 536 enddo 537 do i=1,nshl 538c 539c done with the first and third columns_block-9 for columns AND rows of nshl 540c 541 xKebe(iel,irem1,i,inod) = zero 542 xKebe(iel,irem2,i,inod) = zero 543 xKebe(iel,irem3,i,inod) = zero 544 545 xKebe(iel,ire21,i,inod) = zero 546 xKebe(iel,ire22,i,inod) = zero 547 xKebe(iel,ire23,i,inod) = zero 548 enddo 549c 550c now adjust the 2nd row_block-9 for EACH row nshl for EACH element 551c 552 553 iadj1=4 554 iadj2=iadj1+1 555 iadj3=iadj2+1 556 irem1=1 557 irem2=irem1+1 558 irem3=irem2+1 559 ire21=7 560 ire22=ire21+1 561 ire23=ire22+1 562 do i = 1, nshl 563 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 564 & - BC(in,4) * xKebe(iel,irem1,inod,i) 565 & - BC(in,6) * xKebe(iel,ire21,inod,i) 566 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 567 & - BC(in,4) * xKebe(iel,irem2,inod,i) 568 & - BC(in,6) * xKebe(iel,ire22,inod,i) 569 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 570 & - BC(in,4) * xKebe(iel,irem3,inod,i) 571 & - BC(in,6) * xKebe(iel,ire23,inod,i) 572 573 enddo 574 do i=1,nshl 575 xKebe(iel,irem1,inod,i) = zero 576 xKebe(iel,irem2,inod,i) = zero 577 xKebe(iel,irem3,inod,i) = zero 578 579 xKebe(iel,ire21,inod,i) = zero 580 xKebe(iel,ire22,inod,i) = zero 581 xKebe(iel,ire23,inod,i) = zero 582 583 enddo 584 xKebe(iel,1,inod,inod)=one 585 xKebe(iel,9,inod,inod)=one 586 endif 587c 588c.... x2-velocity and x3-velocity 589c 590 if ( ibits(iBC(in),3,3) .eq. 6 ) then 591c 592! See comment above. Now we are eliminating the 2nd and 3rd column then 593! same rows of 594! the block-9 matrix 595! 1 2 3 596! 4 5 6 597! 7 8 9 598c 599c adjusting the 3rd column for the eventual removal of the first and second 600c column of the block-9 submatrix 601c 602 irem1=2 603 irem2=irem1+3 604 irem3=irem2+3 605 606 ire21=3 607 ire22=ire21+3 608 ire23=ire22+3 609 610 iadj1=1 611 iadj2=iadj1+3 612 iadj3=iadj2+3 613 do i = 1, nshl 614 xKebe(iel,iadj1,i,inod) = xKebe(iel,iadj1,i,inod) 615 & - BC(in,4) * xKebe(iel,irem1,i,inod) 616 & - BC(in,6) * xKebe(iel,ire21,i,inod) 617 xKebe(iel,iadj2,i,inod) = xKebe(iel,iadj2,i,inod) 618 & - BC(in,4) * xKebe(iel,irem2,i,inod) 619 & - BC(in,6) * xKebe(iel,ire22,i,inod) 620 xKebe(iel,iadj3,i,inod) = xKebe(iel,iadj3,i,inod) 621 & - BC(in,4) * xKebe(iel,irem3,i,inod) 622 & - BC(in,6) * xKebe(iel,ire23,i,inod) 623 enddo 624 do i=1,nshl 625 626c 627c done with the first and second columns_block-9 for columns AND rows of nshl 628c 629 xKebe(iel,irem1,i,inod) = zero 630 xKebe(iel,irem2,i,inod) = zero 631 xKebe(iel,irem3,i,inod) = zero 632 633 xKebe(iel,ire21,i,inod) = zero 634 xKebe(iel,ire22,i,inod) = zero 635 xKebe(iel,ire23,i,inod) = zero 636 637 enddo 638c 639c now adjust the 3rd row_block-9 for EACH row nshl for EACH element 640c 641 642 iadj1=1 643 iadj2=iadj1+1 644 iadj3=iadj2+1 645 irem1=4 646 irem2=irem1+1 647 irem3=irem2+1 648 ire21=7 649 ire22=ire21+1 650 ire23=ire22+1 651 do i = 1, nshl 652 xKebe(iel,iadj1,inod,i) = xKebe(iel,iadj1,inod,i) 653 & - BC(in,4) * xKebe(iel,irem1,inod,i) 654 xKebe(iel,iadj2,inod,i) = xKebe(iel,iadj2,inod,i) 655 & - BC(in,4) * xKebe(iel,irem2,inod,i) 656 xKebe(iel,iadj3,inod,i) = xKebe(iel,iadj3,inod,i) 657 & - BC(in,4) * xKebe(iel,irem3,inod,i) 658 & - BC(in,6) * xKebe(iel,ire23,inod,i) 659 660 enddo 661 do i=1,nshl 662 xKebe(iel,irem1,inod,i) = zero 663 xKebe(iel,irem2,inod,i) = zero 664 xKebe(iel,irem3,inod,i) = zero 665 666c 667 xKebe(iel,ire21,inod,i) = zero 668 xKebe(iel,ire22,inod,i) = zero 669 xKebe(iel,ire23,inod,i) = zero 670 671 enddo 672 xKebe(iel,5,inod,inod)=one 673 xKebe(iel,9,inod,inod)=one 674 endif 675 676 5000 continue 677 678c 679c.... end loop over shape functions (nodes) 680c 681 enddo 682c 683c.... end loop over elements 684c 685 enddo 686c 687c These elements should assemble to a matrix with the rows and columns 688c associated with the Dirichlet nodes zeroed out. Note that BC3 Diag 689c 690c 691c.... return 692c 693 return 694 end 695