Télécharger la plaquette Cast3M


fr en

Télécharger inclu4.eso

Retour à la liste

Numérotation des lignes : 

inclu4

  1. C INCLU4    SOURCE    PV        20/03/30    21:19:57     10567          
  2.       SUBROUTINE INCLU4(IPT1,IPT2,IPEX,XCRIT)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8(A-H,O-Z)
  5. -INC CCREEL
  6.  
  7. -INC PPARAM
  8. -INC CCOPTIO
  9. -INC SMELEME
  10. -INC SMCOORD
  11. *
  12. * IL EST SUPPOSER QUE LE IPT2 ne contiennent quer des TET4 et que ipt1
  13. * soit des poi1,
  14. *
  15. * IPEX(I)=1 veut dire que le noeud I est interne à un element de ipt2
  16. *
  17. * POUR DECIDER SI UN POINT EST A L'INTERIEUR D'UN ELEMENT ON CALCULE
  18. * LES COORDONNEES BARYCENTRIQUES DU POINT ET IL FAUT QU'ELLES SOIENT
  19. * TOUTES POSITIVES OU QUE CELLES QUI SOIENT NEGATIVES SOIENT D'UN
  20. * ORDRE DE GRANDEUR TRES INFERIEUR AUX AUTRES ( 0.0001 FOIS ).
  21. *
  22.       SEGMENT ISEG1
  23.          REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL)
  24.       ENDSEGMENT
  25.       SEGMENT ISEG3
  26.          INTEGER NIZO(NZO+1)
  27.       ENDSEGMENT
  28.       SEGMENT ISEG4
  29.          INTEGER NUMZO(NZO)
  30.       ENDSEGMENT
  31.       SEGMENT ISEG5
  32.         INTEGER NNMEL(ILON),IDEJ(NZO)
  33.       ENDSEGMENT
  34.       SEGMENT ISEG6
  35.       REAL*8 AM(4,4),P(4),AL(4),A(4),B(4),C(4),D(4)
  36.       REAL*8 XA(3,4),XPU(3)
  37.       ENDSEGMENT
  38.       SEGMENT IPEX(nbpts)
  39. *
  40. *  ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS
  41. * AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE
  42. * LA DIMENSION HORS TOUT DU MAILLAGE ET ON COMMENCE A CREER
  43. * UNE NUMEROTATION LOCALE DES POINTS DU MAILLAGE IPT
  44. *
  45.       IDIM1=IDIM+1
  46.       XPREC=-XCRIT
  47.       MELEME= IPT2
  48.       SEGACT MELEME
  49.       IF(ITYPEL.NE.23) THEN
  50.          CALL ERREUR(16)
  51.          RETURN
  52.       ENDIF
  53.       NBEL = NUM(/2)
  54.       NBNN=NUM(/1)
  55. *      WRITE(6,FMT='(''  NBEL NBNN '',2I6)') NBEL,NBNN
  56.       SEGINI ISEG1
  57.       ILOC=0
  58.       XZO=0.D0
  59.       YZO=0.D0
  60.       ZZO=0.D0
  61.       XZA=XGRAND
  62.       YZA=XGRAND
  63.       ZZA=XGRAND
  64.       XTOMI=XGRAND
  65.       XTOMA=-XGRAND
  66.       YTOMI=XGRAND
  67.       YTOMA=-XGRAND
  68.       ZTOMI=XGRAND
  69.       ZTOMA=-XGRAND
  70.       DO 1 I1=1,NBEL
  71.          XMI=XGRAND
  72.          YMI=XGRAND
  73.          ZMI=XGRAND
  74.          YMA=-XGRAND
  75.          XMA=-XGRAND
  76.          ZMA=-XGRAND
  77.          DO 2 I2 = 1,NBNN
  78.             IB=NUM(I2,I1)
  79.             IA=(IB-1)*IDIM1
  80.             IF(XCOOR(IA+1).LT.XMI) XMI=XCOOR(IA+1)
  81.             IF(XCOOR(IA+1).GT.XMA) XMA=XCOOR(IA+1)
  82.             IF(XCOOR(IA+2).LT.YMI) YMI=XCOOR(IA+2)
  83.             IF(XCOOR(IA+2).GT.YMA) YMA=XCOOR(IA+2)
  84. *            IF( IDIM.EQ.3 ) THEN
  85.                IF(XCOOR(IA+3).LT.ZMI) ZMI=XCOOR(IA+3)
  86.                IF(XCOOR(IA+3).GT.ZMA) ZMA=XCOOR(IA+3)
  87. *            ENDIF
  88.     2    CONTINUE
  89.          XLIM(1,I1)=XMI
  90.          XLIM(2,I1)=XMA
  91.          YLIM(1,I1)=YMI
  92.          YLIM(2,I1)=YMA
  93.          XZO=MAX (XZO,XMA-XMI)
  94.          YZO=MAX (YZO,YMA-YMI)
  95.          XZA=MIN(XZA,XMA-XMI)
  96.          YZA=MIN(YZA,YMA-YMI)
  97.          IF(XMI.LT.XTOMI) XTOMI=XMI
  98.          IF(XMA.GT.XTOMA) XTOMA=XMA
  99.          IF(YMI.LT.YTOMI) YTOMI=YMI
  100.          IF(YMA.GT.YTOMA) YTOMA=YMA
  101. *         IF(IDIM.EQ.3) THEN
  102.              ZLIM(1,I1)=ZMI
  103.              ZLIM(2,I1)=ZMA
  104.              ZZO=MAX(ZZO,ZMA-ZMI)
  105.              ZZA=MIN(ZZA,ZMA-ZMI)
  106.              IF(ZMI.LT.ZTOMI) ZTOMI=ZMI
  107.              IF(ZMA.GT.ZTOMA) ZTOMA=ZMA
  108. *         ENDIF
  109.     1 CONTINUE
  110. *      WRITE(6,FMT='(''XZO YZO  '',4E12.5)') XZO,YZO
  111. *      WRITE(6,FMT='(''XZA YZA  '',4E12.5)') XZA,YZA
  112. *      WRITE(6,FMT='(''XTOMI XTOMA  '',4E12.5)') XTOMI,XTOMA
  113. *      WRITE(6,FMT='(''YTOMI YTOMA  '',4E12.5)') YTOMI,YTOMA
  114.       XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4)
  115.       YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4)
  116. C      XZO=XZO*1.1
  117. C      YZO=YZO*1.1
  118.       XZA=XZA*0.97
  119.       YZA=YZA*0.97
  120.       XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4
  121.       XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4
  122.       YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4
  123.       YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4
  124. C      XZO=MIN ( XZO, XTOMA-XTOMI)
  125. C      YZO=MIN ( YZO, YTOMA-YTOMI)
  126.       XZA=MIN ( XZA, XTOMA-XTOMI)
  127.       YZA=MIN ( YZA, YTOMA-YTOMI)
  128.       NXZO=INT((XTOMA-XTOMI)/XZA) + 1
  129.       NYZO=INT((YTOMA-YTOMI)/YZA) + 1
  130.       XZO=XZA
  131.       YZO=YZA
  132.       NZZO=1
  133. *      WRITE(6,FMT='('' NXZO NYZO'',2I7)') NXZO,NYZO
  134. *      IF(IDIM.EQ.3) THEN
  135.          ZPR=MIN(ZZO*1.D-2,(ZTOMA-ZTOMI)/2.D+4)
  136. C         ZZO=ZZO*1.1
  137.          ZZA=ZZA*0.97
  138.          ZTOMI= ZTOMI - (ZTOMA-ZTOMI)/1.D+4
  139.          ZTOMA= ZTOMA + (ZTOMA-ZTOMI)/1.D+4
  140. C         ZZO=MIN ( ZZO, ZTOMA-ZTOMI)
  141.          ZZA=MIN ( ZZA, ZTOMA-ZTOMI)
  142.          NZZO=INT((ZTOMA-ZTOMI)/ZZA)+ 1
  143.          ZZO=ZZA
  144. *      WRITE(6,FMT='('' zz0,zzA,ztomi,ztoma'',4e12.5)')
  145. *     $ xzo,xza,ztomi,ztoma
  146. *      ENDIF
  147. *      WRITE(6,FMT='('' XTOMI XTOMA YTOMI YTOMA '',4E12.5 )')
  148. *     $ XTOMI, XTOMA, YTOMI ,YTOMA
  149.        NXDEP=MIN(NXZO,10)
  150.        NYDEP=MIN(NYZO,10)
  151. *      IF(IDIM.EQ.2) THEN
  152. *         IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  153. *          XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90
  154. *          NXZO=MAX(INT(NXZO/XY),NXDEP)
  155. *          NYZO=MAX(INT(NYZO/XY),NYDEP)
  156. *          IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  157. *            XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60
  158. *            NXZO=MAX(INT(NXZO/XY),NXDEP)
  159. *            NYZO=MAX(INT(NYZO/XY),NYDEP)
  160. *          ENDIF
  161. *          XZO=(XTOMA-XTOMI)/NXZO
  162. *          YZO=(YTOMA-YTOMI)/NYZO
  163. *          NXZO=(XTOMA-XTOMI)/XZO +1
  164. *          NYZO=(YTOMA-YTOMI)/YZO +1
  165. *         ENDIF
  166. C        WRITE(6,FMT='(''  XZO NXZO YZO NYZO '' , E12.5,I5,E12.5,I5)')
  167. C     $  XZO ,NXZO, YZO, NYZO
  168. *      ELSE
  169.        NZDEP=MIN(NZZO,10)
  170. *        WRITE(6,FMT='(''  XZO NXZO YZO NYZO ZZO NZZO'' , E12.5,I7,/,
  171. *     $   E12.5,I7,E12.5,I7)')
  172. C     $  XZO ,NXZO, YZO, NYZO,ZZO,NZZO
  173.       IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO  ''
  174.      $,4I7) ') NXZO,NYZO,NZZO
  175.          IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.25000.) THEN
  176.             XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25.
  177.             NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP)
  178.             NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP)
  179.             NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP)
  180.       IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO  ''
  181.      $,4I7) ') NXZO,NYZO,NZZO
  182.          IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.20000.) THEN
  183.             XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25.
  184.                NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP)
  185.                NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP)
  186.                NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP)
  187.       IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO  ''
  188.      $,4I7) ') NXZO,NYZO,NZZO
  189.          IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.20000.) THEN
  190.             XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25.
  191.                   NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP)
  192.                   NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP)
  193.                   NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP)
  194.                ENDIF
  195.             ENDIF
  196.             XZO=(XTOMA-XTOMI)/FLOAT(NXZO)
  197.             YZO=(YTOMA-YTOMI)/FLOAT(NYZO)
  198.             ZZO=(ZTOMA-ZTOMI)/FLOAT(NZZO)
  199.             NXZO=INT((XTOMA-XTOMI)/XZO)+1
  200.             NYZO=INT((YTOMA-YTOMI)/YZO)+1
  201.             NZZO=INT((ZTOMA-ZTOMI)/ZZO)+1
  202.          ENDIF
  203. *      ENDIF
  204. *
  205. *  ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE
  206. *  POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT
  207. *  CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES
  208. *  PAR CHAQUE ELEMENT ET LEUR NOMBRE
  209. *
  210.  
  211.       NZO=NXZO*NYZO*NZZO
  212.       IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO  ''
  213.      $,4I7) ') NZO,NXZO,NYZO,NZZO
  214.       NXYZO=NXZO*NYZO
  215. *      IDI=4
  216. *      IF(IDIM.EQ.3) THEN
  217. *        IDI=8
  218. *      ENDIF
  219.       SEGINI ISEG3
  220.       SEGINI ISEG4
  221.       DO 3 I1=1,NBEL
  222.           NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  223.           NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  224.           NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  225.           NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  226. *          IF(IDIM.EQ.3) THEN
  227.                NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1
  228.                NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1
  229.                DO 200 L3=NIZ1Z,NIZ2Z
  230.                DO 200 L1=NIZ1Y,NIZ2Y
  231.                DO 200 L2=NIZ1X,NIZ2X
  232.                    NIZA = L2 + ( L1-1) * NXZO + ( L3-1)*NXYZO
  233.                    NUMZO(NIZA) = NUMZO(NIZA) +1
  234.   200          CONTINUE
  235. *           ELSE
  236. *               DO 201 L1=NIZ1Y,NIZ2Y
  237. *               DO 201 L2=NIZ1X,NIZ2X
  238. *                   NIZA = L2 + ( L1-1) * NXZO
  239. *                   NUMZO(NIZA) = NUMZO(NIZA) +1
  240. *  201          CONTINUE
  241. *           ENDIF
  242.     3 CONTINUE
  243. *
  244. *  CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
  245. *  ELEMENTS CONCERNEES PAR UNE ZONE
  246. *
  247.       ILON=0
  248.       NIZO(1)=1
  249.       DO 202 L1=1,NZO
  250.       NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
  251.       ILON=ILON+ NUMZO(L1)
  252.  202  CONTINUE
  253. *      WRITE(6,FMT='(''  ILON '',I5)') ILON
  254. *      WRITE(6,109) (KKK,NUMZO(KKK),(NELZO(KI,KKK),KI=1,4),KKK=1,NBEL)
  255. *  109 FORMAT(I6,I5,4I5)
  256. *      WRITE(6,110)( NIZO(KI),KI=1,NZO+1)
  257.  110  FORMAT(16I5)
  258.       SEGINI ISEG5
  259.       DO 5 I1=1,NBEL
  260.           NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  261.           NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  262.           NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  263.           NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  264. *          IF(IDIM.EQ.3) THEN
  265.                NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1
  266.                NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1
  267.                DO 205 L3=NIZ1Z,NIZ2Z
  268.                DO 205 L1=NIZ1Y,NIZ2Y
  269.                DO 205 L2=NIZ1X,NIZ2X
  270.                    NIZA = L2 + ( L1-1) * NXZO + ( L3-1)*NXYZO
  271.                    IAD=NIZO(NIZA)+IDEJ(NIZA)
  272.                    NNMEL(IAD)=I1
  273.                    IDEJ(NIZA)=IDEJ(NIZA)+1
  274.   205          CONTINUE
  275. *           ELSE
  276. *               DO 203 L1=NIZ1Y,NIZ2Y
  277. *               DO 203 L2=NIZ1X,NIZ2X
  278. *                   NIZA = L2 + ( L1-1) * NXZO
  279. *                   IAD=NIZO(NIZA)+IDEJ(NIZA)
  280. *                   NNMEL(IAD)=I1
  281. *                   IDEJ(NIZA)=IDEJ(NIZA)+1
  282. *  203          CONTINUE
  283. *           ENDIF
  284.     5  CONTINUE
  285. *
  286. *  IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE
  287. * POINT DE L'OBJETR MAILLAGE IPT1, ON COMMENCE PAR LE METTRE SOUS
  288. * FORME D'ELEMENTS DE TYPE POI1
  289. *
  290.       SEGSUP ISEG1,ISEG4
  291.       SEGACT IPT1
  292.       IF(IPT1.ITYPEL.NE.1) THEN
  293.       CALL CHANGE(IPT1,1)
  294.       ENDIF
  295.       SEGINI IPEX
  296.       SEGINI ISEG6
  297. C      WRITE(6,FMT='(''  AVANT BOUCLE 10'')')
  298.       DO 10 I=1,IPT1.NUM(/2)
  299.         IP=IPT1.NUM(1,I)
  300.         XPU(1)=XCOOR((IP-1)*IDIM1+1)
  301.         XPU(2)=XCOOR((IP-1)*IDIM1+2)
  302.         XPU(3)=XCOOR((IP-1)*IDIM1+3)
  303. *        write(6,fmt='('' point x y '',i5,2e12.5)') iP,XPU(1),xpu(2)
  304.         IF(XPU(1).LT.XTOMI.OR.XPU(1).GT.XTOMA) GO TO 10
  305.         IF(XPU(2).LT.YTOMI.OR.XPU(2).GT.YTOMA) GO TO 10
  306. *        IF(IDIM.EQ.3) THEN
  307.           IF(XPU(3).LT.ZTOMI.OR.XPU(3).GT.ZTOMA) GO TO 10
  308. *        ENDIF
  309.         INDZO=INT((XPU(1)-XTOMI)/XZO)+ 1 +INT((XPU(2)-YTOMI)/YZO)*NXZO
  310.      #        +INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO
  311. *        IF(IDIM.EQ.3) INDZO=INDZO+INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO
  312.         IDEB=NIZO(INDZO)
  313.         IFIN=NIZO(INDZO+1)-1
  314. *        write(6,fmt='('' ideb ifin'',2i5)') ideb,ifin
  315.         IF(IDEB.GT.IFIN) GO TO 10
  316.         IEL=0
  317.         DO 11 KK=IDEB,IFIN
  318. *
  319. *  ON CALCULE LES COORDONNEES BARYCENTRIQUES ( AU PLUS 4 )
  320. *
  321.           K=NNMEL(KK)
  322. *          write(6,fmt='('' k '',i5)') k
  323.           J1=NUM(1,K)
  324.           J2=NUM(2,K)
  325.           J3=NUM(3,K)
  326.           J1IDIM=(J1-1)*IDIM1
  327.           J2IDIM=(J2-1)*IDIM1
  328.           J3IDIM=(J3-1)*IDIM1
  329. *          IF(IDIM.EQ.3) THEN
  330.              J4=NUM(4,K)
  331.              J4IDIM=(J4-1)*IDIM1
  332. *          ENDIF
  333.           P(1)=1.D0
  334.           DO 12 K1=1,NBNN
  335.             AM(1,K1)=1.D0
  336.    12     CONTINUE
  337.           DO 13 K1=1,IDIM
  338.             P(K1+1)=XPU(K1)
  339.             AM(K1+1,1)=XCOOR(J1IDIM+K1)
  340.             AM(K1+1,2)=XCOOR(J2IDIM+K1)
  341.             AM(K1+1,3)=XCOOR(J3IDIM+K1)
  342. *            IF(IDIM.EQ.3) THEN
  343.             AM(K1+1,4)=XCOOR(J4IDIM+K1)
  344. *            ENDIF
  345.    13     CONTINUE
  346. *          IF(IDIM.EQ.2) THEN
  347. *             X1=AM(2,1)
  348. *             X2=AM(2,2)
  349. *             X3=AM(2,3)
  350. *             Y1=AM(3,1)
  351. *             Y2=AM(3,2)
  352. *             Y3=AM(3,3)
  353. *             X=P(2)
  354. *             Y=P(3)
  355. *             DETAM=X1*Y2+X2*Y3+X3*Y1-Y1*X2-Y2*X3-Y3*X1
  356. *             A(1)=X2*Y3-X3*Y2
  357. *             A(2)=X3*Y1-X1*Y3
  358. *             A(3)=X1*Y2-X2*Y1
  359. *             B(1)=Y2-Y3
  360. *             B(2)=Y3-Y1
  361. *             B(3)=Y1-Y2
  362. *             C(1)=X3-X2
  363. *             C(2)=X1-X3
  364. *             C(3)=X2-X1
  365. *             DO 14 IK=1,NBNN
  366. *               AL(IK)=(A(IK)+B(IK)*X+C(IK)*Y)/DETAM
  367. *   14        CONTINUE
  368. *             AL(4)=1.D0
  369. *          ELSE
  370.              X1=AM(2,1)
  371.              X2=AM(2,2)
  372.              X3=AM(2,3)
  373.              X4=AM(2,4)
  374.              Y1=AM(3,1)
  375.              Y2=AM(3,2)
  376.              Y3=AM(3,3)
  377.              Y4=AM(3,4)
  378.              Z1=AM(4,1)
  379.              Z2=AM(4,2)
  380.              Z3=AM(4,3)
  381.              Z4=AM(4,4)
  382.              X=P(2)
  383.              Y=P(3)
  384.              Z=P(4)
  385.       DETAM=X2*Y3*Z4+X3*Y4*Z2+X4*Y2*Z3-X4*Y3*Z2-X2*Y4*Z3-X3*Y2*Z4-X1*Y3*
  386.      1Z4-X3*Y4*Z1-X4*Y1*Z3+X4*Y3*Z1+X3*Y1*Z4+X1*Y4*Z3+X1*Y2*Z4+X4*Y1*Z2+
  387.      2X2*Y4*Z1-X4*Y2*Z1-X2*Y1*Z4-X1*Y4*Z2-X1*Y2*Z3-X3*Y1*Z2-X2*Y3*Z1+X3*
  388.      3Y2*Z1+X2*Y1*Z3+X1*Y3*Z2
  389.              A(1)=X2*Y3*Z4+X3*Y4*Z2+X4*Y2*Z3-X4*Y3*Z2-X2*Y4*Z3-X3*Y2*Z4
  390.              A(2)=X4*Y3*Z1+X3*Y1*Z4+X1*Y4*Z3-X1*Y3*Z4-X3*Y4*Z1-X4*Y1*Z3
  391.              A(3)=X1*Y2*Z4+X4*Y1*Z2+X2*Y4*Z1-X4*Y2*Z1-X2*Y1*Z4-X1*Y4*Z2
  392.              A(4)=X3*Y2*Z1+X2*Y1*Z3+X1*Y3*Z2-X1*Y2*Z3-X3*Y1*Z2-X2*Y3*Z1
  393.              B(1)=Y4*Z3-Y3*Z4+Y2*Z4-Y4*Z2+Y3*Z2-Y2*Z3
  394.              B(2)=Y3*Z4-Y4*Z3+Y4*Z1-Y1*Z4+Y1*Z3-Y3*Z1
  395.              B(3)=Y4*Z2-Y2*Z4+Y1*Z4-Y4*Z1+Y2*Z1-Y1*Z2
  396.              B(4)=Y2*Z3-Y3*Z2+Y3*Z1-Y1*Z3+Y1*Z2-Y2*Z1
  397.              C(1)=X3*Z4-X4*Z3+X4*Z2-X2*Z4+X2*Z3-X3*Z2
  398.              C(2)=X4*Z3-X3*Z4+X1*Z4-X4*Z1+X3*Z1-X1*Z3
  399.              C(3)=X2*Z4-X4*Z2+X4*Z1-X1*Z4+X1*Z2-X2*Z1
  400.              C(4)=X3*Z2-X2*Z3+X1*Z3-X3*Z1+X2*Z1-X1*Z2
  401.              D(1)=X4*Y3-X3*Y4+X2*Y4-X4*Y2+X3*Y2-X2*Y3
  402.              D(2)=X3*Y4-X4*Y3+X4*Y1-X1*Y4+X1*Y3-X3*Y1
  403.              D(3)=X4*Y2-X2*Y4+X1*Y4-X4*Y1+X2*Y1-X1*Y2
  404.              D(4)=X2*Y3-X3*Y2+X3*Y1-X1*Y3+X1*Y2-X2*Y1
  405.              DO 15 II=1,NBNN
  406.               AL(II)=(A(II)+B(II)*X+C(II)*Y+D(II)*Z)/DETAM
  407.    15        CONTINUE
  408. *          ENDIF
  409. *         write(6,fmt='(''K  al '',I5,4e12.3)')K,
  410. *     $  al(1),al(2),al(3),al(4)
  411.           IF( AL(1).GT.XPREC.AND.AL(2).GT.XPREC.AND.AL(3).GT.XPREC.
  412.      $    AND.AL(4).GT.XPREC) THEN
  413. *
  414. *  LE POINT EST INTERNE A L'ELEMENT
  415. *
  416. *             IEL=IEL+1
  417.              IPEX(IP)=1
  418. *             WRITE(6,*) IP
  419.              GO TO 10
  420.           ENDIF
  421.    11   CONTINUE
  422.    10 CONTINUE
  423.       SEGSUP ISEG5,ISEG6,ISEG3
  424.       SEGDES IPT1,MELEME
  425.       RETURN
  426.       END
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales