1 subroutine bc3LHS (iBC, BC, iens, EGmass ) 2c 3c---------------------------------------------------------------------- 4c 5c This routine satisfies the BC of LHS mass matrix for a single 6c element. 7c 8c input: 9c iBC (nshg) : boundary condition code 10c BC (nshg,11) : Dirichlet BC constraint parameters 11c ien (npro,nshl) : ien array for this element 12c EGmass(npro,nedof,nedof) : element consistent mass matrix before BC 13c 14c output: 15c EGmass(npro,nedof,nedof): 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---------------------------------------------------------------------- 21c 22 include "common.h" 23c 24 dimension iBC(nshg), ien(npro,nshl), 25 & BC(nshg,11), EGmass(npro,nedof,nedof) 26 integer iens(npro,nshl) 27c 28c prefer to show explicit absolute value needed for cubic modes and 29c higher rather than inline abs on pointer as in past versions 30c iens is the signed ien array ien is unsigned 31c 32 ien=abs(iens) 33c 34c.... loop over elements 35c 36 do iel = 1, npro 37c 38c.... loop over number of shape functions for this element 39c 40 do inod = 1, nshl 41c 42c.... set up parameters 43c 44 in = ien(iel,inod) 45 if (iBC(in) .eq. 0) goto 5000 46 ioff = (inod - 1) * nflow 47 i1 = ioff + 1 48 i2 = ioff + 2 49 i3 = ioff + 3 50 i4 = ioff + 4 51 i5 = ioff + 5 52c 53c.... pressure 54c 55 if ( btest(iBC(in),2) ) then 56 57 do i = 1, nedof 58 EGmass(iel,i,i1) = zero 59 EGmass(iel,i1,i) = zero 60 enddo 61 EGmass(iel,i1,i1) = one 62 endif 63c 64c.... density ( not activated yet for LHS matrix ) 65c 66 67c 68c.... velocities 69c 70c 71c.... x1-velocity 72c 73 if ( ibits(iBC(in),3,3) .eq. 1) then 74 do i = 1, nedof 75 EGmass(iel,i3,i) = EGmass(iel,i3,i) 76 & - BC(in,4) * EGmass(iel,i2,i) 77 EGmass(iel,i4,i) = EGmass(iel,i4,i) 78 & - BC(in,5) * EGmass(iel,i2,i) 79 enddo 80 do i = 1, nedof 81 EGmass(iel,i,i3) = EGmass(iel,i,i3) 82 & - BC(in,4) * EGmass(iel,i,i2) 83 EGmass(iel,i,i4) = EGmass(iel,i,i4) 84 & - BC(in,5) * EGmass(iel,i,i2) 85 enddo 86 do i = 1, nedof 87 EGmass(iel,i,i2) = zero 88 EGmass(iel,i2,i) = zero 89 enddo 90 EGmass(iel,i2,i2) = one 91 endif 92c 93c.... x2-velocity 94c 95 if ( ibits(iBC(in),3,3) .eq. 2 ) then 96 do i = 1, nedof 97 EGmass(iel,i2,i) = EGmass(iel,i2,i) 98 & - BC(in,4) * EGmass(iel,i3,i) 99 EGmass(iel,i4,i) = EGmass(iel,i4,i) 100 & - BC(in,5) * EGmass(iel,i3,i) 101 enddo 102 do i = 1, nedof 103 EGmass(iel,i,i2) = EGmass(iel,i,i2) 104 & - BC(in,4) * EGmass(iel,i,i3) 105 EGmass(iel,i,i4) = EGmass(iel,i,i4) 106 & - BC(in,5) * EGmass(iel,i,i3) 107 enddo 108 109 do i = 1, nedof 110 EGmass(iel,i,i3) = zero 111 EGmass(iel,i3,i) = zero 112 enddo 113 EGmass(iel,i3,i3) = one 114 endif 115c 116c.... x1-velocity and x2-velocity 117c 118 if ( ibits(iBC(in),3,3) .eq. 3 ) then 119 do i = 1, nedof 120 EGmass(iel,i4,i) = EGmass(iel,i4,i) 121 & - BC(in,4) * EGmass(iel,i2,i) 122 & - BC(in,6) * EGmass(iel,i3,i) 123 enddo 124 do i = 1, nedof 125 EGmass(iel,i,i4) = EGmass(iel,i,i4) 126 & - BC(in,4) * EGmass(iel,i,i2) 127 & - BC(in,6) * EGmass(iel,i,i3) 128 enddo 129 130 do i = 1, nedof 131 EGmass(iel,i,i2) = zero 132 EGmass(iel,i2,i) = zero 133 EGmass(iel,i,i3) = zero 134 EGmass(iel,i3,i) = zero 135 enddo 136 EGmass(iel,i2,i2) = one 137 EGmass(iel,i3,i3) = one 138 endif 139c 140c.... x3-velocity 141c 142 if ( ibits(iBC(in),3,3) .eq. 4 ) then 143 do i = 1, nedof 144 EGmass(iel,i2,i) = EGmass(iel,i2,i) 145 & - BC(in,4) * EGmass(iel,i4,i) 146 EGmass(iel,i3,i) = EGmass(iel,i3,i) 147 & - BC(in,5) * EGmass(iel,i4,i) 148 enddo 149 do i = 1, nedof 150 EGmass(iel,i,i2) = EGmass(iel,i,i2) 151 & - BC(in,4) * EGmass(iel,i,i4) 152 EGmass(iel,i,i3) = EGmass(iel,i,i3) 153 & - BC(in,5) * EGmass(iel,i,i4) 154 enddo 155 156 do i = 1, nedof 157 EGmass(iel,i,i4) = zero 158 EGmass(iel,i4,i) = zero 159 enddo 160 EGmass(iel,i4,i4) = one 161 endif 162c 163c.... x1-velocity and x3-velocity 164c 165 if ( ibits(iBC(in),3,3) .eq. 5 ) then 166 do i = 1, nedof 167 EGmass(iel,i3,i) = EGmass(iel,i3,i) 168 & - BC(in,4) * EGmass(iel,i2,i) 169 & - BC(in,6) * EGmass(iel,i4,i) 170 enddo 171 do i = 1, nedof 172 EGmass(iel,i,i3) = EGmass(iel,i,i3) 173 & - BC(in,4) * EGmass(iel,i,i2) 174 & - BC(in,6) * EGmass(iel,i,i4) 175 enddo 176 177 do i = 1, nedof 178 EGmass(iel,i ,i2) = zero 179 EGmass(iel,i2,i ) = zero 180 EGmass(iel,i ,i4) = zero 181 EGmass(iel,i4,i ) = zero 182 enddo 183 EGmass(iel,i2,i2) = one 184 EGmass(iel,i4,i4) = one 185 endif 186c 187c.... x2-velocity and x3-velocity 188c 189 if ( ibits(iBC(in),3,3) .eq. 6 ) then 190 do i = 1, nedof 191 EGmass(iel,i2,i) = EGmass(iel,i2,i) 192 & - BC(in,4) * EGmass(iel,i3,i) 193 & - BC(in,6) * EGmass(iel,i4,i) 194 enddo 195 do i = 1, nedof 196 EGmass(iel,i,i2) = EGmass(iel,i,i2) 197 & - BC(in,4) * EGmass(iel,i,i3) 198 & - BC(in,6) * EGmass(iel,i,i4) 199 enddo 200 201 do i = 1, nedof 202 EGmass(iel,i ,i3) = zero 203 EGmass(iel,i3,i ) = zero 204 EGmass(iel,i ,i4) = zero 205 EGmass(iel,i4,i ) = zero 206 enddo 207 EGmass(iel,i3,i3) = one 208 EGmass(iel,i4,i4) = one 209 endif 210c 211c.... x1-velocity, x2-velocity, and x3-velocity 212c 213 if ( ibits(iBC(in),3,3) .eq. 7 ) then 214 do i = 1, nedof 215 EGmass(iel,i ,i2) = zero 216 EGmass(iel,i2,i ) = zero 217 EGmass(iel,i ,i3) = zero 218 EGmass(iel,i3,i ) = zero 219 EGmass(iel,i ,i4) = zero 220 EGmass(iel,i4,i ) = zero 221 enddo 222 EGmass(iel,i2,i2) = one 223 EGmass(iel,i3,i3) = one 224 EGmass(iel,i4,i4) = one 225 endif 226c 227c.... temperature 228c 229 if ( btest(iBC(in),1) ) then 230 do i = 1, nedof 231 EGmass(iel,i,i5) = zero 232 EGmass(iel,i5,i) = zero 233 enddo 234 EGmass(iel,i5,i5) = one 235 endif 236c 237c Elaine 238c.... scaled plane extraction boundary conditions 239c 240 if(intpres.eq.1) then 241 if ( btest(iBC(in),11) ) then 242 do i = 1, nedof 243 EGmass(iel,i ,i1) = zero 244 EGmass(iel,i1,i ) = zero 245 EGmass(iel,i ,i2) = zero 246 EGmass(iel,i2,i ) = zero 247 EGmass(iel,i ,i3) = zero 248 EGmass(iel,i3,i ) = zero 249 EGmass(iel,i ,i4) = zero 250 EGmass(iel,i4,i ) = zero 251 EGmass(iel,i,i5) = zero 252 EGmass(iel,i5,i) = zero 253 enddo 254 EGmass(iel,i1,i1) = one 255 EGmass(iel,i2,i2) = one 256 EGmass(iel,i3,i3) = one 257 EGmass(iel,i4,i4) = one 258 EGmass(iel,i5,i5) = one 259 endif 260 else 261 if ( btest(iBC(in),11) ) then 262 do i = 1, nedof 263 EGmass(iel,i ,i2) = zero 264 EGmass(iel,i2,i ) = zero 265 EGmass(iel,i ,i3) = zero 266 EGmass(iel,i3,i ) = zero 267 EGmass(iel,i ,i4) = zero 268 EGmass(iel,i4,i ) = zero 269 EGmass(iel,i,i5) = zero 270 EGmass(iel,i5,i) = zero 271 enddo 272 EGmass(iel,i2,i2) = one 273 EGmass(iel,i3,i3) = one 274 EGmass(iel,i4,i4) = one 275 EGmass(iel,i5,i5) = one 276 endif 277 endif 278 279 5000 continue 280 281c 282c.... end loop over shape functions (nodes) 283c 284 enddo 285c 286c.... end loop over elements 287c 288 enddo 289c 290c.... return 291c 292 return 293 end 294c 295c 296c 297 subroutine bc3LHSSclr (iBC, iens, EGmasst ) 298c 299c---------------------------------------------------------------------- 300c 301c This routine satisfies the BC of LHS mass matrix for a single 302c element. 303c 304c input: 305c iBC (nshg) : boundary condition code 306c BC (nshg,11) : Dirichlet BC constraint parameters 307c ien (npro,nshl) : ien array for this element 308c EGmasst(npro,nshape,nshape) : element consistent mass matrix before BC 309c 310c output: 311c EGmasst(npro,nshape,nshape): LHS mass matrix after BC is satisfied 312c 313c 314c Farzin Shakib, Winter 1987. 315c Zdenek Johan, Spring 1990. (Modified for general divariant gas) 316c---------------------------------------------------------------------- 317c 318 include "common.h" 319c 320 dimension iBC(nshg), ien(npro,nshl), 321 & EGmasst(npro,nshape,nshape) 322c 323 integer iens(npro,nshl) 324c 325c prefer to show explicit absolute value needed for cubic modes and 326c higher rather than inline abs on pointer as in past versions 327c iens is the signed ien array ien is unsigned 328c 329 ien=abs(iens) 330c 331 id = isclr+5 332c.... loop over elements 333c 334 do iel = 1, npro 335c 336c.... loop over number of shape functions for this element 337c 338 do inod = 1, nshl 339c 340c.... set up parameters 341c 342 in = ien(iel,inod) 343 if (iBC(in) .eq. 0) goto 5000 344 345c 346c.... scalar 347c 348 if ( btest(iBC(in),id) ) then 349 350 do i = 1, nshl 351 EGmasst(iel,i,inod) = zero 352 EGmasst(iel,inod,i) = zero 353 enddo 354 EGmasst(iel,inod,inod) = one 355 endif 356 357c 358 359 5000 continue 360 361c 362c.... end loop over shape functions (nodes) 363c 364 enddo 365c 366c.... end loop over elements 367c 368 enddo 369c 370c.... return 371c 372 return 373 end 374 375