1 subroutine e3StsLhs( xl, lStsVec ) 2c----------------------------------------------------------------------- 3c 4c compute the terms needed for the left hand side matrices 5c needed for the conservative projection 6c 7c----------------------------------------------------------------------- 8 use stats 9 include "common.h" 10 11 integer i 12 real*8 lDir(npro,nshl,3), lStsVec(npro,nshl,nResDims), 13 & xl(npro,nenl,3) 14 15 call e3StsDir( xl, lDir ) 16 17 do i = 1, nshl 18 lStsVec(:,i,1) = lDir(:,i,1) * lDir(:,i,1) 19 lStsVec(:,i,2) = lDir(:,i,2) * lDir(:,i,2) 20 lStsVec(:,i,3) = lDir(:,i,3) * lDir(:,i,3) 21 22 lStsVec(:,i,4) = lDir(:,i,1) * lDir(:,i,2) 23 lStsVec(:,i,5) = lDir(:,i,2) * lDir(:,i,3) 24 lStsVec(:,i,6) = lDir(:,i,3) * lDir(:,i,1) 25 26 lStsVec(:,i,7) = 0.0 27 lStsVec(:,i,8) = 0.0 28 lStsVec(:,i,9) = 0.0 29 lStsVec(:,i,10) = 0.0 30 lStsVec(:,i,11) = 0.0 31 enddo 32 33 return 34 end 35 36 37 subroutine e3StsRes( xl, rl, lStsVec ) 38c----------------------------------------------------------------------- 39c 40c compute the residual terms for the consistent projection 41c 42c----------------------------------------------------------------------- 43 use stats 44 include "common.h" 45 46 real*8 xl(npro,nenl,3), rl(npro,nshl,ndof) 47 real*8 lDir(npro,nshl,3), lStsVec(npro,nshl,nResDims) 48 49 call e3StsDir( xl, lDir ) 50 51 do i = 1, nshl 52 lStsVec(:,i,1) = lDir(:,i,1) * rl(:,i,4) 53 lStsVec(:,i,2) = lDir(:,i,2) * rl(:,i,4) 54 lStsVec(:,i,3) = lDir(:,i,3) * rl(:,i,4) 55 56 lStsVec(:,i,4) = lDir(:,i,1) * rl(:,i,1) 57 lStsVec(:,i,5) = lDir(:,i,2) * rl(:,i,2) 58 lStsVec(:,i,6) = lDir(:,i,3) * rl(:,i,3) 59 60 lStsVec(:,i,7) = lDir(:,i,1) * rl(:,i,2) 61 & + lDir(:,i,2) * rl(:,i,1) 62 lStsVec(:,i,8) = lDir(:,i,2) * rl(:,i,3) 63 & + lDir(:,i,3) * rl(:,i,2) 64 lStsVec(:,i,9) = lDir(:,i,3) * rl(:,i,1) 65 & + lDir(:,i,1) * rl(:,i,3) 66 lStsVec(:,i,10) = 0 67 lStsVec(:,i,11) = 0 68 enddo 69 70 71 return 72 end 73 74 subroutine e3StsDir( xl, lDir ) 75c----------------------------------------------------------------------- 76c 77c compute the normal to each of the nodes 78c 79c----------------------------------------------------------------------- 80 include "common.h" 81 82 real*8 xl(npro,nenl,3), lDir(npro,nshl,3) 83 integer e 84 85c 86c.... linear tets 87c 88 if (nshl .eq. 4 ) then 89 fct = 1.d0 / 6.d0 90 do e = 1, npro 91 92 x12 = xl(e,2,1) - xl(e,1,1) 93 x13 = xl(e,3,1) - xl(e,1,1) 94 x14 = xl(e,4,1) - xl(e,1,1) 95 x23 = xl(e,3,1) - xl(e,2,1) 96 x24 = xl(e,4,1) - xl(e,2,1) 97 x34 = xl(e,4,1) - xl(e,3,1) 98 99 y12 = xl(e,2,2) - xl(e,1,2) 100 y13 = xl(e,3,2) - xl(e,1,2) 101 y14 = xl(e,4,2) - xl(e,1,2) 102 y23 = xl(e,3,2) - xl(e,2,2) 103 y24 = xl(e,4,2) - xl(e,2,2) 104 y34 = xl(e,4,2) - xl(e,3,2) 105 106 z12 = xl(e,2,3) - xl(e,1,3) 107 z13 = xl(e,3,3) - xl(e,1,3) 108 z14 = xl(e,4,3) - xl(e,1,3) 109 z23 = xl(e,3,3) - xl(e,2,3) 110 z24 = xl(e,4,3) - xl(e,2,3) 111 z34 = xl(e,4,3) - xl(e,3,3) 112 113c 114c.. The calculation of the direction of a vertex is based on the average of 115c.. the normals of the neighbor faces(3); And the calculation of the direction 116c.. of the edge is based on the neighbor faces(2); 117c 118 lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14) 119 & + y13*(-z12 + z14)) 120 lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14) 121 & + x12*(-z13 + z14)) 122 lDir(e,1,3) = fct * ( x14*(y12 - y13) 123 & + x12*(y13 - y14) + x13*(-y12 + y14)) 124c 125 lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13 126 & - y12*z14 + y24*z23 - y23*z24) 127 lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14 128 & - x24*z23 + x23*z24 ) 129 lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13 130 & - x12*y14 + x24*y23 - x23*y24) 131c 132 lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14) 133 & + y24*z23 - y23*z24) 134 lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14) 135 & - x24*z23 + x23*z24) 136 lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14) 137 & + x24*y23 - x23*y24) 138c 139 lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14 140 & + y24*z23 - y23*z24) 141 lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14 142 & - x24*z23 + x23*z24) 143 lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14 144 & + x24*y23 - x23*y24) 145c 146 enddo 147c 148c.... quadratic tets 149c 150 else if (nshl .eq. 10 ) then 151 fct = 1.d0 / 6.d0 152 do e = 1, npro 153 154 x12 = xl(e,2,1) - xl(e,1,1) 155 x13 = xl(e,3,1) - xl(e,1,1) 156 x14 = xl(e,4,1) - xl(e,1,1) 157 x23 = xl(e,3,1) - xl(e,2,1) 158 x24 = xl(e,4,1) - xl(e,2,1) 159 x34 = xl(e,4,1) - xl(e,3,1) 160 161 y12 = xl(e,2,2) - xl(e,1,2) 162 y13 = xl(e,3,2) - xl(e,1,2) 163 y14 = xl(e,4,2) - xl(e,1,2) 164 y23 = xl(e,3,2) - xl(e,2,2) 165 y24 = xl(e,4,2) - xl(e,2,2) 166 y34 = xl(e,4,2) - xl(e,3,2) 167 168 z12 = xl(e,2,3) - xl(e,1,3) 169 z13 = xl(e,3,3) - xl(e,1,3) 170 z14 = xl(e,4,3) - xl(e,1,3) 171 z23 = xl(e,3,3) - xl(e,2,3) 172 z24 = xl(e,4,3) - xl(e,2,3) 173 z34 = xl(e,4,3) - xl(e,3,3) 174 175c 176c.... vertex modes 177c 178 lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14) 179 & + y13*(-z12 + z14)) 180 lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14) 181 & + x12*(-z13 + z14)) 182 lDir(e,1,3) = fct * ( x14*(y12 - y13) 183 & + x12*(y13 - y14) + x13*(-y12 + y14)) 184c 185 lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13 186 & - y12*z14 + y24*z23 - y23*z24) 187 lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14 188 & - x24*z23 + x23*z24 ) 189 lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13 190 & - x12*y14 + x24*y23 - x23*y24) 191c 192 lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14) 193 & + y24*z23 - y23*z24) 194 lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14) 195 & - x24*z23 + x23*z24) 196 lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14) 197 & + x24*y23 - x23*y24) 198c 199 lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14 200 & + y24*z23 - y23*z24) 201 lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14 202 & - x24*z23 + x23*z24) 203 lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14 204 & + x24*y23 - x23*y24) 205c 206c.... edge modes (quadratic) 207c 208 lDir(e,5,1) = pt25*(-(y13*z12) + y14*z12 + y12*(z13-z14)) 209 lDir(e,5,2) = pt25*(x13*z12 - x14*z12 + x12*(-z13 + z14)) 210 lDir(e,5,3) = pt25*(-(x13*y12) + x14*y12 + x12*(y13 - y14)) 211 212 lDir(e,6,1) = pt25*(-(y13*z12) + y12*z13 + y24*z23-y23*z24) 213 lDir(e,6,2) = pt25*(x13*z12 - x12*z13 - x24*z23 + x23*z24) 214 lDir(e,6,3) = pt25*(-(x13*y12) + x12*y13 + x24*y23-x23*y24) 215 216 lDir(e,7,1) = pt25*((y12 - y14)*z13 + y13*(-z12 + z14)) 217 lDir(e,7,2) = pt25*((-x12 + x14)*z13 + x13*(z12 - z14)) 218 lDir(e,7,3) = pt25*((x12 - x14)*y13 + x13*(-y12 + y14)) 219 220 lDir(e,8,1) = pt25*(y14*(z12 - z13) + (-y12 + y13)*z14) 221 lDir(e,8,2) = pt25*(x14*(-z12 + z13) + (x12 - x13)*z14) 222 lDir(e,8,3) = pt25*(x14*(y12 - y13) + (-x12 + x13)*y14) 223 224 lDir(e,9,1) = pt25*(y14*z12 - y12*z14 + y24*z23 - y23*z24) 225 lDir(e,9,2) = pt25*(-(x14*z12) + x12*z14 - x24*z23+x23*z24) 226 lDir(e,9,3) = pt25*(x14*y12 - x12*y14 + x24*y23 - x23*y24) 227 228 lDir(e,10,1) = pt25*(-(y14*z13) + y13*z14+y24*z23-y23*z24) 229 lDir(e,10,2) = pt25*(x14*z13 - x13*z14-x24*z23 + x23*z24) 230 lDir(e,10,3) = pt25*(-(x14*y13) + x13*y14+x24*y23-x23*y24) 231 232 233 enddo 234c 235c.... cubic tets 236c 237 else if (nshl .eq. 20 ) then 238 fct = 1.d0 / 6.d0 239 do e = 1, npro 240 241 x12 = xl(e,2,1) - xl(e,1,1) 242 x13 = xl(e,3,1) - xl(e,1,1) 243 x14 = xl(e,4,1) - xl(e,1,1) 244 x23 = xl(e,3,1) - xl(e,2,1) 245 x24 = xl(e,4,1) - xl(e,2,1) 246c$$$ x34 = xl(e,4,1) - xl(e,3,1) 247 248 y12 = xl(e,2,2) - xl(e,1,2) 249 y13 = xl(e,3,2) - xl(e,1,2) 250 y14 = xl(e,4,2) - xl(e,1,2) 251 y23 = xl(e,3,2) - xl(e,2,2) 252 y24 = xl(e,4,2) - xl(e,2,2) 253c$$$ y34 = xl(e,4,2) - xl(e,3,2) 254 255 z12 = xl(e,2,3) - xl(e,1,3) 256 z13 = xl(e,3,3) - xl(e,1,3) 257 z14 = xl(e,4,3) - xl(e,1,3) 258 z23 = xl(e,3,3) - xl(e,2,3) 259 z24 = xl(e,4,3) - xl(e,2,3) 260c$$$ z34 = xl(e,4,3) - xl(e,3,3) 261 262c 263c.... vertex modes 264c 265 lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14) 266 & + y13*(-z12 + z14)) 267 lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14) 268 & + x12*(-z13 + z14)) 269 lDir(e,1,3) = fct * ( x14*(y12 - y13) 270 & + x12*(y13 - y14) + x13*(-y12 + y14)) 271c 272 lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13 273 & - y12*z14 + y24*z23 - y23*z24) 274 lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14 275 & - x24*z23 + x23*z24 ) 276 lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13 277 & - x12*y14 + x24*y23 - x23*y24) 278c 279 lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14) 280 & + y24*z23 - y23*z24) 281 lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14) 282 & - x24*z23 + x23*z24) 283 lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14) 284 & + x24*y23 - x23*y24) 285c 286 lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14 287 & + y24*z23 - y23*z24) 288 lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14 289 & - x24*z23 + x23*z24) 290 lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14 291 & + x24*y23 - x23*y24) 292c 293c.... edge modes (quadratic and cubic) 294c 295 lDir(e,5,1) = pt25*(-(y13*z12) + y14*z12 + y12*(z13-z14)) 296 lDir(e,5,2) = pt25*(x13*z12 - x14*z12 + x12*(-z13 + z14)) 297 lDir(e,5,3) = pt25*(-(x13*y12) + x14*y12 + x12*(y13 - y14)) 298 lDir(e,6,1) = lDir(e,5,1) 299 lDir(e,6,2) = lDir(e,5,2) 300 lDir(e,6,3) = lDir(e,5,3) 301 302 lDir(e,7,1) = pt25*(-(y13*z12) + y12*z13 + y24*z23-y23*z24) 303 lDir(e,7,2) = pt25*(x13*z12 - x12*z13 - x24*z23 + x23*z24) 304 lDir(e,7,3) = pt25*(-(x13*y12) + x12*y13 + x24*y23-x23*y24) 305 lDir(e,8,1) = lDir(e,7,1) 306 lDir(e,8,2) = lDir(e,7,2) 307 lDir(e,8,3) = lDir(e,7,3) 308 309 lDir(e,9,1) = pt25*((y12 - y14)*z13 + y13*(-z12 + z14)) 310 lDir(e,9,2) = pt25*((-x12 + x14)*z13 + x13*(z12 - z14)) 311 lDir(e,9,3) = pt25*((x12 - x14)*y13 + x13*(-y12 + y14)) 312 lDir(e,10,1) = lDir(e,9,1) 313 lDir(e,10,2) = lDir(e,9,2) 314 lDir(e,10,3) = lDir(e,9,3) 315 316 lDir(e,11,1) = pt25*(y14*(z12 - z13) + (-y12 + y13)*z14) 317 lDir(e,11,2) = pt25*(x14*(-z12 + z13) + (x12 - x13)*z14) 318 lDir(e,11,3) = pt25*(x14*(y12 - y13) + (-x12 + x13)*y14) 319 lDir(e,12,1) = lDir(e,11,1) 320 lDir(e,12,2) = lDir(e,11,2) 321 lDir(e,12,3) = lDir(e,11,3) 322 323 324 lDir(e,13,1) = pt25*(y14*z12 - y12*z14 + y24*z23 - y23*z24) 325 lDir(e,13,2) = pt25*(-(x14*z12) + x12*z14-x24*z23+x23*z24) 326 lDir(e,13,3) = pt25*(x14*y12 - x12*y14 + x24*y23 - x23*y24) 327 lDir(e,14,1) = lDir(e,13,1) 328 lDir(e,14,2) = lDir(e,13,2) 329 lDir(e,14,3) = lDir(e,13,3) 330 331 lDir(e,15,1) = pt25*(-(y14*z13) + y13*z14+y24*z23-y23*z24) 332 lDir(e,15,2) = pt25*(x14*z13 - x13*z14-x24*z23 + x23*z24) 333 lDir(e,15,3) = pt25*(-(x14*y13) + x13*y14+x24*y23-x23*y24) 334 lDir(e,16,1) = lDir(e,15,1) 335 lDir(e,16,2) = lDir(e,15,2) 336 lDir(e,16,3) = lDir(e,15,3) 337c 338c.... face modes (cubic) 339c 340 lDir(e,17,1) = pt5*(-(y13*z12) + y12*z13) 341 lDir(e,17,2) = pt5*(x13*z12 - x12*z13) 342 lDir(e,17,3) = pt5*(-(x13*y12) + x12*y13) 343 344 lDir(e,18,1) = pt5*(y14*z12 - y12*z14) 345 lDir(e,18,2) = pt5*(-(x14*z12) + x12*z14) 346 lDir(e,18,3) = pt5*(x14*y12 - x12*y14) 347 348 lDir(e,19,1) = pt5*(y24*z23 - y23*z24) 349 lDir(e,19,2) = pt5*(-(x24*z23) + x23*z24) 350 lDir(e,19,3) = pt5*(x24*y23 - x23*y24) 351 352 lDir(e,20,1) = pt5*(-(y14*z13) + y13*z14) 353 lDir(e,20,2) = pt5*(x14*z13 - x13*z14) 354 lDir(e,20,3) = pt5*(-(x14*y13) + x13*y14) 355 356 enddo 357c 358c.... hexes 359c 360 else if (nenl .eq. 8) then 361 fct = 1.d0 / 12.d0 362 do e = 1, npro 363 x13 = xl(e,1,1) - xl(e,3,1) 364 x16 = xl(e,1,1) - xl(e,6,1) 365 x18 = xl(e,1,1) - xl(e,8,1) 366 x24 = xl(e,2,1) - xl(e,4,1) 367 x25 = xl(e,2,1) - xl(e,5,1) 368 x27 = xl(e,2,1) - xl(e,7,1) 369 x36 = xl(e,3,1) - xl(e,6,1) 370 x38 = xl(e,3,1) - xl(e,8,1) 371 x45 = xl(e,4,1) - xl(e,5,1) 372 x47 = xl(e,4,1) - xl(e,7,1) 373 x57 = xl(e,5,1) - xl(e,7,1) 374 x68 = xl(e,6,1) - xl(e,8,1) 375c 376 y13 = xl(e,1,2) - xl(e,3,2) 377 y16 = xl(e,1,2) - xl(e,6,2) 378 y18 = xl(e,1,2) - xl(e,8,2) 379 y24 = xl(e,2,2) - xl(e,4,2) 380 y25 = xl(e,2,2) - xl(e,5,2) 381 y27 = xl(e,2,2) - xl(e,7,2) 382 y36 = xl(e,3,2) - xl(e,6,2) 383 y38 = xl(e,3,2) - xl(e,8,2) 384 y45 = xl(e,4,2) - xl(e,5,2) 385 y47 = xl(e,4,2) - xl(e,7,2) 386 y57 = xl(e,5,2) - xl(e,7,2) 387 y68 = xl(e,6,2) - xl(e,8,2) 388c 389 z13 = xl(e,1,3) - xl(e,3,3) 390 z16 = xl(e,1,3) - xl(e,6,3) 391 z18 = xl(e,1,3) - xl(e,8,3) 392 z24 = xl(e,2,3) - xl(e,4,3) 393 z25 = xl(e,2,3) - xl(e,5,3) 394 z27 = xl(e,2,3) - xl(e,7,3) 395 z36 = xl(e,3,3) - xl(e,6,3) 396 z38 = xl(e,3,3) - xl(e,8,3) 397 z45 = xl(e,4,3) - xl(e,5,3) 398 z47 = xl(e,4,3) - xl(e,7,3) 399 z57 = xl(e,5,3) - xl(e,7,3) 400 z68 = xl(e,6,3) - xl(e,8,3) 401c 402 x31= -x13 403 x61= -x16 404 x81= -x18 405 x42= -x24 406 x52= -x25 407 x72= -x27 408 x63= -x36 409 x83= -x38 410 x54= -x45 411 x74= -x47 412 x75= -x57 413 x86= -x68 414 y31= -y13 415 y61= -y16 416 y81= -y18 417 y42= -y24 418 y52= -y25 419 y72= -y27 420 y63= -y36 421 y83= -y38 422 y54= -y45 423 y74= -y47 424 y75= -y57 425 y86= -y68 426 z31= -z13 427 z61= -z16 428 z81= -z18 429 z42= -z24 430 z52= -z25 431 z72= -z27 432 z63= -z36 433 z83= -z38 434 z54= -z45 435 z74= -z47 436 z75= -z57 437 z86= -z68 438 439 lDir(e,1,1) = fct * (-y24 * z45 + y36 * z24 - y68 * z45 440 1 + z24 * y45 - z36 * y24 + z68 * y45 ) 441 lDir(e,2,1) = fct * (-y16 * z63 + y54 * z16 - y47 * z63 442 1 + z16 * y63 - z54 * y16 + z47 * y63 ) 443 lDir(e,3,1) = fct * (-y42 * z27 + y18 * z42 - y86 * z27 444 1 + z42 * y27 - z18 * y42 + z86 * y27 ) 445 lDir(e,4,1) = fct * (-y38 * z81 + y72 * z38 - y25 * z81 446 1 + z38 * y81 - z72 * y38 + z25 * y81 ) 447 lDir(e,5,1) = fct * (-y61 * z18 + y27 * z61 - y74 * z18 448 1 + z61 * y18 - z27 * y61 + z74 * y18 ) 449 lDir(e,6,1) = fct * (-y57 * z72 + y81 * z57 - y13 * z72 450 1 + z57 * y72 - z81 * y57 + z13 * y72 ) 451 lDir(e,7,1) = fct * (-y83 * z36 + y45 * z83 - y52 * z36 452 1 + z83 * y36 - z45 * y83 + z52 * y36 ) 453 lDir(e,8,1) = fct * (-y75 * z54 + y63 * z75 - y31 * z54 454 1 + z75 * y54 - z63 * y75 + z31 * y54 ) 455c 456 lDir(e,1,2) = fct * (-z24 * x45 + z36 * x24 - z68 * x45 457 1 + x24 * z45 - x36 * z24 + x68 * z45 ) 458 lDir(e,2,2) = fct * (-z16 * x63 + z54 * x16 - z47 * x63 459 1 + x16 * z63 - x54 * z16 + x47 * z63 ) 460 lDir(e,3,2) = fct * (-z42 * x27 + z18 * x42 - z86 * x27 461 1 + x42 * z27 - x18 * z42 + x86 * z27 ) 462 lDir(e,4,2) = fct * (-z38 * x81 + z72 * x38 - z25 * x81 463 1 + x38 * z81 - x72 * z38 + x25 * z81 ) 464 lDir(e,5,2) = fct * (-z61 * x18 + z27 * x61 - z74 * x18 465 1 + x61 * z18 - x27 * z61 + x74 * z18 ) 466 lDir(e,6,2) = fct * (-z57 * x72 + z81 * x57 - z13 * x72 467 1 + x57 * z72 - x81 * z57 + x13 * z72 ) 468 lDir(e,7,2) = fct * (-z83 * x36 + z45 * x83 - z52 * x36 469 1 + x83 * z36 - x45 * z83 + x52 * z36 ) 470 lDir(e,8,2) = fct * (-z75 * x54 + z63 * x75 - z31 * x54 471 1 + x75 * z54 - x63 * z75 + x31 * z54 ) 472c 473 lDir(e,1,3) = fct * (-x24 * y45 + x36 * y24 - x68 * y45 474 1 + y24 * x45 - y36 * x24 + y68 * x45 ) 475 lDir(e,2,3) = fct * (-x16 * y63 + x54 * y16 - x47 * y63 476 1 + y16 * x63 - y54 * x16 + y47 * x63 ) 477 lDir(e,3,3) = fct * (-x42 * y27 + x18 * y42 - x86 * y27 478 1 + y42 * x27 - y18 * x42 + y86 * x27 ) 479 lDir(e,4,3) = fct * (-x38 * y81 + x72 * y38 - x25 * y81 480 1 + y38 * x81 - y72 * x38 + y25 * x81 ) 481 lDir(e,5,3) = fct * (-x61 * y18 + x27 * y61 - x74 * y18 482 1 + y61 * x18 - y27 * x61 + y74 * x18 ) 483 lDir(e,6,3) = fct * (-x57 * y72 + x81 * y57 - x13 * y72 484 1 + y57 * x72 - y81 * x57 + y13 * x72 ) 485 lDir(e,7,3) = fct * (-x83 * y36 + x45 * y83 - x52 * y36 486 1 + y83 * x36 - y45 * x83 + y52 * x36 ) 487 lDir(e,8,3) = fct * (-x75 * y54 + x63 * y75 - x31 * y54 488 1 + y75 * x54 - y63 * x75 + y31 * x54 ) 489c 490 enddo 491 else 492 write(*,*) 'Error in e3sts: elt type not impl.' 493 stop 494 endif 495 496 return 497 end 498 499 500 501