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

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