Télécharger la plaquette Cast3M


fr en

Télécharger permab.eso

Retour à la liste

Numérotation des lignes : 

permab

  1. C PERMAB    SOURCE    OF166741  25/02/21    21:18:07     12166          
  2.       SUBROUTINE PERMAB (MODORI,MCHORI,IPRIGI,IRET)
  3. *______________________________________________________________________
  4. *
  5. *      OPERATEUR PERMEABILITE (MILIEUX POREUX) APPELE PAR PERMEA
  6. *
  7. *  Entrees :
  8. *  ---------
  9. *
  10. *     MODORI   Pointeur sur un MMODEL
  11. *     MCHORI   Pointeur sur un MCHAML de materiau
  12. *
  13. *  Sorties :
  14. *  ---------
  15. *
  16. *     IPRIGI   Pointeur sur un objet RIGIDITE de permeabilite
  17. *     IRET     =1 ou 0 suivant succes ou non
  18. *
  19. *  Passage aux nouveaux CHAMELEMs par jm CAMPENON le 07/91
  20. *______________________________________________________________________
  21. *
  22.       IMPLICIT INTEGER(I-N)
  23.       IMPLICIT REAL*8(A-H,O-Z)
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC CCGEOME
  28. -INC CCHAMP
  29.  
  30. -INC SMRIGID
  31. -INC SMCHAML
  32. -INC SMELEME
  33. -INC SMCOORD
  34. -INC SMINTE
  35. -INC SMMODEL
  36.  
  37. -INC TMPTVAL
  38.  
  39.       SEGMENT WRK1
  40.          REAL*8 REL(LRE,LRE),XE(3,NBBB)
  41.       ENDSEGMENT
  42. *
  43.       SEGMENT WRK2
  44.          REAL*8 SHPWRK(6,NBNO),BGENE(NSTB,LRE)
  45.       ENDSEGMENT
  46. *
  47.       SEGMENT WRK3
  48.          REAL*8 XGENE(NSTN,LRN)
  49.       ENDSEGMENT
  50. *
  51.       SEGMENT WRK4
  52.          REAL*8 XLOC(3,3),XGLOB(3,3),TXR(IDIM,IDIM)
  53.          REAL*8 VALMAT(NMATT)
  54.          REAL*8 PMAT(NSTB,NSTB),PMAT1(IDIM,IDIM),PMAT2(IDIM,IDIM)
  55.       ENDSEGMENT
  56. *
  57.       SEGMENT WRK5
  58.        REAL*8 BPSS(3,3),XEL(3,NBBB)
  59.        REAL*8 XNTH(LPP,LPP),XNTB(LPP,LPP),XNTT(LPP)
  60.       ENDSEGMENT
  61. *
  62.       SEGMENT WRK6
  63.        REAL*8 PKK(NSTPK,NSTPK)
  64.       ENDSEGMENT
  65. *
  66.       SEGMENT INFO
  67.        INTEGER INFELL(JG)
  68.       ENDSEGMENT
  69. *
  70.       SEGMENT NOTYPE
  71.         CHARACTER*16 TYPE(NBTYPE)
  72.       ENDSEGMENT
  73.  
  74.       INTEGER OOOVAL
  75.  
  76.       CHARACTER*8 CMATE
  77.       CHARACTER*(NCONCH) CONM
  78.  
  79. * INTTYP correspond au type de points d'integration utilise
  80.       PARAMETER ( INTTYP=3 )
  81.  
  82.       PARAMETER (NINF=3)
  83.       INTEGER INFOS(NINF)
  84.       LOGICAL lsupfo,lsupdp
  85. *
  86.       NHRM=NIFOUR
  87.       IRET = 0
  88.       IPRIGI = 0
  89. *
  90. *     Reduction du modele a la formulation poreuse
  91. *
  92.       MMODE1 = MODORI
  93.       SEGINI,MMODEL=MMODE1
  94.       NSOUS = KMODEL(/1)
  95.       N1 = 0
  96.       DO isous = 1, NSOUS
  97.         IMODEL = KMODEL(isous)
  98.         SEGACT,IMODEL
  99.         IF (FORMOD(1).EQ.'POREUX') THEN
  100.           N1 = N1 + 1
  101.           KMODEL(N1) = IMODEL
  102.         ELSE
  103.           SEGDES,IMODEL
  104.         ENDIF
  105.       ENDDO
  106.       IF (N1.NE.NSOUS) SEGADJ,MMODEL
  107.       IPMODL = MMODEL
  108.       NSOUS = N1
  109.       IF (NSOUS.LE.0) THEN
  110.         MOTERR(1:8) = 'MMODEL  '
  111.         INTERR(1) = MODORI
  112.         CALL ERREUR(356)
  113.         GOTO 9991
  114.       ENDIF
  115.  
  116. *
  117. *     Reduction du champ au modele precedemment reduit
  118. *
  119.       MCHELM = MCHORI
  120.       SEGACT,MCHELM
  121.       IF (TITCHE(1:8).NE.'CARACTER') THEN
  122.         SEGDES,MCHELM
  123.         MOTERR(1:16) = 'CARACTERISTIQUES'
  124.         CALL ERREUR(291)
  125.         GOTO 9991
  126.       ENDIF
  127.       CALL REDUAF(MCHORI,IPMODL,IPCHE1,0,IRET,KERRE)
  128.       SEGDES,MCHELM
  129.       IF (IRET.NE.1) THEN
  130.         CALL ERREUR(KERRE)
  131.         GOTO 9991
  132.       ENDIF
  133. *
  134. *     Verification du lieu support du MCHAML de materiau
  135. *
  136.       ISUP=0
  137.       CALL  QUESUP(IPMODL,IPCHE1,INTTYP,0,ISUP,IRETMA)
  138.       IF (ISUP.GT.1) GO TO 9991
  139. *
  140. *     Activation du MMODEL
  141. *
  142.       MMODEL=IPMODL
  143.       SEGACT MMODEL
  144.       NSOUS=KMODEL(/1)
  145. *
  146. *  ON FABRIQUE LES MATRICES UNIQUEMENT POUR LES ZONES
  147. *  DE MILIEU POREUX.
  148. *
  149.        NRIGEL=NSOUS
  150. *
  151. *     Initialisation du chapeau de l'objet RIGIDITE
  152. *
  153.       SEGINI MRIGID
  154.       ICHOLE=0
  155.       IMGEO1=0
  156.       IMGEO2=0
  157.       IFORIG=IFOUR
  158.       MTYMAT='PERMEABI'
  159. *
  160. *      BOUCLE SUR LES SOUS ZONES DU MODELE
  161. *
  162.       ISORI=0
  163.       DO 500 ISOUS=1,NSOUS
  164. *
  165. *        On recupere l'information generale
  166. *
  167.          IMODEL=KMODEL(ISOUS)
  168.          SEGACT IMODEL
  169.  
  170.          IF(FORMOD(1).NE.'POREUX') THEN
  171.             CALL ERREUR(19)
  172.             GO TO 9999
  173.          ENDIF
  174. *
  175. *        Traitement du modele
  176. *
  177.          MELE=NEFMOD
  178.          IPMAIL=IMAMOD
  179.          CONM  =CONMOD
  180.  
  181. *
  182.          CALL IDENT(IPMAIL,CONM,IPCHE1,0,INFOS,IRTD)
  183.          IF (IRTD.EQ.0) GOTO 9999
  184. *
  185. *        Nature du materiau
  186. *
  187.          CMATE = CMATEE
  188.          MATE  = IMATEE
  189.          INAT  = INATUU
  190. *
  191. *        Information sur l'element fini
  192. *
  193.         IF (INFMOD(/1).LT.2+INTTYP) THEN
  194.          CALL ELQUOI(MELE,0,INTTYP,IPINF,IMODEL)
  195.          IF (IERR.NE.0) GO TO 9999
  196.          INFO=IPINF
  197.          MFR   =INFELL(13)
  198.          IELE  =INFELL(14)
  199.          IPORE=INFELL(8)
  200.          MINTE=INFELL(11)
  201.          segsup info
  202.          else
  203.          MFR   =INFELE(13)
  204.          IELE  =INFELE(14)
  205.          IPORE=INFELE(8)
  206.          MINTE=infmod(5)
  207.          endif
  208.          IPMINT=MINTE
  209. *
  210. *  Si necessaire PARTITIONNEMENT du segment XMATRI
  211. *
  212.         IPT1=IPMAIL
  213.         SEGACT,IPT1
  214.         NBNN1 =IPT1.NUM(/1)
  215.         NBELE1=IPT1.NUM(/2)
  216. *
  217.         LASYM=0
  218.         IF(MFR.EQ.33) THEN
  219.             IDECAP=1
  220.         ELSE IF(MFR.EQ.57) THEN
  221.             IDECAP=2
  222.             LASYM=2
  223.         ELSE IF(MFR.EQ.59) THEN
  224.             IDECAP=3
  225.             LASYM=2
  226.         ENDIF
  227.         LR1=NBNNE(IELE)
  228.         LRE=LR1*IDECAP
  229. *
  230.         LTRK=OOOVAL(1,4)
  231.         IF (LTRK.EQ.0) LTRK=OOOVAL(1,1)
  232.       LTRK=MAX(LTRK,2**24)
  233. * Ajout a la taille en mots de la matrice des infos du segment
  234.         LSEG=LRE*LRE*NBELE1 + 16
  235.         NBLPRT=(LSEG-1)/LTRK+1
  236.         NBLMAX=(NBELE1-1)/NBLPRT+1
  237.         NBLPRT=(NBELE1-1)/NBLMAX+1
  238. *        write(ioimp,*) ' PERMAB nblprt nblmax ',NBLPRT,NBLMAX,NBELE1
  239.         MELEME=IPT1
  240.  
  241. *   BOUCLE DE PARTITIONNEMENT DU SEGMENT XMATRI
  242.  
  243.       DO 5000 IPRT = 1,NBLPRT
  244.         ISORI= ISORI+1
  245.         IF (ISORI.GT.IRIGEL(/2)) THEN
  246.           NRIGEL=ISORI
  247.           SEGADJ,MRIGID
  248.         ENDIF
  249.         IF (NBLPRT.GT.1) THEN
  250.           JPRT=(IPRT-1)*NBLMAX
  251.           SEGACT,IPT1
  252.           NBSOUS=0
  253.           NBREF=0
  254.           NBNN=NBNN1
  255.           NBELEM=MIN(NBLMAX,NBELE1-JPRT)
  256. *          write(ioimp,*) ' creation segment ',nbnn,nbelem
  257.           SEGINI,MELEME
  258.           ITYPEL=IPT1.ITYPEL
  259.           DO I=1,NBELEM
  260.             IB=I+JPRT
  261.             DO J=1,NBNN
  262.               NUM(J,I)=IPT1.NUM(J,IB)
  263.             ENDDO
  264.             ICOLOR(I)=IPT1.ICOLOR(I)
  265.           ENDDO
  266.         ENDIF
  267.         IPMAIL=MELEME
  268. *  Fin du traitement particulier en cas de PARTITIONNEMENT du XMATRI
  269. *  Quelques initialisations suite au partionnement
  270.         IPDES  = 0
  271. *
  272.         NMATR  = 0
  273.         NMATF  = 0
  274.         IVAMAT = 0
  275.         NCARA  = 0
  276.         NCARF  = 0
  277.         IVACAR = 0
  278.  
  279. *  Activation du MELEME support des rigidites
  280.         MELEME=IPMAIL
  281.         SEGACT,MELEME
  282.         NBNN  =NUM(/1)
  283.         NBELEM=NUM(/2)
  284.  
  285.          NLIGRP = LRE
  286.          NLIGRD = LRE
  287.          SEGINI DESCR
  288.          if(lnomid(1).ne.0) then
  289.            nomid=lnomid(1)
  290.            segact nomid
  291.            modepl=nomid
  292.            ndepl=lesobl(/2)
  293.            ndum=lesfac(/2)
  294.            lsupdp=.false.
  295.          else
  296.            lsupdp=.true.
  297.            CALL IDPRIM(IMODEL,MFR,MODEPL,NDEPL,NDUM)
  298.          endif
  299.          if(lnomid(2).ne.0) then
  300.           nomid=lnomid(2)
  301.           segact nomid
  302.           moforc=nomid
  303.           nforc=lesobl(/2)
  304.           lsupfo=.false.
  305.          else
  306.           lsupfo=.true.
  307.           CALL IDDUAL(IMODEL,MFR,MOFORC,NFORC,NDUM)
  308.          endif
  309. *
  310.          ID=1
  311.          NOMID=MODEPL
  312.          SEGACT NOMID
  313.          NCP=LESOBL(/2)
  314.          NOMID=MOFORC
  315.          SEGACT NOMID
  316. *
  317.          IF (MFR.EQ.33) THEN
  318.  
  319.            DO 4005 IB=1,NBSOM(IELE)
  320.               NOMID=MODEPL
  321.               LISINC(ID)=LESOBL(NCP)
  322.               NOMID=MOFORC
  323.               LISDUA(ID)=LESOBL(NCP)
  324.               NOELEP(ID)=IBSOM(NSPOS(IELE)+IB-1)
  325.               NOELED(ID)=IBSOM(NSPOS(IELE)+IB-1)
  326.               ID=ID+1
  327.  4005      CONTINUE
  328. *
  329.            IF (MELE.GE.108.AND.MELE.LE.110) THEN
  330. *
  331.                LR1=(3*LRE-IPORE)/2
  332.                DO 4008 INOEUD=LR1+1,LRE
  333.                  NOMID=MODEPL
  334.                  LISINC(ID)=LESOBL(NCP)
  335.                  NOMID=MOFORC
  336.                  LISDUA(ID)=LESOBL(NCP)
  337.                  NOELEP(ID)=INOEUD
  338.                  NOELED(ID)=INOEUD
  339.                  ID=ID+1
  340. 4008           CONTINUE
  341.            END IF
  342. *
  343.            DO 4006 IB=1,LR1
  344.               DO 4007 INSOM=1,NBSOM(IELE)
  345.                 IF(IB.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4006
  346. 4007         CONTINUE
  347.              NOMID=MODEPL
  348.              LISINC(ID)=LESOBL(NCP)
  349.              NOMID=MOFORC
  350.              LISDUA(ID)=LESOBL(NCP)
  351.              NOELEP(ID)=IB
  352.              NOELED(ID)=IB
  353.              ID=ID+1
  354. 4006       CONTINUE
  355. *
  356.          ELSE IF (MFR.EQ.57.OR.MFR.EQ.59) THEN
  357.  
  358.            DO 4205 IPR=1,IDECAP
  359.               NCPDEC=NCP-IDECAP+IPR
  360. *
  361.               DO 4105 IB=1,NBSOM(IELE)
  362.                  NOMID=MODEPL
  363.                  LISINC(ID)=LESOBL(NCPDEC)
  364.                  NOMID=MOFORC
  365.                  LISDUA(ID)=LESOBL(NCPDEC)
  366.                  NOELEP(ID)=IBSOM(NSPOS(IELE)+IB-1)
  367.                  NOELED(ID)=IBSOM(NSPOS(IELE)+IB-1)
  368.                  ID=ID+1
  369.  4105         CONTINUE
  370. *
  371.            IF (MELE.GE.185.AND.MELE.LE.190) THEN
  372. *
  373.                LR1=(3*NBNNE(IELE)-IPORE)/2
  374.                DO 4108 INOEUD=LR1+1,NBNNE(IELE)
  375.                  NOMID=MODEPL
  376.                  LISINC(ID)=LESOBL(NCPDEC)
  377.                  NOMID=MOFORC
  378.                  LISDUA(ID)=LESOBL(NCPDEC)
  379.                  NOELEP(ID)=INOEUD
  380.                  NOELED(ID)=INOEUD
  381.                  ID=ID+1
  382. 4108           CONTINUE
  383.            END IF
  384. * FIN NEW
  385. *
  386.              DO 4106 IB=1,LR1
  387.                 DO 4107 INSOM=1,NBSOM(IELE)
  388.                   IF(IB.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4106
  389. 4107            CONTINUE
  390.                 NOMID=MODEPL
  391.                 LISINC(ID)=LESOBL(NCPDEC)
  392.                 NOMID=MOFORC
  393.                 LISDUA(ID)=LESOBL(NCPDEC)
  394.                 NOELEP(ID)=IB
  395.                 NOELED(ID)=IB
  396.                 ID=ID+1
  397. 4106         CONTINUE
  398.  
  399. 4205       CONTINUE
  400. *
  401.          ENDIF
  402.  
  403.          NOMID =MOFORC
  404.          if(lsupfo)SEGSUP NOMID
  405.          NOMID =MODEPL
  406.          if(lsupdp)SEGSUP NOMID
  407.          IPDES=DESCR
  408.          SEGDES DESCR
  409. *
  410. *        Initialisation de MINTE
  411. *
  412.          SEGACT MINTE
  413.          NBPGAU=POIGAU(/1)
  414. *
  415.  
  416. CCCCCC    LVAL=(LRE*(LRE+1))/2
  417.          NELRIG=NBELEM
  418.          SEGINI xMATRI
  419. *
  420. *        Verification de la presence des composantes pour le materiau
  421. *
  422.          NBROBL=0
  423.          NBRFAC=0
  424. *                                         cas isotrope
  425.             IF (MATE.EQ.1) THEN
  426. *
  427.               IF (MELE.GE.79.AND.MELE.LE.83) THEN
  428.                 NBROBL=2
  429.                 SEGINI NOMID
  430.                 LESOBL(1)='PERM'
  431.                 LESOBL(2)='VISC'
  432.               ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN
  433.                 NBROBL=4
  434.                 SEGINI NOMID
  435.                 LESOBL(1)='PERT'
  436.                 LESOBL(2)='PERH'
  437.                 LESOBL(3)='PERB'
  438.                 LESOBL(4)='VISC'
  439.               ELSE IF (MELE.GE.173.AND.MELE.LE.177) THEN
  440.                 NBROBL=4
  441.                 SEGINI NOMID
  442.                 LESOBL(1)='PK11'
  443.                 LESOBL(2)='PK12'
  444.                 LESOBL(3)='PK21'
  445.                 LESOBL(4)='PK22'
  446.               ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN
  447.                 NBROBL=9
  448.                 SEGINI NOMID
  449.                 LESOBL(1)='PK11'
  450.                 LESOBL(2)='PK12'
  451.                 LESOBL(3)='PK13'
  452.                 LESOBL(4)='PK21'
  453.                 LESOBL(5)='PK22'
  454.                 LESOBL(6)='PK23'
  455.                 LESOBL(7)='PK31'
  456.                 LESOBL(8)='PK32'
  457.                 LESOBL(9)='PK33'
  458.               ELSE IF (MELE.GE.185.AND.MELE.LE.187) THEN
  459.                 NBROBL=12
  460.                 SEGINI NOMID
  461.                 LESOBL(1)='PT11'
  462.                 LESOBL(2)='PH11'
  463.                 LESOBL(3)='PB11'
  464.                 LESOBL(4)='PT12'
  465.                 LESOBL(5)='PH12'
  466.                 LESOBL(6)='PB12'
  467.                 LESOBL(7)='PT21'
  468.                 LESOBL(8)='PH21'
  469.                 LESOBL(9)='PB21'
  470.                 LESOBL(10)='PT22'
  471.                 LESOBL(11)='PH22'
  472.                 LESOBL(12)='PB22'
  473.               ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN
  474.                 NBROBL=27
  475.                 SEGINI NOMID
  476.                 LESOBL(1)='PT11'
  477.                 LESOBL(2)='PH11'
  478.                 LESOBL(3)='PB11'
  479.                 LESOBL(4)='PT12'
  480.                 LESOBL(5)='PH12'
  481.                 LESOBL(6)='PB12'
  482.                 LESOBL(7)='PT13'
  483.                 LESOBL(8)='PH13'
  484.                 LESOBL(9)='PB13'
  485.                 LESOBL(10)='PT21'
  486.                 LESOBL(11)='PH21'
  487.                 LESOBL(12)='PB21'
  488.                 LESOBL(13)='PT22'
  489.                 LESOBL(14)='PH22'
  490.                 LESOBL(15)='PB22'
  491.                 LESOBL(16)='PT23'
  492.                 LESOBL(17)='PH23'
  493.                 LESOBL(18)='PB23'
  494.                 LESOBL(19)='PT31'
  495.                 LESOBL(20)='PH31'
  496.                 LESOBL(21)='PB31'
  497.                 LESOBL(22)='PT32'
  498.                 LESOBL(23)='PH32'
  499.                 LESOBL(24)='PB32'
  500.                 LESOBL(25)='PT33'
  501.                 LESOBL(26)='PH33'
  502.                 LESOBL(27)='PB33'
  503.               ENDIF
  504. *                                          cas orthotrope
  505.             ELSE IF (MATE.EQ.2) THEN
  506.               IF (IDIM.EQ.3) THEN
  507.                 NBROBL=10
  508.                 SEGINI NOMID
  509.                 LESOBL(1)='PER1'
  510.                 LESOBL(2)='PER2'
  511.                 LESOBL(3)='PER3'
  512.                 LESOBL(4)='VISC'
  513.                 LESOBL(5)='V1X '
  514.                 LESOBL(6)='V1Y '
  515.                 LESOBL(7)='V1Z '
  516.                 LESOBL(8)='V2X '
  517.                 LESOBL(9)='V2Y '
  518.                 LESOBL(10)='V2Z '
  519.               ELSE IF(IDIM.EQ.2) THEN
  520.                 IF (IFOUR.LE.0) THEN
  521.                   NBROBL=5
  522.                   SEGINI NOMID
  523.                   LESOBL(1)='PER1'
  524.                   LESOBL(2)='PER2'
  525.                   LESOBL(3)='VISC'
  526.                   LESOBL(4)='V1X '
  527.                   LESOBL(5)='V1Y '
  528.                 ELSE IF (IFOUR.EQ.1) THEN
  529.                   NBROBL=6
  530.                   SEGINI NOMID
  531.                   LESOBL(1)='PER1'
  532.                   LESOBL(2)='PER2'
  533.                   LESOBL(3)='PER3'
  534.                   LESOBL(4)='VISC'
  535.                   LESOBL(5)='V1X '
  536.                   LESOBL(6)='V1Y '
  537.                 ENDIF
  538.               ENDIF
  539. *                                          cas anisotrope
  540.             ELSE IF (MATE.EQ.3)THEN
  541.               IF(IDIM.EQ.3)THEN
  542.                 NBROBL=13
  543.                 SEGINI NOMID
  544.                 LESOBL(1)='PER1'
  545.                 LESOBL(2)='PER2'
  546.                 LESOBL(3)='PER3'
  547.                 LESOBL(4)='PE12'
  548.                 LESOBL(5)='PE13'
  549.                 LESOBL(6)='PE23'
  550.                 LESOBL(7)='VISC'
  551.                 LESOBL(8)='V1X '
  552.                 LESOBL(9)='V1Y '
  553.                 LESOBL(10)='V1Z '
  554.                 LESOBL(11)='V2X '
  555.                 LESOBL(12)='V2Y '
  556.                 LESOBL(13)='V2Z '
  557.               ELSE IF (IDIM.EQ.2) THEN
  558.                 IF (IFOUR.LE.0) THEN
  559.                   NBROBL=6
  560.                   SEGINI NOMID
  561.                   LESOBL(1)='PER1'
  562.                   LESOBL(2)='PER2'
  563.                   LESOBL(3)='PE12'
  564.                   LESOBL(4)='VISC'
  565.                   LESOBL(5)='V1X '
  566.                   LESOBL(6)='V1Y '
  567.                 ELSE IF (IFOUR.EQ.1) THEN
  568.                   NBROBL=7
  569.                   SEGINI NOMID
  570.                   LESOBL(1)='PER1'
  571.                   LESOBL(2)='PER2'
  572.                   LESOBL(3)='PE12'
  573.                   LESOBL(4)='PER3'
  574.                   LESOBL(5)='VISC'
  575.                   LESOBL(6)='V1X '
  576.                   LESOBL(7)='V1Y '
  577.                 ENDIF
  578.               ENDIF
  579. *                                         cas unidirectionnel
  580.             ELSE IF (MATE.EQ.4) THEN
  581.              IF (IDIM.EQ.3) THEN
  582.               NBROBL=8
  583.               SEGINI NOMID
  584.               LESOBL(1)='PERM'
  585.               LESOBL(2)='VISC'
  586.               LESOBL(3)='V1X '
  587.               LESOBL(4)='V1Y '
  588.               LESOBL(5)='V1Z '
  589.               LESOBL(6)='V2X '
  590.               LESOBL(7)='V2Y '
  591.               LESOBL(8)='V2Z '
  592.              ELSE
  593.               NBROBL=4
  594.               SEGINI NOMID
  595.               LESOBL(1)='PERM'
  596.               LESOBL(2)='VISC'
  597.               LESOBL(3)='V1X '
  598.               LESOBL(4)='V1Y '
  599.              ENDIF
  600.             ENDIF
  601. *
  602.          NMATR=NBROBL
  603.          NMATF=NBRFAC
  604.          NMATT = NMATR+NMATF
  605.          MOMATR=NOMID
  606. *
  607.          NBTYPE=1
  608.          SEGINI NOTYPE
  609.          MOTYPE=NOTYPE
  610.          TYPE(1)='REAL*8'
  611. *
  612.          CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT)
  613.          SEGSUP NOTYPE
  614.          IF (IERR.NE.0) GOTO 9992
  615.          IF (ISUP.EQ.1) THEN
  616.             CALL VALCHE(IVAMAT,NMATR,IPMINT,IPPORE,MOMATR,MELE)
  617.          ENDIF
  618. C
  619. C  TRAITEMENT DES CHAMPS DE CARACTERISTIQUES                   *
  620. C
  621.            NBROBL=0
  622.            NBRFAC=0
  623.            MOCARA=0
  624.            IVECT=0
  625. *
  626. * EPAISSEUR DANS LE CAS DES CONTRAINTES PLANES
  627. *
  628.          IF(IFOUR.EQ.-2.AND.((MELE.GE.79.AND.MELE.LE.83)
  629.      &                 .OR.(MELE.GE.173.AND.MELE.LE.182)))THEN
  630. *
  631.            NBROBL=0
  632.            NBRFAC=1
  633.            SEGINI NOMID
  634.            MOCARA=NOMID
  635.            LESFAC(1)='DIM3'
  636. *
  637.            NBTYPE=1
  638.            SEGINI NOTYPE
  639.            TYPE(1)='REAL*8'
  640.          ENDIF
  641. *
  642.          NCARA=NBROBL
  643.          NCARF=NBRFAC
  644.          NCARR=NCARA+NCARF
  645. *
  646.          IF (MOCARA.NE.0)  THEN
  647.            MOTYPE=NOTYPE
  648.                CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOCARA,MOTYPE,1,INFOS,3,
  649.      $                     IVACAR)
  650.                SEGSUP NOTYPE
  651.                IF (IERR.NE.0) GOTO 9992
  652. *
  653.            IF (ISUP.EQ.1) THEN
  654.               CALL VALCHE(IVACAR,NCARR,IPMINT,IPPORE,MOCARA,MELE)
  655.            ENDIF
  656.          ENDIF
  657. *
  658. *_____________________________________________________________________
  659. *
  660. *     NUMERO DES ETIQUETTES      :
  661. *     ETIQUETTES DE 1 A 98 POUR TRAITEMENT SPECIFIQUE A L ELEMENT
  662. *     DANS LA ZONE SPECIFIQUE A CHAQUE ELEMENT COMMENCANT PAR :
  663. *     5  CONTINUE
  664. *     ELEMENT 5   ETIQUETTES 1005 2005 3005 4005   ...
  665. *     44 CONTINUE
  666. *     ELEMENT 44  ETIQUETTES 1044 2044 3044 4044   ...
  667. *_____________________________________________________________________
  668.       GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  669.      1      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  670.      2      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  671.      3      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,79,79,
  672.      4      79,79,79,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  673.      5      99,99,99,99,99,99,99,80,80,80),MELE
  674. *
  675.       IF(MELE.GE.173.AND.MELE.LE.182) GO TO 173
  676.       IF(MELE.GE.185.AND.MELE.LE.190) GO TO 185
  677. *
  678.    99 CONTINUE
  679.       SEGSUP xMATRI
  680.       MOTERR(1:4)=NOMTP(MELE)
  681.       MOTERR(5:12)='PERMEABI'
  682.       CALL ERREUR(86)
  683.       GOTO 9990
  684. *_______________________________________________________________________
  685. *
  686. *     MILIEU POREUX
  687. *_______________________________________________________________________
  688. *
  689.   79  CONTINUE
  690. *
  691. *     Pour ces elements NBBB : Nb de noeuds
  692. *                       NBNO : Nb de fonctions de forme
  693. *
  694.       DIM3=1.D0
  695.       NBNO=IPORE
  696.       NBBB=NBNN
  697.       NSTN=1
  698.       LRN=NBNO-NBBB
  699.       NSTB=2
  700.       IF(IFOUR.GT.0) NSTB=3
  701. *
  702. *     CAS NON ISOTROPES
  703. *     RECUPERATION DES FONCTIONS DE FORME ET LEURS DERIVEES
  704. *     AU CENTRE DE L'ELEMENT POUR LE CALCUL DES AXES LOCAUX
  705. *
  706.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  707.         CALL RESHPT(1,NBNO,IELE,MELE,0,MINT,IRT1)
  708.         MINTE2=MINT
  709.         SEGACT MINTE2
  710.         SEGINI WRK4
  711.       ENDIF
  712. *
  713.       SEGINI WRK1,WRK2,WRK3
  714.       I195=0
  715.       I259=0
  716.       I367=0
  717.       DO 3079  IB=1,NBELEM
  718. *
  719. *        On cherche les coordonnees des noeuds de l'element IB
  720. *
  721.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  722. *
  723. *       calcul des axes locaux dans les cas non isotropes
  724. *
  725.          IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  726.            NBSH=MINTE2.SHPTOT(/2)
  727.            CALL RLOCAL (XE,MINTE2.SHPTOT,NBSH,NBNN,TXR)
  728.         if (nbsh.eq.-1) then
  729.          call erreur(525)
  730.          return
  731.         endif
  732.          ENDIF
  733. *
  734.          CALL ZERO (REL,LRE,LRE)
  735. *
  736. *       boucle sur les points de Gauss
  737. *
  738.          ISDJC=0
  739.          DO 4079  IGAU=1,NBPGAU
  740.  
  741. *        PRINT *, ' POINT DE GAUSS ',IGAU
  742. C
  743. C   RECUPERATION DE L'EPAISSEUR
  744. C
  745.             IF (IFOUR.EQ.-2)THEN
  746.              MPTVAL=IVACAR
  747.              IF (IVACAR.NE.0) THEN
  748.               MELVAL=IVAL(1)
  749.               IF (MELVAL.NE.0) THEN
  750.                 IGMN=MIN(IGAU,VELCHE(/1))
  751.                 IBMN=MIN(IB,VELCHE(/2))
  752.                 DIM3=VELCHE(IGMN,IBMN)
  753.               ELSE
  754.                 DIM3=1.D0
  755.               ENDIF
  756.              ENDIF
  757.             ENDIF
  758. *
  759.             CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  760.      &               DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,2)
  761.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  762.             IF(DJAC.EQ.0.) I259 =IB
  763.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  764.  
  765. *           PRINT *,' MATRICE B LIGNE PAR LIGNE '
  766. *            DO 3367 IPZ = 1,NSTB
  767. **           PRINT *,' LIGNE ',IPZ
  768. *            WRITE(6,3368) (BGENE(IPZ,JPZ), JPZ=1,LRE)
  769. *3368        FORMAT(8(1X,1PE10.3)/)
  770. *3367        CONTINUE
  771.  
  772.             EREF=1.D0
  773.             MPTVAL=IVAMAT
  774. *
  775. *           le cas isotrope
  776. *
  777.             IF (MATE.EQ.1) THEN
  778.               MELVAL=IVAL(1)
  779.               IGMN=MIN(IGAU,VELCHE(/1))
  780.               IBMN=MIN(IB  ,VELCHE(/2))
  781.               XK  =VELCHE(IGMN,IBMN)
  782. *
  783.               MELVAL=IVAL(2)
  784.               IGMN=MIN(IGAU,VELCHE(/1))
  785.               IBMN=MIN(IB  ,VELCHE(/2))
  786.               XMU =VELCHE(IGMN,IBMN)
  787.               IF(XMU.EQ.0.D0) THEN
  788.                 I367=IB
  789.                 GO TO 4079
  790.               ENDIF
  791.               COMJAC=DJAC*EREF*EREF*XK/XMU
  792.               DO 4279 I=1,LRN
  793.                  DO 4279 J=1,I
  794.                     DO 4279 K=1,NSTB
  795.                       REL(I,J)=REL(I,J)+COMJAC*BGENE(K,I)*BGENE(K,J)
  796.  4279         CONTINUE
  797. *
  798. *           les cas non isotropes
  799. *
  800.             ELSE IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  801. *
  802.               DO 4379 IM=1,NMATT
  803.                 IF (IVAL(IM).NE.0) THEN
  804.                   MELVAL=IVAL(IM)
  805.                   IBMN=MIN(IB  ,VELCHE(/2))
  806.                   IGMN=MIN(IGAU,VELCHE(/1))
  807.                   VALMAT(IM)=VELCHE(IGMN,IBMN)
  808.                 ELSE
  809.                   VALMAT(IM)=0.D0
  810.                 ENDIF
  811.  4379         CONTINUE
  812. *
  813.               CALL PERMAO(WRK4,IFOUR,MATE,EREF,KERRE)
  814.               IF(KERRE.EQ.1) GO TO 99
  815.               IF(KERRE.EQ.2) THEN
  816.                 I367=IB
  817.                 GO TO 4079
  818.               ENDIF
  819. *
  820.               CALL BDBST(BGENE,DJAC,PMAT,LRE,NSTB,REL)
  821. *
  822. *           les cas non prevus
  823. *
  824.             ELSE
  825.               GO TO 99
  826.             ENDIF
  827. *
  828.  4079    CONTINUE
  829.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  830. *         SEGINI XMATRI
  831. *         IMATTT(IB)=XMATRI
  832. *
  833. *        Remplissage de XMATRI
  834. *
  835.          CALL REMPMT(REL,LRE,RE(1,1,ib))
  836. *         SEGDES XMATRI
  837.  3079 CONTINUE
  838. *
  839.       SEGSUP WRK1,WRK2,WRK3
  840.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  841.         SEGDES MINTE2
  842.         SEGSUP WRK4
  843.       ENDIF
  844. *
  845.       IF(I195.NE.0) THEN
  846.          INTERR(1)=I195
  847.          CALL ERREUR(195)
  848.          GOTO 9990
  849.       ELSE IF(I259.NE.0) THEN
  850.          INTERR(1)=I259
  851.          CALL ERREUR(259)
  852.          GOTO 9990
  853.       ELSE IF(I367.NE.0) THEN
  854.          INTERR(1)=I367
  855.          CALL ERREUR(367)
  856.          GOTO 9990
  857.       ENDIF
  858. *
  859.       SEGDES xMATRI
  860. *
  861.       GOTO 610
  862. *_______________________________________________________________________
  863. *
  864. *     JOINTS EN FORMULATION MILIEUX POREUX
  865. *_______________________________________________________________________
  866. *
  867.   80  CONTINUE
  868. *
  869. *     Pour ces elements NBBB : Nb de noeuds
  870. *                       NBNO : Nb de fonctions de forme
  871. *
  872.       NBNO=IPORE
  873.       NBBB=NBNN
  874.       NSTN=1
  875.       LPP=(NBNO-NBBB)*3/2
  876.       LRN=LPP
  877.       NSTB=2
  878.       IF(IFOUR.EQ.2) NSTB=3
  879.  
  880. *     PRINT *,' NBNO=', NBNO
  881. *     PRINT *,' NBBB=', NBBB
  882. *     PRINT *,' NSTN=', NSTN
  883. *     PRINT *,' LRN =', LRN
  884. *     PRINT *,' LRE =', LRE
  885. *     PRINT *,' NSTB =', NSTB
  886. *
  887.       SEGINI WRK1,WRK2,WRK3,WRK5
  888.       I195=0
  889.       I259=0
  890.       I367=0
  891.       DO 3080  IB=1,NBELEM
  892. *
  893. *        On cherche les coordonnees des noeuds de l'element IB
  894. *
  895.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  896. *
  897. *       calcul des axes locaux
  898. *
  899.       CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  900. *
  901.          CALL ZERO (REL,LRE,LRE)
  902. *
  903.       CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE)
  904. *
  905. *       boucle sur les points de Gauss
  906. *
  907.          ISDJC=0
  908.          DO 4080  IGAU=1,NBPGAU
  909. *
  910.       CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,XE,XEL,
  911.      .            SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,3)
  912.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  913.             IF(DJAC.EQ.0.) I259 =IB
  914.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  915. *
  916.             EREF=1.D0
  917.             MPTVAL=IVAMAT
  918. *
  919. *           le cas isotrope (le seul)
  920. *
  921.               MELVAL=IVAL(1)
  922.               IGMN=MIN(IGAU,VELCHE(/1))
  923.               IBMN=MIN(IB  ,VELCHE(/2))
  924.               XKT =VELCHE(IGMN,IBMN)
  925. *
  926.               MELVAL=IVAL(2)
  927.               IGMN=MIN(IGAU,VELCHE(/1))
  928.               IBMN=MIN(IB  ,VELCHE(/2))
  929.               XKNH =VELCHE(IGMN,IBMN)
  930. *
  931.               MELVAL=IVAL(3)
  932.               IGMN=MIN(IGAU,VELCHE(/1))
  933.               IBMN=MIN(IB  ,VELCHE(/2))
  934.               XKNB =VELCHE(IGMN,IBMN)
  935. *
  936.               MELVAL=IVAL(4)
  937.               IGMN=MIN(IGAU,VELCHE(/1))
  938.               IBMN=MIN(IB  ,VELCHE(/2))
  939.               XMU =VELCHE(IGMN,IBMN)
  940.               IF(XMU.EQ.0.D0) THEN
  941.                 I367=IB
  942.                 GO TO 4080
  943.               ENDIF
  944.               COMJAT=DJAC*EREF*EREF*XKT/XMU
  945.               COMJNH=DJAC*EREF*EREF*XKNH/XMU
  946.               COMJNB=DJAC*EREF*EREF*XKNB/XMU
  947.               DO 4280 I=1,LRN
  948.                  DO 4280 J=1,I
  949.                   REL(I,J)=REL(I,J)+COMJAT*BGENE(1,I)*BGENE(1,J)
  950.      .                     *XNTT(I)*XNTT(J)
  951.      .           +COMJNH*XGENE(1,I)*XGENE(1,J)*XNTH(I,J)
  952.      .           +COMJNB*XGENE(1,I)*XGENE(1,J)*XNTB(I,J)
  953.                    IF(IFOUR.EQ.2)THEN
  954.                   REL(I,J)=REL(I,J)+COMJAT*BGENE(2,I)*BGENE(2,J)
  955.      .                     *XNTT(I)*XNTT(J)
  956.                    ENDIF
  957.  4280        CONTINUE
  958. *
  959.  4080    CONTINUE
  960.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  961. *         SEGINI XMATRI
  962. *         IMATTT(IB)=XMATRI
  963. *
  964. *        Remplissage de XMATRI
  965. *
  966.          CALL REMPMT(REL,LRE,RE(1,1,ib))
  967. *         SEGDES XMATRI
  968.  3080 CONTINUE
  969. *
  970.       SEGSUP WRK1,WRK2,WRK3,WRK5
  971. *
  972.       IF(I195.NE.0) THEN
  973.          INTERR(1)=I195
  974.          CALL ERREUR(195)
  975.          GOTO 9990
  976.       ELSE IF(I259.NE.0) THEN
  977.          INTERR(1)=I259
  978.          CALL ERREUR(259)
  979.          GOTO 9990
  980.       ELSE IF(I367.NE.0) THEN
  981.          INTERR(1)=I367
  982.          CALL ERREUR(367)
  983.          GOTO 9990
  984.       ENDIF
  985. *
  986.       SEGDES xMATRI
  987. *
  988.       GOTO 610
  989. *_______________________________________________________________________
  990. *
  991. *     MILIEU POREUX  - SUITE
  992. *_______________________________________________________________________
  993. *
  994.  173  CONTINUE
  995. *
  996. *     Pour ces elements NBBB : Nb de noeuds
  997. *                       NBNO : Nb de fonctions de forme
  998. *
  999.       DIM3=1.D0
  1000.       NBNO=IPORE
  1001.       NBBB=NBNN
  1002.       NSTN=IDECAP
  1003.       LPP=NBNO-NBBB
  1004.       LRN=IDECAP*LPP
  1005.       NSTBE=2
  1006.       IF(IFOUR.GT.0) NSTBE=3
  1007.       NSTB=NSTBE*IDECAP
  1008.  
  1009. *     PRINT *,'NSTBE=',NSTBE
  1010. *     PRINT *,'NSTB=',NSTB
  1011. *     PRINT *,'IDECAP=',IDECAP
  1012. *     PRINT *,'LRE   =',LRE
  1013.  
  1014. *
  1015. *     CAS NON ISOTROPES
  1016. *     NON PREVU
  1017. *
  1018.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  1019.         CALL ERREUR(251)
  1020.         GO TO 9990
  1021.       ENDIF
  1022. *
  1023.       NSTPK= NSTB
  1024.       SEGINI WRK1,WRK2,WRK3,WRK6
  1025.       I195=0
  1026.       I259=0
  1027.       I367=0
  1028.       DO 3173  IB=1,NBELEM
  1029.  
  1030. *     PRINT *,'ELEMENT ' , IB
  1031. *
  1032. *        On cherche les coordonnees des noeuds de l'element IB
  1033. *
  1034.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1035. *
  1036.          CALL ZERO (REL,LRE,LRE)
  1037. *
  1038. *       boucle sur les points de Gauss
  1039. *
  1040.          ISDJC=0
  1041.          DO 4173  IGAU=1,NBPGAU
  1042.  
  1043. *        PRINT *, ' POINT DE GAUSS ',IGAU
  1044.  
  1045.  
  1046. C
  1047. C   RECUPERATION DE L'EPAISSEUR
  1048. C
  1049.             IF (IFOUR.EQ.-2)THEN
  1050.              MPTVAL=IVACAR
  1051.              IF (IVACAR.NE.0) THEN
  1052.               MELVAL=IVAL(1)
  1053.               IF (MELVAL.NE.0) THEN
  1054.                 IGMN=MIN(IGAU,VELCHE(/1))
  1055.                 IBMN=MIN(IB,VELCHE(/2))
  1056.                 DIM3=VELCHE(IGMN,IBMN)
  1057.               ELSE
  1058.                 DIM3=1.D0
  1059.               ENDIF
  1060.              ENDIF
  1061.             ENDIF
  1062. *
  1063.             LHOO=NSTB
  1064.             CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  1065.      &         DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOO,2)
  1066.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  1067.             IF(DJAC.EQ.0.) I259 =IB
  1068.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  1069.  
  1070. *           PRINT *,' MATRICE B LIGNE PAR LIGNE '
  1071. *           DO 1367 IPZ = 1,NSTB
  1072. *           PRINT *,' LIGNE ',IPZ
  1073. *           WRITE(6,1368) (BGENE(IPZ,JPZ), JPZ=1,LRE)
  1074. *1368        FORMAT(8(1X,1PE10.3)/)
  1075. *1367        CONTINUE
  1076.  
  1077. *
  1078.             EREF=1.D0
  1079.             MPTVAL=IVAMAT
  1080. *
  1081. *           le cas isotrope
  1082. *
  1083.             IF (MATE.EQ.1) THEN
  1084.  
  1085.               ICO=1
  1086.               DO 1731 ICD = 1,IDECAP
  1087.                 ICDA =(ICD -1) * NSTBE
  1088.                 DO 1732 JCD = 1,IDECAP
  1089.                   JCDA =(JCD -1) * NSTBE
  1090.                   MELVAL=IVAL(ICO)
  1091.                   IGMN=MIN(IGAU,VELCHE(/1))
  1092.                   IBMN=MIN(IB  ,VELCHE(/2))
  1093.                   DO 1733 KCD = 1,NSTBE
  1094.                     PKK(ICDA+KCD,JCDA+KCD) =VELCHE(IGMN,IBMN)
  1095. 1733            CONTINUE
  1096.                  ICO=ICO+1
  1097. 1732            CONTINUE
  1098. 1731          CONTINUE
  1099. *
  1100.  
  1101. *             PRINT *,' MATRICE PKK'
  1102. *              IF (IDECAP.EQ.2) THEN
  1103. *              WRITE (6,1342) ((PKK(I,J),J=1,NSTB),I=1,NSTB)
  1104. *1342          FORMAT(4(1X,1PE12.5)/)
  1105. *
  1106. *              ELSE IF (IDECAP.EQ.3) THEN
  1107. *              WRITE (6,1343) ((PKK(I,J),J=1,NSTB),I=1,NSTB)
  1108. *1343          FORMAT(6(1X,1PE12.5)/)
  1109. *              ENDIF
  1110.  
  1111.               COMJAC=DJAC*EREF*EREF
  1112.               CALL BDBSTS(BGENE,COMJAC,PKK,LRE,NSTB,REL)
  1113. *
  1114. *           les cas non prevus
  1115. *
  1116.             ELSE
  1117.               GO TO 99
  1118.             ENDIF
  1119. *
  1120.  4173    CONTINUE
  1121.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  1122. *         SEGINI XMATRI
  1123. *         IMATTT(IB)=XMATRI
  1124. *
  1125. *        Remplissage de XMATRI
  1126. *
  1127.          CALL REMPMS(REL,LRE,RE(1,1,ib))
  1128. *        SEGDES XMATRI
  1129.  3173 CONTINUE
  1130. *
  1131.       SEGSUP WRK1,WRK2,WRK3,WRK6
  1132.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  1133.         SEGDES MINTE2
  1134.         SEGSUP WRK4
  1135.       ENDIF
  1136. *
  1137.       IF(I195.NE.0) THEN
  1138.          INTERR(1)=I195
  1139.          CALL ERREUR(195)
  1140.          GOTO 9990
  1141.       ELSE IF(I259.NE.0) THEN
  1142.          INTERR(1)=I259
  1143.          CALL ERREUR(259)
  1144.          GOTO 9990
  1145.       ELSE IF(I367.NE.0) THEN
  1146.          INTERR(1)=I367
  1147.          CALL ERREUR(367)
  1148.          GOTO 9990
  1149.       ENDIF
  1150. *
  1151.       SEGDES xMATRI
  1152. *
  1153.       GOTO 610
  1154. *_______________________________________________________________________
  1155. *
  1156. *     JOINTS EN FORMULATION MILIEUX POREUX - SUITE
  1157. *_______________________________________________________________________
  1158. *
  1159.  185  CONTINUE
  1160. *
  1161. *     Pour ces elements NBBB : Nb de noeuds
  1162. *                       NBNO : Nb de fonctions de forme
  1163. *
  1164.       NBNO=IPORE
  1165.       NBBB=NBNN
  1166.       NSTN=IDECAP
  1167.       LPP=(NBNO-NBBB)*3/2
  1168.       LRN=IDECAP*LPP
  1169.       NSTBE=2
  1170.       IF(IFOUR.EQ.2) NSTBE=3
  1171.       NSTB=NSTBE*IDECAP
  1172.       NSTPKE=3
  1173.       NSTPK=NSTPKE*IDECAP
  1174.  
  1175. *     PRINT *,' NBNO=', NBNO
  1176. *     PRINT *,' NBBB=', NBBB
  1177. *     PRINT *,' NSTN=', NSTN
  1178. *     PRINT *,' LPP =', LPP
  1179. *     PRINT *,' LRN =', LRN
  1180. *     PRINT *,' LRE =', LRE
  1181. *     PRINT *,' NSTBE=', NSTBE
  1182. *     PRINT *,' NSTB =', NSTB
  1183. *     PRINT *,' NSTPKE =', NSTPKE
  1184. *     PRINT *,' NSTPK =', NSTPK
  1185. *
  1186.       SEGINI WRK1,WRK2,WRK3,WRK5,WRK6
  1187.       I195=0
  1188.       I259=0
  1189.       I367=0
  1190.       DO 3185  IB=1,NBELEM
  1191. *
  1192. *        On cherche les coordonnees des noeuds de l'element IB
  1193. *
  1194.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1195. *
  1196. *       calcul des axes locaux
  1197. *
  1198.       CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1199. *
  1200.          CALL ZERO (REL,LRE,LRE)
  1201. *
  1202.       CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE)
  1203. *
  1204. *       boucle sur les points de Gauss
  1205. *
  1206.       ISDJC=0
  1207.       DO 4185  IGAU=1,NBPGAU
  1208. *
  1209.       LHOO=NSTB
  1210.       CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOO,NSTN,XE,XEL,
  1211.      .     SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,3)
  1212.       IF(DJAC.LT.0.) ISDJC=ISDJC+1
  1213.       IF(DJAC.EQ.0.) I259 =IB
  1214.       DJAC=ABS(DJAC)*POIGAU(IGAU)
  1215. *
  1216.       EREF=1.D0
  1217.       MPTVAL=IVAMAT
  1218. *
  1219. *           le cas isotrope (le seul)
  1220. *
  1221.       IF(MATE.EQ.1) THEN
  1222.  
  1223.         FAC = DJAC*EREF*EREF
  1224.         IE=0
  1225.         DO 2185 IPR=1,IDECAP
  1226.           IPR1 = (IPR-1) * NSTPKE
  1227.           DO 2185 JPR=1,IDECAP
  1228.             JPR1 = (JPR-1) * NSTPKE
  1229.             DO 2186 I=1,NSTPKE
  1230.               II = I + IPR1
  1231.               JJ = I + JPR1
  1232.               IE=IE+1
  1233.               MELVAL=IVAL(IE)
  1234.               IGMN=MIN(IGAU,VELCHE(/1))
  1235.               IBMN=MIN(IB  ,VELCHE(/2))
  1236.               PKK(II,JJ)=VELCHE(IGMN,IBMN)*FAC
  1237. 2186        CONTINUE
  1238. 2185    CONTINUE
  1239.  
  1240. *
  1241.         DO 8985 IPR=1,IDECAP
  1242.           IPR1 = (IPR-1)*NSTPKE
  1243.           IPR2 = 2*IPR
  1244.           IPPDEC=(IPR-1)*LPP
  1245.           IRRDEC=(IPR-1)*NBBB
  1246.           DO 8985 JPR=1,IDECAP
  1247.             JPR1 = (JPR-1)*NSTPKE
  1248.             JPR2 = 2*JPR
  1249.             JPPDEC=(JPR-1)*LPP
  1250.             JRRDEC=(JPR-1)*NBBB
  1251. *
  1252.             COMJAT=PKK(IPR1+1,JPR1+1)
  1253.             COMJNH=PKK(IPR1+2,JPR1+2)
  1254.             COMJNB=PKK(IPR1+3,JPR1+3)
  1255. *
  1256.             IF(IFOUR.LE.0) THEN
  1257.               DO 4285 I=1,LPP
  1258.                 II =I+IPPDEC
  1259.                 IR =I+IRRDEC
  1260.                 DO 4285 J=1,LPP
  1261.                   JJ =J+JPPDEC
  1262.                   JR =J+JRRDEC
  1263.                   REL(IR,JR)=REL(IR,JR)
  1264.      .       +COMJAT*BGENE(IPR,II)*BGENE(JPR,JJ)
  1265.      .                   *XNTT(I)*XNTT(J)
  1266.      .       +COMJNH*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTH(I,J)
  1267.      .       +COMJNB*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTB(I,J)
  1268.  4285     CONTINUE
  1269. *
  1270.             ELSE
  1271.               DO 4385 I=1,LPP
  1272.                 II =I+IPPDEC
  1273.                 DO 4385 J=1,LPP
  1274.                   JJ =J+JPPDEC
  1275.                   REL(IR,JR)=REL(IR,JR)
  1276.      .            +COMJAT*XNTT(I)*XNTT(J)*
  1277.      .             (BGENE(IPR2-1,II)*BGENE(JPR2-1,JJ)
  1278.      .              + BGENE(IPR2,II)*BGENE(JPR2,JJ))
  1279.      .           +COMJNH*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTH(I,J)
  1280.      .           +COMJNB*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTB(I,J)
  1281.  4385         CONTINUE
  1282.           ENDIF
  1283.  
  1284.  8985   CONTINUE
  1285. *
  1286.       ELSE
  1287.         GO TO 9990
  1288.       ENDIF
  1289.  
  1290. *
  1291.  4185    CONTINUE
  1292.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  1293. *
  1294. *        Remplissage de XMATRI
  1295. *
  1296.          CALL REMPMS(REL,LRE,RE(1,1,ib))
  1297.  
  1298.  3185 CONTINUE
  1299. *
  1300.       SEGSUP WRK1,WRK2,WRK3,WRK5,WRK6
  1301. *
  1302.       IF(I195.NE.0) THEN
  1303.          INTERR(1)=I195
  1304.          CALL ERREUR(195)
  1305.          GOTO 9990
  1306.       ELSE IF(I259.NE.0) THEN
  1307.          INTERR(1)=I259
  1308.          CALL ERREUR(259)
  1309.          GOTO 9990
  1310.       ELSE IF(I367.NE.0) THEN
  1311.          INTERR(1)=I367
  1312.          CALL ERREUR(367)
  1313.          GOTO 9990
  1314.       ENDIF
  1315. *
  1316.       GOTO 610
  1317. *
  1318. *     Desactivation des segment propres a la geometrie ISOUS
  1319. *
  1320.   610 CONTINUE
  1321. *
  1322.       SEGDES MELEME
  1323.       SEGDES MINTE
  1324. *
  1325.       IF (ISUP.EQ.1) THEN
  1326.         CALL DTMVAL(IVAMAT,3)
  1327.       ELSE
  1328.         CALL DTMVAL(IVAMAT,1)
  1329.       ENDIF
  1330. *
  1331.       NOMID=MOMATR
  1332.       SEGSUP NOMID
  1333.       IF (lsupdp) THEN
  1334.         NOMID=MODEPL
  1335.         SEGSUP,NOMID
  1336.       ENDIF
  1337.       IF (lsupfo) THEN
  1338.         NOMID=MOFORC
  1339.         SEGSUP,NOMID
  1340.       ENDIF
  1341. *
  1342.       COERIG(ISORI)  = 1.D0
  1343.       IRIGEL(1,ISORI)=IPMAIL
  1344.       IRIGEL(2,ISORI)=0
  1345.       IRIGEL(3,ISORI)=IPDES
  1346.       IRIGEL(4,ISORI)=xMATRI
  1347.       IRIGEL(5,ISORI)=NHRM
  1348.       IRIGEL(6,ISORI)=0
  1349.       IRIGEL(7,ISORI)=LASYM
  1350.       xmatri.symre=lasym
  1351.       SEGDES xMATRI
  1352.       IRIGEL(8,ISORI)=0
  1353.  
  1354. *  Fin de la boucle de PARTITIONNEMENT du segment XMATRI
  1355.  5000 CONTINUE
  1356.  
  1357.       SEGDES IMODEL
  1358.  
  1359.   500 CONTINUE
  1360.       IF(ISORI.NE.NRIGEL) GO TO 9999
  1361.       SEGDES MRIGID
  1362. C      SEGSUP MMODEL
  1363.       IRET = 1
  1364.       IPRIGI = MRIGID
  1365.       RETURN
  1366. *
  1367. *     Erreur dans une sous zone desactivation et retour
  1368. *
  1369.  9990 CONTINUE
  1370.       IF (ISUP.EQ.1) THEN
  1371.         CALL DTMVAL(IVAMAT,3)
  1372.       ELSE
  1373.         CALL DTMVAL(IVAMAT,1)
  1374.       ENDIF
  1375. *
  1376.  9992 CONTINUE
  1377.  
  1378.       SEGSUP DESCR
  1379.       SEGSUP xMATRI
  1380.       SEGDES MELEME
  1381.       SEGDES MINTE
  1382. *
  1383.       NOMID=MOMATR
  1384.       SEGSUP NOMID
  1385.  9999 CONTINUE
  1386.       SEGSUP MRIGID
  1387.       IRET = 0
  1388.       IPRIGI = 0
  1389.  
  1390.  9991 CONTINUE
  1391.  
  1392.       MMODEL = IPMODL
  1393.       DO isous = 1, NSOUS
  1394.         IMODEL = KMODEL(isous)
  1395.         SEGDES,IMODEL
  1396.       ENDDO
  1397. C      SEGSUP,MMODEL
  1398.  
  1399.       RETURN
  1400.       END
  1401.  
  1402.  
  1403.  

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