Télécharger la plaquette Cast3M


fr en

Télécharger biosav.eso

Retour à la liste

Numérotation des lignes : 

biosav

  1. C BIOSAV    SOURCE    FD218221  24/02/09    21:15:01     11837          
  2.       SUBROUTINE BIOSAV
  3. *----------------------------------------------------------
  4. *
  5. *           CALCUL DES CHAMPS DE BIOT ET SAVART
  6. *
  7. *----------------------------------------------------------
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8(A-H,O-Z)
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC CCREEL
  14. -INC SMELEME
  15. -INC SMCOORD
  16. -INC TMTRAV
  17. -INC SMCHPOI
  18.       INTEGER ENT
  19.       DIMENSION C(3),PN1(3),RESP(3)
  20.       DIMENSION VX(3),VZ(3),VY(3),VI(3),PP3(3)
  21.       DIMENSION ZF1(1),ZALF(1),ZDENS(1),ZBET(1)
  22.       DIMENSION P1(3),P2(3),P3(3),ROTA(3,3)
  23.       SEGMENT ICPR(nbpts)
  24.       SEGMENT SANGLE
  25.         REAL*8 TETM(NT4)
  26.         REAL*8 TETI(NT4)
  27.         REAL*8 DTE(NT4)
  28.       ENDSEGMENT
  29. *
  30.       CHARACTER*4 MCLE(4),CHSECT(1),CHOIX(2)
  31.       LOGICAL INDU,POTE
  32. *
  33.       DATA EPSIL3/1.D-12/,EPSIL2/1.D-8/
  34. *
  35.       DATA MCLE/'CERC','ARC ','BARR','FIL '/
  36.       DATA CHSECT/'TRAP'/
  37.       DATA CHOIX/'INDU','POTE'/
  38. *
  39.       segact mcoord
  40.       CALL LIROBJ('CHPOINT  ',MCHPO1,0,IRETOU)
  41.       IF (IERR.NE.0) RETURN
  42.       IF (IRETOU.EQ.1) THEN
  43.         CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  44.         IF (IERR.NE.0) RETURN
  45.  
  46. * OPTION CALCUL PAR INTEGRALE ELLIPTIQUE (PHR 2002)
  47.  
  48.         CALL BIOKAM (MCHPO1,IPT1,ENT)
  49.         IF (IERR.NE.0) RETURN
  50.         MCHPOI = ENT
  51.         GOTO 9990
  52.       ENDIF
  53. *   CHOIX DE L'OPTION DE CALCUL
  54. *   (INDUCTION PAR DEFAUT)
  55. *
  56.       INDU=.TRUE.
  57.       POTE=.FALSE.
  58.       CALL LIRMOT(CHOIX,2,IMC,0)
  59.       IF (IERR.NE.0) RETURN
  60.       IF(IMC.EQ.2) THEN
  61.          POTE=.TRUE.
  62.          INDU=.FALSE.
  63.       ENDIF
  64. *
  65. *   LECTURE D'UN MAILLAGE
  66. *
  67.       CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU)
  68.       IF (IERR.NE.0) RETURN
  69. *
  70. *   LECTURE D'UN MOT CLE
  71. *
  72.       CALL LIRMOT(MCLE,4,IMC,1)
  73.       IF (IERR.NE.0) RETURN
  74. *
  75. *         IMC = 1  ON A LU CERCLE
  76. *         IMC = 2  ON A LU ARC
  77. *         IMC = 3  ON A LU BARRE
  78. *         IMC = 4  ON A LU FIL
  79. *
  80. *
  81. *   LECTURE DE 2 ou 3 POINTS
  82. *
  83.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  84.       CALL LIROBJ('POINT ',IP2,1,IRETOU)
  85.       IP3=IP2
  86.       IF(IMC.LT.4) THEN
  87.          CALL LIROBJ('POINT ',IP3,1,IRETOU)
  88.       ELSE
  89. C        CAS DU FIL RECTILIGNE
  90.          GO TO 10
  91.       ENDIF
  92.       IF (IERR.NE.0) RETURN
  93. *
  94. *  LECTURE DE 4 OU 5 FLOTTANTS SELON LES CAS
  95. *
  96.       CALL LIRREE(F1,1,IRETOU)
  97.       CALL LIRREE(F2,1,IRETOU)
  98.       F3=0.D0
  99.       IF(IMC.LT.3) CALL LIRREE(F3,1,IRETOU)
  100. *
  101. *  LECTURE DU TYPE DE SECTION
  102. *
  103.       ISECT=0
  104.       CALL LIRMOT(CHSECT,1,ISECT,0)
  105.       IF(ISECT.EQ.1) THEN
  106.          CALL LIRREE(PENT1,1,IRETOU)
  107.          CALL LIRREE(PENT2,1,IRETOU)
  108.          IF(IERR.NE.0) RETURN
  109.       ELSE
  110.          ISECT=3
  111.       ENDIF
  112.  10   CONTINUE
  113.       CALL LIRREE(DENS,1,IRETOU)
  114.       CALL LIRREE(XIU,1,IRETOU)
  115. C   ON VIENT DE LIRE  MU0 DANS L UNITE UTILISATEUR
  116.       RUMKSA= 4.D0 * XPI * 1.D-7
  117.       RAP=(XIU / RUMKSA)
  118.       IU= 1
  119.       IF (IERR.NE.0) RETURN
  120.       ZF1(1)=F1
  121.       ZDENS(1)=DENS
  122. *
  123. *  EXEMPLE DE RECUPERATION DES COORDONNEES DU POINT P1
  124. *
  125.       IREF=(IDIM+1)*IP1
  126.       P1(1)=XCOOR(IREF-IDIM)
  127.       P1(2)=XCOOR(IREF-IDIM+1)
  128.       P1(3)=0.D0
  129.       IF (IDIM.EQ.3) P1(3)=XCOOR(IREF-IDIM+2)
  130.       IREF=(IDIM+1)*IP2
  131.       P2(1)=XCOOR(IREF-IDIM)
  132.       P2(2)=XCOOR(IREF-IDIM+1)
  133.       P2(3)=0.D0
  134.       IF (IDIM.EQ.3) P2(3)=XCOOR(IREF-IDIM+2)
  135.       IF(IMC.LT.4) THEN
  136.          IREF=(IDIM+1)*IP3
  137.          P3(1)=XCOOR(IREF-IDIM)
  138.          P3(2)=XCOOR(IREF-IDIM+1)
  139.          P3(3)=0.D0
  140.          IF (IDIM.EQ.3) P3(3)=XCOOR(IREF-IDIM+2)
  141.       ENDIF
  142. *
  143. *                -----------------------
  144. *                |   LES CALCULS ....  |
  145. *                -----------------------
  146. *
  147. *    ON COMMENCE PAR EXPLORER LE MAILLAGE
  148. *
  149.       SEGACT MELEME
  150.       IF (ITYPEL.NE.1.OR.LISOUS(/1).NE.0) CALL CHANGE(MELEME,1)
  151. *
  152. *    ON VERIFIE QUE CHAQUE POINT N'APPARAIT QU'UNE SEULE FOIS
  153. *
  154.       NNNOE=NUM(/2)
  155.       SEGINI ICPR
  156.       ICON=0
  157.       DO I=1,NNNOE
  158.          IKI=NUM(1,I)
  159.          IF(ICPR(IKI).EQ.0) THEN
  160.             ICON=ICON+1
  161.             ICPR(IKI)=ICON
  162.          ENDIF
  163.       ENDDO
  164.       SEGSUP ICPR
  165. *
  166.       IF(ICON.NE.NNNOE) THEN
  167.          CALL ERREUR(75)
  168.          RETURN
  169.       ENDIF
  170. *
  171. *    ON INITIALISE LES NOMS DES COMPOSANTES ET LES HARMONIQUES
  172. *           POUR L'INSTANT  LIMITE AU 3D ET 2D  PLAN
  173. *
  174.       NNIN=IDIM
  175.       SEGINI MTRAV
  176.       NHRM=NIFOUR
  177.       IF(NIFOUR.EQ.1) NHRM=IFOUR
  178.       INCO(1)='BX  '
  179.       IF(POTE) INCO(1)='AX  '
  180.       INCO(2)='BY  '
  181.       IF(POTE) INCO(2)='AY  '
  182.       NHAR(1)=NHRM
  183.       NHAR(2)=NHRM
  184.       IF(IDIM.EQ.3) THEN
  185.          INCO(3)='BZ  '
  186.          IF(POTE) INCO(3)='AZ  '
  187.          NHAR(3)=NHRM
  188.       ENDIF
  189. C
  190.       GOTO (100,100,200,300),IMC
  191.  100  CONTINUE
  192. C          BOBINE CIRCULAIRE
  193.  
  194. C                  ETABLISSEMENT DU REPERE LOCAL
  195. C   AXE DES Z   (P2 - P1 ) * (P3 - P1)
  196.       DO I=1,3
  197.          VX(I)=P2(I)- P1(I)
  198.          VI(I)=P3(I)- P1(I)
  199.       ENDDO
  200.       CALL PROVEC(VX,VI,VZ)
  201.       CALL PROVEC(VZ,VX,VY)
  202.       CALL ROTREP (ROTA,VX,VY,VZ)
  203.       IF(IERR.NE.0) THEN
  204.          SEGSUP MTRAV
  205.          RETURN
  206.       ENDIF
  207.       NB=1
  208.       IF(F1.EQ.0.D0)  RETURN
  209.       ALF=F2/F1
  210.       BET=F3/(2.D0*F1)
  211.       ZALF(1)=ALF
  212.       ZBET(1)=BET
  213.       IC= 3
  214.       RAYO2 = EPSIL3
  215.       PROD  = 0.D0
  216.       DO J=1,3
  217.          PP3(J)= P3(J) - P1(J)
  218.          RAYO2 = RAYO2 + VX(J) * VX(J)
  219.          PROD  = PROD  + VX(J) * VI(J)
  220.       ENDDO
  221.       TM=ACOS(PROD/RAYO2)
  222.       GOTO  350
  223.  200  CONTINUE
  224. C   CAS   DES BARRES ON VA TOURNER / OX POUR AMENER OY //  P1P2
  225.       GY=0.D0
  226.       DO  803 I=1,3
  227.          VX(I)=P2(I)- P1(I)
  228.          VI(I)=P3(I)- P1(I)
  229.          GY=GY +( VX(I)*VX(I))
  230.          C(I)= (P2(I)+ P1(I)) * 0.5D0
  231.  803  CONTINUE
  232.       GY= SQRT( GY)
  233.       CALL PROVEC(VX,VI,VZ)
  234.       CALL PROVEC(VZ,VX,VY)
  235.       CALL ROTREP (ROTA,VX,VY,VZ)
  236.       IF(IERR.NE.0) THEN
  237.          SEGSUP MTRAV
  238.          RETURN
  239.       ENDIF
  240.       DO 804 I= 1,3
  241.          VZ(I)=0.D0
  242.          DO 8041 J= 1,3
  243.             VZ(I)= VZ(I) + ROTA(I,J) * C(J)
  244.  8041    CONTINUE
  245.  804  CONTINUE
  246. *      CALL ROTATI (C,1,1,1)
  247.       GO TO 350
  248. C
  249. C     CAS DU FIL
  250.  300  CONTINUE
  251. C
  252.  350  CONTINUE
  253.       UPSI =1.D-9
  254. *
  255. *    ON BOUCLE SUR LES NOEUDS
  256. *
  257.       NT4=8000
  258.       SEGINI SANGLE
  259.       DO 30 IPT=1,NNNOE
  260.          L=NUM(1,IPT)
  261.          IREF= (IDIM+1)*L
  262. *   L EST LE NUMERO DU NOEUD  DANS LE TABLEAU DES COORDONNEES
  263.          DO  805 J=1,3
  264.             PN1(J)=XCOOR(IREF-IDIM-1+J)
  265.  805     CONTINUE
  266. *
  267.          GOTO (400,400,500,600),IMC
  268.  400     CONTINUE
  269. *   CAS DES CERCLES ET DES ARC
  270. C   MESURE DU POINT DANS LE REPERE LOCAL
  271.          DO  806 J=1,3
  272.             PN1(J)=PN1(J) - P1(J)
  273.  806     CONTINUE
  274.          DO 807 I= 1,3
  275.             VZ(I)=0.D0
  276.             VY(I)=0.D0
  277.             DO 8071 J= 1,3
  278.                VY(I)= VY(I) + ROTA(I,J) * PN1(J)
  279.                VZ(I)= VZ(I) + ROTA(I,J) * PP3(J)
  280.  8071       CONTINUE
  281.  807     CONTINUE
  282. *      CALL ROTATI (PP3,1,1,1)
  283. *      CALL ROTATI (PN1,1,1,1)
  284. C                       ON  CHANGE DE RERERENTIEL
  285. C     MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS
  286.  
  287.          ZP= VY(3)
  288.          RP= SQRT ( VY(1)**2 + VY(2)**2 )
  289.          IF(RP.EQ.0.D0) THEN
  290.             TMIN=0.D0
  291.             COST=1.D0
  292.             SINT=0.D0
  293.             RP=EPSIL2
  294.          ELSE
  295.             TMIN=ACOS(VY(1)/(RP+EPSIL3))
  296.             COST=COS(TMIN)
  297.             SINT=VY(2)/(RP+EPSIL3)
  298.          ENDIF
  299.          DO 808 J= 1,3
  300.             C(J)=0.D0
  301.  808     CONTINUE
  302. C
  303.          IF(BET.EQ.0.D0.AND.ALF.EQ.1.D0) THEN
  304. C       CAS DU FIL CIRCULAIRE
  305.             RAYON=F1
  306.             COUR=DENS
  307. C       ON CONSTRUIT LE REPERE LOCAL
  308.             sigt=1.d0
  309.             if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT)
  310.             TMIN=SIGT*(TMIN-XPI)+XPI
  311.             IF(IMC.EQ.1) TM=2.D0*XPI
  312.             TMAX=TMIN+TM
  313.             CALL TEGRAL(SANGLE)
  314.             IF(INDU) THEN
  315.                CALL BFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ)
  316.             ENDIF
  317.             IF(POTE) THEN
  318.                CALL AFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ)
  319.             ENDIF
  320.             RESP(1)=BX*COST-BY*SINT
  321.             RESP(2)=BX*SINT+BY*COST
  322.             RESP(3)=BZ
  323.          ELSE IF(BET.EQ.0.D0.AND.ISECT.EQ.3) THEN
  324. C       CAS DE LA PLAQUE CIRCULAIRE
  325.             RINT=F1
  326.             REXT=F2
  327.             DR=REXT-RINT
  328.             COUR=DENS*DR
  329. C       ON CONSTRUIT LE REPERE LOCAL
  330.             sigt=1.d0
  331.             if( abs(sint).gt.1.e-6)  SIGT=SIGN(1.D0,-SINT)
  332.             TMIN=SIGT*(TMIN-XPI)+XPI
  333.             IF(IMC.EQ.1) TM=2.D0*XPI
  334.             TMAX=TMIN+TM
  335.             CALL TEGRAL(SANGLE)
  336.             IF(INDU) THEN
  337.                CALL BPLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,
  338.      &              BX,BY,BZ)
  339.             ENDIF
  340.             IF(POTE) THEN
  341.                CALL APLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,
  342.      &              BX,BY,BZ)
  343.             ENDIF
  344.             RESP(1)=BX*COST-BY*SINT
  345.             RESP(2)=BX*SINT+BY*COST
  346.             RESP(3)=BZ
  347.          ELSE IF(ALF.EQ.1D0) THEN
  348. C       CAS DU CYLINDRE
  349.             RAYON=F1
  350.             H=F3
  351.             COUR=DENS*H
  352. C       ON CONSTRUIT LE REPERE LOCAL
  353.             sigt=1.d0
  354.             if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT)
  355.             TMIN=SIGT*(TMIN-XPI)+XPI
  356.             IF(IMC.EQ.1) TM=2.D0*XPI
  357.             TMAX=TMIN+TM
  358.             CALL TEGRAL(SANGLE)
  359.             IF(INDU) THEN
  360.                CALL BCYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,
  361.      &              BX,BY,BZ)
  362.             ENDIF
  363.             IF(POTE) THEN
  364.                CALL ACYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,
  365.      &              BX,BY,BZ)
  366.             ENDIF
  367.             RESP(1)=BX*COST-BY*SINT
  368.             RESP(2)=BX*SINT+BY*COST
  369.             RESP(3)=BZ
  370.          ELSE
  371.             IF (IMC.EQ.1.AND.ISECT.EQ.3) THEN
  372.                IF(INDU) THEN
  373.                   CALL CHAMBO(IU,IC,NB,C(3),ZF1,ZALF,ZBET,ZDENS,ZP,RP,BZ
  374.      $                 ,BR,BP)
  375.                   IF (RP.LE.UPSI) THEN
  376.                      RESP(1)=0.D0
  377.                      RESP(2)=0.D0
  378.                   ELSE
  379.                      RESP(1)=BR * VY(1) / RP
  380.                      RESP(2)=BR * VY(2) / RP
  381.                   END IF
  382.                   RESP(3)=BZ
  383.                   GO TO 700
  384.                ENDIF
  385.                IF(POTE) THEN
  386.                   RINT=F1
  387.                   REXT=F2
  388.                   HCEN=F3
  389.                   WIDTH=F2-F1
  390.                   AIRE=F3*WIDTH
  391.                   COUR=DENS*AIRE
  392.                   DELTA=0.D0
  393. C         ON CONSTRUIT LE REPERE LOCAL
  394.                   SIGT=1.D0
  395.                   if(ABS(SINT).gt.1.e-8)SIGT=SIGN(1.D0,-SINT)
  396.                   TMIN=SIGT*(TMIN-XPI)+XPI
  397.                   TM=2.D0*XPI
  398.                   TMAX=TMIN+TM
  399.                   PENT1=0.D0
  400.                   PENT2=0.D0
  401.                   CALL TEGRAL(SANGLE)
  402.                   CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT,
  403.      *                 PENT1,PENT2,TMIN,TMAX,
  404.      *                 RP,ZP,BX,BY,BZ)
  405.                   RESP(1)=BX*COST-BY*SINT
  406.                   RESP(2)=BX*SINT+BY*COST
  407.                   RESP(3)=BZ
  408.                   GO TO 700
  409.                ENDIF
  410.             ENDIF
  411.             IF(ISECT.EQ.3) THEN
  412. C       CAS DES SECTIONS RECTANGULAIRES
  413. C       ON CALCULE POUR UN ARC  ---> ANGLES
  414.                A1=0.D0
  415.                IF (ABS(VZ(1)).LT.UPSI) THEN
  416.                   A2=XPI/2.D0
  417.                ELSE
  418.                   AA= VZ(2)/VZ(1)
  419.                   A2=ATAN ( AA )
  420.                   IF ( A2.LT.0.D0) A2 = XPI + A2
  421.                END IF
  422.                IF(INDU) THEN
  423.                   CALL CHAMAR(IU,IC,C,A1,A2,F1,ALF,BET,DENS,
  424.      *                 VY(1),VY(2),VY(3),BX,BY,BZ)
  425.                   RESP(1)=BX
  426.                   RESP(2)=BY
  427.                   RESP(3)=BZ
  428.                ENDIF
  429.                IF(POTE) THEN
  430.                   RINT=F1
  431.                   REXT=F2
  432.                   HCEN=F3
  433.                   WIDTH=F2-F1
  434.                   AIRE=F3*WIDTH
  435.                   COUR=DENS*AIRE
  436.                   PENT1=0.D0
  437.                   PENT2=0.D0
  438. C         ON CONSTRUIT LE REPERE LOCAL
  439.                   sigt=1.d0
  440.                   if( abs(sint).gt.1.e-6)SIGT=SIGN(1.d0,-SINT)
  441.                   TMIN=SIGT*(TMIN-XPI)+XPI
  442.                   TMAX=TMIN+TM
  443.                   CALL TEGRAL(SANGLE)
  444.                   CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT,
  445.      *                 PENT1,PENT2,TMIN,TMAX,
  446.      *                 RP,ZP,BX,BY,BZ)
  447.                   RESP(1)=BX*COST-BY*SINT
  448.                   RESP(2)=BX*SINT+BY*COST
  449.                   RESP(3)=BZ
  450.                ENDIF
  451.             ELSEIF(ISECT.EQ.1) THEN
  452. C     CAS DES SECTIONS TRAPEZOIDALES
  453.                RINT=F1
  454.                REXT=F2
  455.                HCEN=F3
  456.                WIDTH=F2-F1
  457.                AIRE=F3*WIDTH
  458.                COUR=DENS*AIRE
  459.                DELTA=0.5D0*WIDTH*(PENT2-PENT1)
  460.                HINT=HCEN-DELTA
  461.                HEXT=HCEN+DELTA
  462. C       ON CONSTRUIT LE REPERE LOCAL
  463.                SIGT=1.d0
  464.                if(abs(sint).gt.1.e-06)SIGT=SIGN(1.D0,-SINT)
  465.                TMIN=SIGT*(TMIN-XPI)+XPI
  466.                IF(IMC.EQ.1) TM=2.D0*XPI
  467.                TMAX=TMIN+TM
  468.                CALL TEGRAL(SANGLE)
  469.                IF(HCEN.EQ.0.D0) THEN
  470. C         CAS DU TRONC DE CONE
  471.                   TETA=ATAN(PENT1)
  472.                   COUR=DENS*WIDTH/COS(TETA)
  473.                   IF(INDU) THEN
  474.                      CALL BTRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX,
  475.      &                    RP,ZP,BX,BY,BZ)
  476.                   ENDIF
  477.                   IF(POTE) THEN
  478.                      CALL ATRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX,
  479.      &                    RP,ZP,BX,BY,BZ)
  480.                   ENDIF
  481.                ELSE
  482.                   IF(INDU) THEN
  483.                      CALL ARCTRA(SANGLE,COUR,HINT,HEXT,RINT,REXT,
  484.      *                    PENT1,PENT2,TMIN,TMAX,
  485.      *                    RP,ZP,BX,BY,BZ)
  486.                   ENDIF
  487.                   IF(POTE) THEN
  488.                      CALL AARCTR(SANGLE,COUR,HINT,HEXT,RINT,REXT,
  489.      *                    PENT1,PENT2,TMIN,TMAX,
  490.      *                    RP,ZP,BX,BY,BZ)
  491.                   ENDIF
  492.                ENDIF
  493.                RESP(1)=BX*COST-BY*SINT
  494.                RESP(2)=BX*SINT+BY*COST
  495.                RESP(3)=BZ
  496.             ENDIF
  497.          END IF
  498.          IF(IERR.NE.0) THEN
  499.             SEGSUP MTRAV
  500.             RETURN
  501.          ENDIF
  502. C
  503.          GOTO   700
  504. *
  505.  500     CONTINUE
  506. C      CAS DES BARRES         ON  CHANGE DE RERERENTIEL
  507.          DO 809 I= 1,3
  508.             VY(I)=0.D0
  509.             DO 8091 J= 1,3
  510.                VY(I)= VY(I) + ROTA(I,J) * PN1(J)
  511.  8091       CONTINUE
  512.  809     CONTINUE
  513. *      CALL ROTATI (PN1,1,1,1)
  514. C     MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS
  515.          IF(ISECT.EQ.3) THEN
  516. C     CAS DES SECTIONS RETANGULAIRES
  517.             IF(F2.EQ.0.D0) THEN
  518. C         CAS DE LA PLAQUE RECTANGULAIRE
  519.                DO 813 I=1,3
  520.                   VY(I)=0.D0
  521.                   DO 8131 J=1,3
  522.                      VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
  523.  8131             CONTINUE
  524.  813           CONTINUE
  525.                DY=F1
  526.                REXT=0.5D0*DY
  527.                RINT=-REXT
  528.                COUR=DENS*DY
  529.                YMIN=-VY(1)
  530.                YMAX=YMIN+GY
  531.                XP=-VY(2)
  532.                ZP=VY(3)
  533.                IF(INDU) THEN
  534.                   CALL BPLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP
  535.      $                 )
  536.                ENDIF
  537.                IF(POTE) THEN
  538.                   CALL APLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP
  539.      $                 )
  540.                ENDIF
  541.                RESP(1)=BYP
  542.                RESP(2)=-BXP
  543.                RESP(3)= BZP
  544.             ELSE
  545.                IF(INDU) THEN
  546.                   CALL CHABAR(IU,VZ,GY,F1,F2,DENS,VY(1),VY(2),VY(3),
  547.      *                 BXP,BYP,BZP)
  548.                   RESP(1)=BXP
  549.                   RESP(2)=BYP
  550.                   RESP(3)=BZP
  551.                ENDIF
  552.                IF(POTE) THEN
  553.                   DO 814 I=1,3
  554.                      VY(I)=0.D0
  555.                      DO 8141 J=1,3
  556.                         VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
  557.  8141                CONTINUE
  558.  814              CONTINUE
  559.                   DY=F1
  560.                   DZ=F2
  561.                   REXT=0.5D0*DY
  562.                   RINT=-REXT
  563.                   AIRE=DY*DZ
  564.                   COUR=DENS*AIRE
  565.                   PENT1=0.D0
  566.                   PENT2=0.D0
  567.                   YMIN=-VY(1)
  568.                   YMAX=YMIN+GY
  569.                   XP=-VY(2)
  570.                   ZP=VY(3)
  571.                   CALL ASEGTR(COUR,DZ,DZ,RINT,REXT,PENT1,PENT2,
  572.      *                 YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
  573.                   RESP(1)=BYP
  574.                   RESP(2)=-BXP
  575.                   RESP(3)= BZP
  576.                ENDIF
  577.             ENDIF
  578.          ELSE IF(ISECT.EQ.1) THEN
  579. C     CAS DES SECTIONS TRAPEZOIDALES
  580.             DO 812 I=1,3
  581.                VY(I)=0.D0
  582.                DO 8121 J=1,3
  583.                   VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
  584.  8121          CONTINUE
  585.  812        CONTINUE
  586.             DY=F1
  587.             DZ=F2
  588.             REXT=0.5D0*DY
  589.             RINT=-REXT
  590.             AIRE=DY*DZ
  591.             COUR=DENS*AIRE
  592.             DELTA=0.5D0*DY*(PENT2-PENT1)
  593.             HINT=DZ-DELTA
  594.             HEXT=DZ+DELTA
  595.             YMIN=-VY(1)
  596.             YMAX=YMIN+GY
  597.             XP=-VY(2)
  598.             ZP=VY(3)
  599.             IF(INDU) THEN
  600.                CALL SEGTRA(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2,
  601.      *              YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
  602.             ENDIF
  603.             IF(POTE) THEN
  604.                CALL ASEGTR(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2,
  605.      *              YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
  606.             ENDIF
  607.             RESP(1)=BYP
  608.             RESP(2)=-BXP
  609.             RESP(3)= BZP
  610.          ENDIF
  611.          GO TO 700
  612. C
  613. C     CAS DU FIL
  614.  600     CONTINUE
  615.          GY=0.D0
  616.          DO  601 I=1,3
  617.             VY(I)=P2(I)-P1(I)
  618.             VI(I)=PN1(I)-P1(I)
  619.             GY=GY+(VY(I)*VY(I))
  620.  601     CONTINUE
  621.          GY=SQRT(GY)
  622.          CALL PROVEC(VI,VY,VZ)
  623.          CALL PROVEC(VY,VZ,VX)
  624.          CALL ROTREP (ROTA,VX,VY,VZ)
  625.          IF(IERR.NE.0) THEN
  626.             SEGSUP MTRAV
  627.             RETURN
  628.          ENDIF
  629.          XP=0.D0
  630.          YMIN=0.D0
  631.          DO 602 J= 1,3
  632.             XP=XP+ROTA(1,J)*VI(J)
  633.             YMIN=YMIN-ROTA(2,J)*VI(J)
  634.  602     CONTINUE
  635.          YMAX=YMIN+GY
  636.          COUR=DENS
  637.          IF(INDU) THEN
  638.             CALL BFILDR(COUR,YMIN,YMAX,XP,BZP)
  639.             RESP(1)=0.D0
  640.             RESP(2)=0.D0
  641.             RESP(3)=BZP
  642.          ENDIF
  643.          IF(POTE) THEN
  644.             CALL AFILDR(COUR,YMIN,YMAX,XP,AYP)
  645.             RESP(1)=0.D0
  646.             RESP(2)=AYP
  647.             RESP(3)=0.D0
  648.          ENDIF
  649. CC
  650.  700     CONTINUE
  651. *
  652. C     RETOUR  DU CHAMP DANS LE REPERE  GENERAL
  653.  
  654.          DO 810 I= 1,3
  655.             VY(I)=0.D0
  656.             DO 8101 J= 1,3
  657.                VY(I)= VY(I) + ROTA(J,I) * RESP(J) * RAP
  658.  8101       CONTINUE
  659.  810     CONTINUE
  660. *      CALL RET_REF   (RESP)
  661. *   REMPLISSAGE DU SEGMENT DE TRAVAIL MTRAV
  662. *
  663.          IGEO(IPT)=L
  664.          DO 811 J= 1,3
  665.             BB(J,IPT) = VY(J)
  666.             IBIN(J,IPT) = 1
  667.  811     CONTINUE
  668. C          FIN BOUCLE SUR LES POINTS
  669. *
  670.  30   CONTINUE
  671.       SEGSUP SANGLE
  672. *
  673.       CALL CRECHP(MTRAV,MCHPOI)
  674.       SEGSUP MTRAV
  675. *
  676. *       ATTRIBUTION D'UNE NATURE DISCRETE AU CHAMP QUI SORT
  677. *
  678.  9990 SEGACT MCHPOI*MOD
  679.       JATTRI(1) = 2
  680. *
  681.       CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  682.       CALL ECROBJ('CHPOINT ',MCHPOI)
  683.  
  684. c      return
  685.       END
  686.  
  687.  
  688.  

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