5 &( t,xm,ppq,lego,xmul,sw,w,h,
6 & surfac,ikle,nelem,nelmax)
85 INTEGER,
INTENT(IN) :: NELEM,NELMAX
86 INTEGER,
INTENT(IN) :: IKLE(nelmax,6),PPQ(6,6)
88 DOUBLE PRECISION,
INTENT(INOUT) :: T(nelmax,6),XM(30,nelmax)
90 DOUBLE PRECISION,
INTENT(IN) :: XMUL
91 DOUBLE PRECISION,
INTENT(IN) :: W(*)
92 DOUBLE PRECISION,
INTENT(IN) :: H(nelmax,6)
94 LOGICAL,
INTENT(IN) :: LEGO
98 TYPE(bief_obj),
INTENT(IN) :: SW
100 DOUBLE PRECISION,
INTENT(IN) :: SURFAC(nelmax)
107 DOUBLE PRECISION XSUR3
108 DOUBLE PRECISION LI0J0,LI0K0,LJ0I0,LJ0K0,LK0I0,LK0J0
109 DOUBLE PRECISION LI3J3,LI3K3,LJ3I3,LJ3K3,LK3I3,LK3J3
110 DOUBLE PRECISION LI3J0,LI3K0,LJ3I0,LJ3K0,LK3I0,LK3J0
111 DOUBLE PRECISION LI3I0,LJ3J0,LK3K0,LI0I3
112 DOUBLE PRECISION ALFA,ALFAJ,ALFAK,ALFA1,ALFA2,BETA,SOM0,SOM3
113 DOUBLE PRECISION ALFAI0,ALFAJ0,ALFAK0,ALFAI3,ALFAJ3,ALFAK3
114 DOUBLE PRECISION ALFAII,ALFAIJ,ALFAIK,ALFAJI,ALFAJJ,ALFAJK
115 DOUBLE PRECISION ALFAKI,ALFAKJ,ALFAKK
116 DOUBLE PRECISION,
PARAMETER :: EPS = 1.d-10
118 INTEGER :: IPLUS1(6),IPLUS3(6),INDIC(0:7)
119 INTEGER IXM,I0,J0,K0,I3,J3,K3
125 parameter( iplus1 = (/ 2 , 3 , 1 , 5 , 6 , 4 /) )
126 parameter( iplus3 = (/ 4 , 5 , 6 , 1 , 2 , 3 /) )
127 parameter( indic = (/ 4 , 4 , 5 , 3 , 6 , 2 , 1 , 1 /) )
166 WRITE(
lu,*)
'MT14PP DISCRETISATION NOT HANDLED' 177 IF (w(ikle(ielem,1)).GT.0.d0) ixm = 1
178 IF (w(ikle(ielem,2)).GT.0.d0) ixm = ixm + 2
179 IF (w(ikle(ielem,3)).GT.0.d0) ixm = ixm + 4
201 alfa = xsur3 * surfac(ielem)
203 li3i0 = alfa * abs(w(ikle(ielem,min(i0,i3))))
205 IF (ixm.GE.1.AND.ixm.LE.6)
THEN 209 lj3j0 = alfa * abs(w(ikle(ielem,min(j0,j3))))
210 xm(ppq(j0,j3),ielem) = 0.d0
211 lk3k0 = alfa * abs(w(ikle(ielem,min(k0,k3))))
212 xm(ppq(k0,k3),ielem) = 0.d0
214 alfai0 = -h(ielem,i0)
215 alfaj0 = -h(ielem,j0)
216 alfak0 = -h(ielem,k0)
217 alfai3 = -h(ielem,i3)
218 alfaj3 = -h(ielem,j3)
219 alfak3 = -h(ielem,k3)
221 lj0i0 = max(0.d0,min(alfai0,-alfaj0))
222 lk0i0 = max(0.d0,min(alfai0,-alfak0))
223 li3j3 = max(0.d0,min(alfaj3,-alfai3))
224 li3k3 = max(0.d0,min(alfak3,-alfai3))
226 alfaj = min(lj3j0,max(li3j3,lj0i0))
227 alfak = min(lk3k0,max(li3k3,lk0i0))
228 alfa = min(li0i3,alfaj+alfak)
229 IF (alfa.GT.eps)
THEN 231 beta = alfa / (alfaj+alfak)
241 alfaj0 = alfaj0+alfaj
242 alfaj3 = alfaj3-alfaj
243 alfak0 = alfak0+alfak
244 alfak3 = alfak3-alfak
246 lj0i0 = max(0.d0,min(alfai0,-alfaj0))
247 lk0i0 = max(0.d0,min(alfai0,-alfak0))
248 li3j3 = max(0.d0,min(alfaj3,-alfai3))
249 li3k3 = max(0.d0,min(alfak3,-alfai3))
253 li0j0 = max(0.d0,min(alfaj0,-alfai0))
254 lk0j0 = max(0.d0,min(alfaj0,-alfak0))
255 li0k0 = max(0.d0,min(alfak0,-alfai0))
256 lj0k0 = max(0.d0,min(alfak0,-alfaj0))
257 lj3i3 = max(0.d0,min(alfai3,-alfaj3))
258 lj3k3 = max(0.d0,min(alfak3,-alfaj3))
259 lk3i3 = max(0.d0,min(alfai3,-alfak3))
260 lk3j3 = max(0.d0,min(alfaj3,-alfak3))
262 xm(ppq(j0,i3),ielem) = 0.d0
263 xm(ppq(k0,i3),ielem) = 0.d0
264 xm(ppq(i0,j3),ielem) = 0.d0
265 xm(ppq(k0,j3),ielem) = 0.d0
266 xm(ppq(i0,k3),ielem) = 0.d0
267 xm(ppq(j0,k3),ielem) = 0.d0
279 alfa1 = min(li0i3,som0)
280 IF (alfa1.GT.eps)
THEN 283 alfaj0 = lj0i0 * beta
284 alfak0 = lk0i0 * beta
285 alfa = max(0.d0,alfa1-som3)
287 lj0i0 = lj0i0 - alfaj0*alfa1
288 lk0i0 = lk0i0 - alfak0*alfa1
289 xm(ppq(j0,i3),ielem) = alfaj0*alfa
290 xm(ppq(k0,i3),ielem) = alfak0*alfa
294 alfa2 = min(li0i3,som3)
295 IF (alfa2.GT.eps)
THEN 298 alfaj3 = li3j3 * beta
299 alfak3 = li3k3 * beta
300 alfa = max(0.d0,alfa2-som0)
302 li3j3 = li3j3 - alfaj3*alfa2
303 li3k3 = li3k3 - alfak3*alfa2
304 xm(ppq(i0,j3),ielem) = alfaj3*alfa
305 xm(ppq(i0,k3),ielem) = alfak3*alfa
309 alfa = min(alfa1,alfa2)
310 IF (alfa.GT.0.d0)
THEN 312 alfaj3 = alfaj3 * alfa
313 alfak3 = alfak3 * alfa
314 alfajj = alfaj3 * alfaj0
315 alfakk = alfak3 * alfak0
316 alfajk = alfaj3 * alfak0
317 alfakj = alfak3 * alfaj0
318 alfa = min(alfajk,alfakj)
320 xm(ppq(j0,j3),ielem) = alfajj + alfa
321 xm(ppq(k0,k3),ielem) = alfakk + alfa
322 xm(ppq(k0,j3),ielem) = alfajk - alfa
323 xm(ppq(j0,k3),ielem) = alfakj - alfa
327 xm(ppq(i0,i3),ielem) = max(0.d0,li0i3-max(som3,som0))
339 alfa1 = min(li3i0,som0)
340 IF (alfa1.GT.eps)
THEN 343 alfaj0 = li0j0 * beta
344 alfak0 = li0k0 * beta
345 alfa = max(0.d0,alfa1-som3)
347 li0j0 = li0j0 - alfaj0*alfa1
348 li0k0 = li0k0 - alfak0*alfa1
349 li3j0 = li3j0 + alfaj0*alfa
350 li3k0 = li3k0 + alfak0*alfa
354 xm(ppq(i0,j0),ielem) = li0j0
355 xm(ppq(i0,k0),ielem) = li0k0
357 alfa2 = min(li3i0,som3)
358 IF (alfa2.GT.eps)
THEN 361 alfaj3 = lj3i3 * beta
362 alfak3 = lk3i3 * beta
363 alfa = max(0.d0,alfa2-som0)
365 lj3i3 = lj3i3 - alfaj3*alfa2
366 lk3i3 = lk3i3 - alfak3*alfa2
367 lj3i0 = lj3i0 + alfaj3*alfa
368 lk3i0 = lk3i0 + alfak3*alfa
372 xm(ppq(j3,i3),ielem) = lj3i3
373 xm(ppq(k3,i3),ielem) = lk3i3
375 alfa = min(alfa1,alfa2)
376 IF (alfa.GT.0.d0)
THEN 378 alfaj0 = alfaj0 * alfa
379 alfak0 = alfak0 * alfa
380 alfajj = alfaj0 * alfaj3
381 alfakk = alfak0 * alfak3
382 alfajk = alfaj0 * alfak3
383 alfakj = alfak0 * alfaj3
384 alfa = min(alfajk,alfakj)
386 lj3j0 = lj3j0 + alfajj + alfa
387 lk3k0 = lk3k0 + alfakk + alfa
388 lk3j0 = lk3j0 + alfajk - alfa
389 lj3k0 = lj3k0 + alfakj - alfa
393 li3i0 = max(0.d0,li3i0-max(som0,som3))
398 alfa1 = min(lj3j0,som0)
399 IF (alfa1.GT.eps)
THEN 402 alfai0 = lj0i0 * beta
403 alfak0 = lj0k0 * beta
404 alfa = max(0.d0,alfa1-som3)
406 lj0i0 = lj0i0 - alfai0*alfa1
407 lj0k0 = lj0k0 - alfak0*alfa1
408 lj3i0 = lj3i0 + alfai0*alfa
409 lj3k0 = lj3k0 + alfak0*alfa
413 xm(ppq(j0,i0),ielem) = lj0i0
414 xm(ppq(j0,k0),ielem) = lj0k0
416 alfa2 = min(lj3j0,som3)
417 IF (alfa2.GT.eps)
THEN 420 alfai3 = li3j3 * beta
421 alfak3 = lk3j3 * beta
422 alfa = max(0.d0,alfa2-som0)
424 li3j3 = li3j3 - alfai3*alfa2
425 lk3j3 = lk3j3 - alfak3*alfa2
426 li3j0 = li3j0 + alfai3*alfa
427 lk3j0 = lk3j0 + alfak3*alfa
431 xm(ppq(i3,j3),ielem) = li3j3
432 xm(ppq(k3,j3),ielem) = lk3j3
434 alfa = min(alfa1,alfa2)
435 IF (alfa.GT.0.d0)
THEN 437 alfai0 = alfai0 * alfa
438 alfak0 = alfak0 * alfa
439 alfaii = alfai0 * alfai3
440 alfakk = alfak0 * alfak3
441 alfaik = alfai0 * alfak3
442 alfaki = alfak0 * alfai3
443 alfa = min(alfaik,alfaki)
445 li3i0 = li3i0 + alfaii + alfa
446 lk3k0 = lk3k0 + alfakk + alfa
447 lk3i0 = lk3i0 + alfaik - alfa
448 li3k0 = li3k0 + alfaki - alfa
452 lj3j0 = max(0.d0,lj3j0-max(som0,som3))
457 alfa1 = min(lk3k0,som0)
458 IF (alfa1.GT.eps)
THEN 461 alfai0 = lk0i0 * beta
462 alfaj0 = lk0j0 * beta
463 alfa = max(0.d0,alfa1-som3)
465 lk0i0 = lk0i0 - alfai0*alfa1
466 lk0j0 = lk0j0 - alfaj0*alfa1
467 lk3i0 = lk3i0 + alfai0*alfa
468 lk3j0 = lk3j0 + alfaj0*alfa
472 xm(ppq(k0,i0),ielem) = lk0i0
473 xm(ppq(k0,j0),ielem) = lk0j0
475 alfa2 = min(lk3k0,som3)
476 IF (alfa2.GT.eps)
THEN 479 alfai3 = li3k3 * beta
480 alfaj3 = lj3k3 * beta
481 alfa = max(0.d0,alfa2-som0)
483 li3k3 = li3k3 - alfai3*alfa2
484 lj3k3 = lj3k3 - alfaj3*alfa2
485 li3k0 = li3k0 + alfai3*alfa
486 lj3k0 = lj3k0 + alfaj3*alfa
490 xm(ppq(i3,k3),ielem) = li3k3
491 xm(ppq(j3,k3),ielem) = lj3k3
493 alfa = min(alfa1,alfa2)
494 IF (alfa.GT.0.d0)
THEN 496 alfai0 = alfai0 * alfa
497 alfaj0 = alfaj0 * alfa
498 alfaii = alfai0 * alfai3
499 alfajj = alfaj0 * alfaj3
500 alfaij = alfai0 * alfaj3
501 alfaji = alfaj0 * alfai3
502 alfa = min(alfaij,alfaji)
504 li3i0 = li3i0 + alfaii + alfa
505 lj3j0 = lj3j0 + alfajj + alfa
506 lj3i0 = lj3i0 + alfaij - alfa
507 li3j0 = li3j0 + alfaji - alfa
511 lk3k0 = max(0.d0,lk3k0-max(som0,som3))
513 xm(ppq(i3,i0),ielem) = li3i0
514 xm(ppq(j3,i0),ielem) = lj3i0
515 xm(ppq(k3,i0),ielem) = lk3i0
516 xm(ppq(i3,j0),ielem) = li3j0
517 xm(ppq(j3,j0),ielem) = lj3j0
518 xm(ppq(k3,j0),ielem) = lk3j0
519 xm(ppq(i3,k0),ielem) = li3k0
520 xm(ppq(j3,k0),ielem) = lj3k0
521 xm(ppq(k3,k0),ielem) = lk3k0
534 t(ielem,1) = -xm(01,ielem)-xm(02,ielem)
535 & -xm(03,ielem)-xm(04,ielem)-xm(05,ielem)
536 t(ielem,2) = -xm(16,ielem)-xm(06,ielem)
537 & -xm(07,ielem)-xm(08,ielem)-xm(09,ielem)
538 t(ielem,3) = -xm(17,ielem)-xm(21,ielem)
539 & -xm(10,ielem)-xm(11,ielem)-xm(12,ielem)
540 t(ielem,4) = -xm(18,ielem)-xm(22,ielem)
541 & -xm(25,ielem)-xm(13,ielem)-xm(14,ielem)
542 t(ielem,5) = -xm(19,ielem)-xm(23,ielem)
543 & -xm(26,ielem)-xm(28,ielem)-xm(15,ielem)
544 t(ielem,6) = -xm(20,ielem)-xm(24,ielem)
545 & -xm(27,ielem)-xm(29,ielem)-xm(30,ielem)
subroutine mt14pp(T, XM, PPQ, LEGO, XMUL, SW, W, H, SURFAC, IKLE, NELEM, NELMAX)