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

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