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    PV090527  26/04/30    21:15:58     12529          
  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.          rigrel=0
  404.          SEGINI xMATRI
  405. *
  406. *        Verification de la presence des composantes pour le materiau
  407. *
  408.          NBROBL=0
  409.          NBRFAC=0
  410. *                                         cas isotrope
  411.             IF (MATE.EQ.1) THEN
  412. *
  413.               IF (MELE.GE.79.AND.MELE.LE.83) THEN
  414.                 NBROBL=2
  415.                 SEGINI NOMID
  416.                 LESOBL(1)='PERM'
  417.                 LESOBL(2)='VISC'
  418.               ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN
  419.                 NBROBL=4
  420.                 SEGINI NOMID
  421.                 LESOBL(1)='PERT'
  422.                 LESOBL(2)='PERH'
  423.                 LESOBL(3)='PERB'
  424.                 LESOBL(4)='VISC'
  425.               ELSE IF (MELE.GE.173.AND.MELE.LE.177) THEN
  426.                 NBROBL=4
  427.                 SEGINI NOMID
  428.                 LESOBL(1)='PK11'
  429.                 LESOBL(2)='PK12'
  430.                 LESOBL(3)='PK21'
  431.                 LESOBL(4)='PK22'
  432.               ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN
  433.                 NBROBL=9
  434.                 SEGINI NOMID
  435.                 LESOBL(1)='PK11'
  436.                 LESOBL(2)='PK12'
  437.                 LESOBL(3)='PK13'
  438.                 LESOBL(4)='PK21'
  439.                 LESOBL(5)='PK22'
  440.                 LESOBL(6)='PK23'
  441.                 LESOBL(7)='PK31'
  442.                 LESOBL(8)='PK32'
  443.                 LESOBL(9)='PK33'
  444.               ELSE IF (MELE.GE.185.AND.MELE.LE.187) THEN
  445.                 NBROBL=12
  446.                 SEGINI NOMID
  447.                 LESOBL(1)='PT11'
  448.                 LESOBL(2)='PH11'
  449.                 LESOBL(3)='PB11'
  450.                 LESOBL(4)='PT12'
  451.                 LESOBL(5)='PH12'
  452.                 LESOBL(6)='PB12'
  453.                 LESOBL(7)='PT21'
  454.                 LESOBL(8)='PH21'
  455.                 LESOBL(9)='PB21'
  456.                 LESOBL(10)='PT22'
  457.                 LESOBL(11)='PH22'
  458.                 LESOBL(12)='PB22'
  459.               ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN
  460.                 NBROBL=27
  461.                 SEGINI NOMID
  462.                 LESOBL(1)='PT11'
  463.                 LESOBL(2)='PH11'
  464.                 LESOBL(3)='PB11'
  465.                 LESOBL(4)='PT12'
  466.                 LESOBL(5)='PH12'
  467.                 LESOBL(6)='PB12'
  468.                 LESOBL(7)='PT13'
  469.                 LESOBL(8)='PH13'
  470.                 LESOBL(9)='PB13'
  471.                 LESOBL(10)='PT21'
  472.                 LESOBL(11)='PH21'
  473.                 LESOBL(12)='PB21'
  474.                 LESOBL(13)='PT22'
  475.                 LESOBL(14)='PH22'
  476.                 LESOBL(15)='PB22'
  477.                 LESOBL(16)='PT23'
  478.                 LESOBL(17)='PH23'
  479.                 LESOBL(18)='PB23'
  480.                 LESOBL(19)='PT31'
  481.                 LESOBL(20)='PH31'
  482.                 LESOBL(21)='PB31'
  483.                 LESOBL(22)='PT32'
  484.                 LESOBL(23)='PH32'
  485.                 LESOBL(24)='PB32'
  486.                 LESOBL(25)='PT33'
  487.                 LESOBL(26)='PH33'
  488.                 LESOBL(27)='PB33'
  489.               ENDIF
  490. *                                          cas orthotrope
  491.             ELSE IF (MATE.EQ.2) THEN
  492.               IF (IDIM.EQ.3) THEN
  493.                 NBROBL=10
  494.                 SEGINI NOMID
  495.                 LESOBL(1)='PER1'
  496.                 LESOBL(2)='PER2'
  497.                 LESOBL(3)='PER3'
  498.                 LESOBL(4)='VISC'
  499.                 LESOBL(5)='V1X '
  500.                 LESOBL(6)='V1Y '
  501.                 LESOBL(7)='V1Z '
  502.                 LESOBL(8)='V2X '
  503.                 LESOBL(9)='V2Y '
  504.                 LESOBL(10)='V2Z '
  505.               ELSE IF(IDIM.EQ.2) THEN
  506.                 IF (IFOUR.LE.0) THEN
  507.                   NBROBL=5
  508.                   SEGINI NOMID
  509.                   LESOBL(1)='PER1'
  510.                   LESOBL(2)='PER2'
  511.                   LESOBL(3)='VISC'
  512.                   LESOBL(4)='V1X '
  513.                   LESOBL(5)='V1Y '
  514.                 ELSE IF (IFOUR.EQ.1) THEN
  515.                   NBROBL=6
  516.                   SEGINI NOMID
  517.                   LESOBL(1)='PER1'
  518.                   LESOBL(2)='PER2'
  519.                   LESOBL(3)='PER3'
  520.                   LESOBL(4)='VISC'
  521.                   LESOBL(5)='V1X '
  522.                   LESOBL(6)='V1Y '
  523.                 ENDIF
  524.               ENDIF
  525. *                                          cas anisotrope
  526.             ELSE IF (MATE.EQ.3)THEN
  527.               IF(IDIM.EQ.3)THEN
  528.                 NBROBL=13
  529.                 SEGINI NOMID
  530.                 LESOBL(1)='PER1'
  531.                 LESOBL(2)='PER2'
  532.                 LESOBL(3)='PER3'
  533.                 LESOBL(4)='PE12'
  534.                 LESOBL(5)='PE13'
  535.                 LESOBL(6)='PE23'
  536.                 LESOBL(7)='VISC'
  537.                 LESOBL(8)='V1X '
  538.                 LESOBL(9)='V1Y '
  539.                 LESOBL(10)='V1Z '
  540.                 LESOBL(11)='V2X '
  541.                 LESOBL(12)='V2Y '
  542.                 LESOBL(13)='V2Z '
  543.               ELSE IF (IDIM.EQ.2) THEN
  544.                 IF (IFOUR.LE.0) THEN
  545.                   NBROBL=6
  546.                   SEGINI NOMID
  547.                   LESOBL(1)='PER1'
  548.                   LESOBL(2)='PER2'
  549.                   LESOBL(3)='PE12'
  550.                   LESOBL(4)='VISC'
  551.                   LESOBL(5)='V1X '
  552.                   LESOBL(6)='V1Y '
  553.                 ELSE IF (IFOUR.EQ.1) THEN
  554.                   NBROBL=7
  555.                   SEGINI NOMID
  556.                   LESOBL(1)='PER1'
  557.                   LESOBL(2)='PER2'
  558.                   LESOBL(3)='PE12'
  559.                   LESOBL(4)='PER3'
  560.                   LESOBL(5)='VISC'
  561.                   LESOBL(6)='V1X '
  562.                   LESOBL(7)='V1Y '
  563.                 ENDIF
  564.               ENDIF
  565. *                                         cas unidirectionnel
  566.             ELSE IF (MATE.EQ.4) THEN
  567.              IF (IDIM.EQ.3) THEN
  568.               NBROBL=8
  569.               SEGINI NOMID
  570.               LESOBL(1)='PERM'
  571.               LESOBL(2)='VISC'
  572.               LESOBL(3)='V1X '
  573.               LESOBL(4)='V1Y '
  574.               LESOBL(5)='V1Z '
  575.               LESOBL(6)='V2X '
  576.               LESOBL(7)='V2Y '
  577.               LESOBL(8)='V2Z '
  578.              ELSE
  579.               NBROBL=4
  580.               SEGINI NOMID
  581.               LESOBL(1)='PERM'
  582.               LESOBL(2)='VISC'
  583.               LESOBL(3)='V1X '
  584.               LESOBL(4)='V1Y '
  585.              ENDIF
  586.             ENDIF
  587. *
  588.          NMATR=NBROBL
  589.          NMATF=NBRFAC
  590.          NMATT = NMATR+NMATF
  591.          MOMATR=NOMID
  592. *
  593.          NBTYPE=1
  594.          SEGINI NOTYPE
  595.          MOTYPE=NOTYPE
  596.          TYPE(1)='REAL*8'
  597. *
  598.          CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT)
  599.          SEGSUP NOTYPE
  600.          IF (IERR.NE.0) GOTO 9992
  601.          IF (ISUP.EQ.1) THEN
  602.             CALL VALCHE(IVAMAT,NMATR,IPMINT,IPPORE,MOMATR,MELE)
  603.          ENDIF
  604. C
  605. C  TRAITEMENT DES CHAMPS DE CARACTERISTIQUES                   *
  606. C
  607.            NBROBL=0
  608.            NBRFAC=0
  609.            MOCARA=0
  610.            IVECT=0
  611. *
  612. * EPAISSEUR DANS LE CAS DES CONTRAINTES PLANES
  613. *
  614.          IF(IFOUR.EQ.-2.AND.((MELE.GE.79.AND.MELE.LE.83)
  615.      &                 .OR.(MELE.GE.173.AND.MELE.LE.182)))THEN
  616. *
  617.            NBROBL=0
  618.            NBRFAC=1
  619.            SEGINI NOMID
  620.            MOCARA=NOMID
  621.            LESFAC(1)='DIM3'
  622. *
  623.            NBTYPE=1
  624.            SEGINI NOTYPE
  625.            TYPE(1)='REAL*8'
  626.          ENDIF
  627. *
  628.          NCARA=NBROBL
  629.          NCARF=NBRFAC
  630.          NCARR=NCARA+NCARF
  631. *
  632.          IF (MOCARA.NE.0)  THEN
  633.            MOTYPE=NOTYPE
  634.                CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOCARA,MOTYPE,1,INFOS,3,
  635.      $                     IVACAR)
  636.                SEGSUP NOTYPE
  637.                IF (IERR.NE.0) GOTO 9992
  638. *
  639.            IF (ISUP.EQ.1) THEN
  640.               CALL VALCHE(IVACAR,NCARR,IPMINT,IPPORE,MOCARA,MELE)
  641.            ENDIF
  642.          ENDIF
  643. *
  644. *_____________________________________________________________________
  645. *
  646. *     NUMERO DES ETIQUETTES      :
  647. *     ETIQUETTES DE 1 A 98 POUR TRAITEMENT SPECIFIQUE A L ELEMENT
  648. *     DANS LA ZONE SPECIFIQUE A CHAQUE ELEMENT COMMENCANT PAR :
  649. *     5  CONTINUE
  650. *     ELEMENT 5   ETIQUETTES 1005 2005 3005 4005   ...
  651. *     44 CONTINUE
  652. *     ELEMENT 44  ETIQUETTES 1044 2044 3044 4044   ...
  653. *_____________________________________________________________________
  654.       GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  655.      1      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  656.      2      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  657.      3      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,79,79,
  658.      4      79,79,79,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  659.      5      99,99,99,99,99,99,99,80,80,80),MELE
  660. *
  661.       IF(MELE.GE.173.AND.MELE.LE.182) GO TO 173
  662.       IF(MELE.GE.185.AND.MELE.LE.190) GO TO 185
  663. *
  664.    99 CONTINUE
  665.       SEGSUP xMATRI
  666.       MOTERR(1:4)=NOMTP(MELE)
  667.       MOTERR(5:12)='PERMEABI'
  668.       CALL ERREUR(86)
  669.       GOTO 9990
  670. *_______________________________________________________________________
  671. *
  672. *     MILIEU POREUX
  673. *_______________________________________________________________________
  674. *
  675.   79  CONTINUE
  676. *
  677. *     Pour ces elements NBBB : Nb de noeuds
  678. *                       NBNO : Nb de fonctions de forme
  679. *
  680.       DIM3=1.D0
  681.       NBNO=IPORE
  682.       NBBB=NBNN
  683.       NSTN=1
  684.       LRN=NBNO-NBBB
  685.       NSTB=2
  686.       IF(IFOUR.GT.0) NSTB=3
  687. *
  688. *     CAS NON ISOTROPES
  689. *     RECUPERATION DES FONCTIONS DE FORME ET LEURS DERIVEES
  690. *     AU CENTRE DE L'ELEMENT POUR LE CALCUL DES AXES LOCAUX
  691. *
  692.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  693.         CALL RESHPT(1,NBNO,IELE,MELE,0,MINT,IRT1)
  694.         MINTE2=MINT
  695.         SEGACT MINTE2
  696.         SEGINI WRK4
  697.       ENDIF
  698. *
  699.       SEGINI WRK1,WRK2,WRK3
  700.       I195=0
  701.       I259=0
  702.       I367=0
  703.       DO 3079  IB=1,NBELEM
  704. *
  705. *        On cherche les coordonnees des noeuds de l'element IB
  706. *
  707.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  708. *
  709. *       calcul des axes locaux dans les cas non isotropes
  710. *
  711.          IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  712.            NBSH=MINTE2.SHPTOT(/2)
  713.            CALL RLOCAL (XE,MINTE2.SHPTOT,NBSH,NBNN,TXR)
  714.         if (nbsh.eq.-1) then
  715.          call erreur(525)
  716.          return
  717.         endif
  718.          ENDIF
  719. *
  720.          CALL ZERO (REL,LRE,LRE)
  721. *
  722. *       boucle sur les points de Gauss
  723. *
  724.          ISDJC=0
  725.          DO 4079  IGAU=1,NBPGAU
  726.  
  727. *        PRINT *, ' POINT DE GAUSS ',IGAU
  728. C
  729. C   RECUPERATION DE L'EPAISSEUR
  730. C
  731.             IF (IFOUR.EQ.-2)THEN
  732.              MPTVAL=IVACAR
  733.              IF (IVACAR.NE.0) THEN
  734.               MELVAL=IVAL(1)
  735.               IF (MELVAL.NE.0) THEN
  736.                 IGMN=MIN(IGAU,VELCHE(/1))
  737.                 IBMN=MIN(IB,VELCHE(/2))
  738.                 DIM3=VELCHE(IGMN,IBMN)
  739.               ELSE
  740.                 DIM3=1.D0
  741.               ENDIF
  742.              ENDIF
  743.             ENDIF
  744. *
  745.             CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  746.      &               DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,2)
  747.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  748.             IF(DJAC.EQ.0.) I259 =IB
  749.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  750.  
  751. *           PRINT *,' MATRICE B LIGNE PAR LIGNE '
  752. *            DO 3367 IPZ = 1,NSTB
  753. **           PRINT *,' LIGNE ',IPZ
  754. *            WRITE(6,3368) (BGENE(IPZ,JPZ), JPZ=1,LRE)
  755. *3368        FORMAT(8(1X,1PE10.3)/)
  756. *3367        CONTINUE
  757.  
  758.             EREF=1.D0
  759.             MPTVAL=IVAMAT
  760. *
  761. *           le cas isotrope
  762. *
  763.             IF (MATE.EQ.1) THEN
  764.               MELVAL=IVAL(1)
  765.               IGMN=MIN(IGAU,VELCHE(/1))
  766.               IBMN=MIN(IB  ,VELCHE(/2))
  767.               XK  =VELCHE(IGMN,IBMN)
  768. *
  769.               MELVAL=IVAL(2)
  770.               IGMN=MIN(IGAU,VELCHE(/1))
  771.               IBMN=MIN(IB  ,VELCHE(/2))
  772.               XMU =VELCHE(IGMN,IBMN)
  773.               IF(XMU.EQ.0.D0) THEN
  774.                 I367=IB
  775.                 GO TO 4079
  776.               ENDIF
  777.               COMJAC=DJAC*EREF*EREF*XK/XMU
  778.               DO 4277 I=1,LRN
  779.                  DO 4278 J=1,I
  780.                     DO 4279 K=1,NSTB
  781.                       REL(I,J)=REL(I,J)+COMJAC*BGENE(K,I)*BGENE(K,J)
  782.  4279         CONTINUE
  783.  4278         CONTINUE
  784.  4277         CONTINUE
  785. *
  786. *           les cas non isotropes
  787. *
  788.             ELSE IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  789. *
  790.               DO 4379 IM=1,NMATT
  791.                 IF (IVAL(IM).NE.0) THEN
  792.                   MELVAL=IVAL(IM)
  793.                   IBMN=MIN(IB  ,VELCHE(/2))
  794.                   IGMN=MIN(IGAU,VELCHE(/1))
  795.                   VALMAT(IM)=VELCHE(IGMN,IBMN)
  796.                 ELSE
  797.                   VALMAT(IM)=0.D0
  798.                 ENDIF
  799.  4379         CONTINUE
  800. *
  801.               CALL PERMAO(WRK4,IFOUR,MATE,EREF,KERRE)
  802.               IF(KERRE.EQ.1) GO TO 99
  803.               IF(KERRE.EQ.2) THEN
  804.                 I367=IB
  805.                 GO TO 4079
  806.               ENDIF
  807. *
  808.               CALL BDBST(BGENE,DJAC,PMAT,LRE,NSTB,REL)
  809. *
  810. *           les cas non prevus
  811. *
  812.             ELSE
  813.               GO TO 99
  814.             ENDIF
  815. *
  816.  4079    CONTINUE
  817.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  818. *         SEGINI XMATRI
  819. *         IMATTT(IB)=XMATRI
  820. *
  821. *        Remplissage de XMATRI
  822. *
  823.          CALL REMPMT(REL,LRE,RE(1,1,ib))
  824. *         SEGDES XMATRI
  825.  3079 CONTINUE
  826. *
  827.       SEGSUP WRK1,WRK2,WRK3
  828.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  829.         SEGDES MINTE2
  830.         SEGSUP WRK4
  831.       ENDIF
  832. *
  833.       IF(I195.NE.0) THEN
  834.          INTERR(1)=I195
  835.          CALL ERREUR(195)
  836.          GOTO 9990
  837.       ELSE IF(I259.NE.0) THEN
  838.          INTERR(1)=I259
  839.          CALL ERREUR(259)
  840.          GOTO 9990
  841.       ELSE IF(I367.NE.0) THEN
  842.          INTERR(1)=I367
  843.          CALL ERREUR(367)
  844.          GOTO 9990
  845.       ENDIF
  846. *
  847.       SEGDES xMATRI
  848. *
  849.       GOTO 610
  850. *_______________________________________________________________________
  851. *
  852. *     JOINTS EN FORMULATION MILIEUX POREUX
  853. *_______________________________________________________________________
  854. *
  855.   80  CONTINUE
  856. *
  857. *     Pour ces elements NBBB : Nb de noeuds
  858. *                       NBNO : Nb de fonctions de forme
  859. *
  860.       NBNO=IPORE
  861.       NBBB=NBNN
  862.       NSTN=1
  863.       LPP=(NBNO-NBBB)*3/2
  864.       LRN=LPP
  865.       NSTB=2
  866.       IF(IFOUR.EQ.2) NSTB=3
  867.  
  868. *     PRINT *,' NBNO=', NBNO
  869. *     PRINT *,' NBBB=', NBBB
  870. *     PRINT *,' NSTN=', NSTN
  871. *     PRINT *,' LRN =', LRN
  872. *     PRINT *,' LRE =', LRE
  873. *     PRINT *,' NSTB =', NSTB
  874. *
  875.       SEGINI WRK1,WRK2,WRK3,WRK5
  876.       I195=0
  877.       I259=0
  878.       I367=0
  879.       DO 3080  IB=1,NBELEM
  880. *
  881. *        On cherche les coordonnees des noeuds de l'element IB
  882. *
  883.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  884. *
  885. *       calcul des axes locaux
  886. *
  887.       CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  888. *
  889.          CALL ZERO (REL,LRE,LRE)
  890. *
  891.       CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE)
  892. *
  893. *       boucle sur les points de Gauss
  894. *
  895.          ISDJC=0
  896.          DO 4080  IGAU=1,NBPGAU
  897. *
  898.       CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,XE,XEL,
  899.      .            SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,3)
  900.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  901.             IF(DJAC.EQ.0.) I259 =IB
  902.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  903. *
  904.             EREF=1.D0
  905.             MPTVAL=IVAMAT
  906. *
  907. *           le cas isotrope (le seul)
  908. *
  909.               MELVAL=IVAL(1)
  910.               IGMN=MIN(IGAU,VELCHE(/1))
  911.               IBMN=MIN(IB  ,VELCHE(/2))
  912.               XKT =VELCHE(IGMN,IBMN)
  913. *
  914.               MELVAL=IVAL(2)
  915.               IGMN=MIN(IGAU,VELCHE(/1))
  916.               IBMN=MIN(IB  ,VELCHE(/2))
  917.               XKNH =VELCHE(IGMN,IBMN)
  918. *
  919.               MELVAL=IVAL(3)
  920.               IGMN=MIN(IGAU,VELCHE(/1))
  921.               IBMN=MIN(IB  ,VELCHE(/2))
  922.               XKNB =VELCHE(IGMN,IBMN)
  923. *
  924.               MELVAL=IVAL(4)
  925.               IGMN=MIN(IGAU,VELCHE(/1))
  926.               IBMN=MIN(IB  ,VELCHE(/2))
  927.               XMU =VELCHE(IGMN,IBMN)
  928.               IF(XMU.EQ.0.D0) THEN
  929.                 I367=IB
  930.                 GO TO 4080
  931.               ENDIF
  932.               COMJAT=DJAC*EREF*EREF*XKT/XMU
  933.               COMJNH=DJAC*EREF*EREF*XKNH/XMU
  934.               COMJNB=DJAC*EREF*EREF*XKNB/XMU
  935.               DO 4281 I=1,LRN
  936.                  DO 4280 J=1,I
  937.                   REL(I,J)=REL(I,J)+COMJAT*BGENE(1,I)*BGENE(1,J)
  938.      .                     *XNTT(I)*XNTT(J)
  939.      .           +COMJNH*XGENE(1,I)*XGENE(1,J)*XNTH(I,J)
  940.      .           +COMJNB*XGENE(1,I)*XGENE(1,J)*XNTB(I,J)
  941.                    IF(IFOUR.EQ.2)THEN
  942.                   REL(I,J)=REL(I,J)+COMJAT*BGENE(2,I)*BGENE(2,J)
  943.      .                     *XNTT(I)*XNTT(J)
  944.                    ENDIF
  945.  4280        CONTINUE
  946.  4281        CONTINUE
  947. *
  948.  4080    CONTINUE
  949.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  950. *         SEGINI XMATRI
  951. *         IMATTT(IB)=XMATRI
  952. *
  953. *        Remplissage de XMATRI
  954. *
  955.          CALL REMPMT(REL,LRE,RE(1,1,ib))
  956. *         SEGDES XMATRI
  957.  3080 CONTINUE
  958. *
  959.       SEGSUP WRK1,WRK2,WRK3,WRK5
  960. *
  961.       IF(I195.NE.0) THEN
  962.          INTERR(1)=I195
  963.          CALL ERREUR(195)
  964.          GOTO 9990
  965.       ELSE IF(I259.NE.0) THEN
  966.          INTERR(1)=I259
  967.          CALL ERREUR(259)
  968.          GOTO 9990
  969.       ELSE IF(I367.NE.0) THEN
  970.          INTERR(1)=I367
  971.          CALL ERREUR(367)
  972.          GOTO 9990
  973.       ENDIF
  974. *
  975.       SEGDES xMATRI
  976. *
  977.       GOTO 610
  978. *_______________________________________________________________________
  979. *
  980. *     MILIEU POREUX  - SUITE
  981. *_______________________________________________________________________
  982. *
  983.  173  CONTINUE
  984. *
  985. *     Pour ces elements NBBB : Nb de noeuds
  986. *                       NBNO : Nb de fonctions de forme
  987. *
  988.       DIM3=1.D0
  989.       NBNO=IPORE
  990.       NBBB=NBNN
  991.       NSTN=IDECAP
  992.       LPP=NBNO-NBBB
  993.       LRN=IDECAP*LPP
  994.       NSTBE=2
  995.       IF(IFOUR.GT.0) NSTBE=3
  996.       NSTB=NSTBE*IDECAP
  997.  
  998. *     PRINT *,'NSTBE=',NSTBE
  999. *     PRINT *,'NSTB=',NSTB
  1000. *     PRINT *,'IDECAP=',IDECAP
  1001. *     PRINT *,'LRE   =',LRE
  1002.  
  1003. *
  1004. *     CAS NON ISOTROPES
  1005. *     NON PREVU
  1006. *
  1007.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  1008.         CALL ERREUR(251)
  1009.         GO TO 9990
  1010.       ENDIF
  1011. *
  1012.       NSTPK= NSTB
  1013.       SEGINI WRK1,WRK2,WRK3,WRK6
  1014.       I195=0
  1015.       I259=0
  1016.       I367=0
  1017.       DO 3173  IB=1,NBELEM
  1018.  
  1019. *     PRINT *,'ELEMENT ' , IB
  1020. *
  1021. *        On cherche les coordonnees des noeuds de l'element IB
  1022. *
  1023.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1024. *
  1025.          CALL ZERO (REL,LRE,LRE)
  1026. *
  1027. *       boucle sur les points de Gauss
  1028. *
  1029.          ISDJC=0
  1030.          DO 4173  IGAU=1,NBPGAU
  1031.  
  1032. *        PRINT *, ' POINT DE GAUSS ',IGAU
  1033.  
  1034.  
  1035. C
  1036. C   RECUPERATION DE L'EPAISSEUR
  1037. C
  1038.             IF (IFOUR.EQ.-2)THEN
  1039.              MPTVAL=IVACAR
  1040.              IF (IVACAR.NE.0) THEN
  1041.               MELVAL=IVAL(1)
  1042.               IF (MELVAL.NE.0) THEN
  1043.                 IGMN=MIN(IGAU,VELCHE(/1))
  1044.                 IBMN=MIN(IB,VELCHE(/2))
  1045.                 DIM3=VELCHE(IGMN,IBMN)
  1046.               ELSE
  1047.                 DIM3=1.D0
  1048.               ENDIF
  1049.              ENDIF
  1050.             ENDIF
  1051. *
  1052.             LHOO=NSTB
  1053.             CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  1054.      &         DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOO,2)
  1055.             IF(DJAC.LT.0.) ISDJC=ISDJC+1
  1056.             IF(DJAC.EQ.0.) I259 =IB
  1057.             DJAC=ABS(DJAC)*POIGAU(IGAU)
  1058.  
  1059. *           PRINT *,' MATRICE B LIGNE PAR LIGNE '
  1060. *           DO 1367 IPZ = 1,NSTB
  1061. *           PRINT *,' LIGNE ',IPZ
  1062. *           WRITE(6,1368) (BGENE(IPZ,JPZ), JPZ=1,LRE)
  1063. *1368        FORMAT(8(1X,1PE10.3)/)
  1064. *1367        CONTINUE
  1065.  
  1066. *
  1067.             EREF=1.D0
  1068.             MPTVAL=IVAMAT
  1069. *
  1070. *           le cas isotrope
  1071. *
  1072.             IF (MATE.EQ.1) THEN
  1073.  
  1074.               ICO=1
  1075.               DO 1731 ICD = 1,IDECAP
  1076.                 ICDA =(ICD -1) * NSTBE
  1077.                 DO 1732 JCD = 1,IDECAP
  1078.                   JCDA =(JCD -1) * NSTBE
  1079.                   MELVAL=IVAL(ICO)
  1080.                   IGMN=MIN(IGAU,VELCHE(/1))
  1081.                   IBMN=MIN(IB  ,VELCHE(/2))
  1082.                   DO 1733 KCD = 1,NSTBE
  1083.                     PKK(ICDA+KCD,JCDA+KCD) =VELCHE(IGMN,IBMN)
  1084. 1733            CONTINUE
  1085.                  ICO=ICO+1
  1086. 1732            CONTINUE
  1087. 1731          CONTINUE
  1088. *
  1089.  
  1090. *             PRINT *,' MATRICE PKK'
  1091. *              IF (IDECAP.EQ.2) THEN
  1092. *              WRITE (6,1342) ((PKK(I,J),J=1,NSTB),I=1,NSTB)
  1093. *1342          FORMAT(4(1X,1PE12.5)/)
  1094. *
  1095. *              ELSE IF (IDECAP.EQ.3) THEN
  1096. *              WRITE (6,1343) ((PKK(I,J),J=1,NSTB),I=1,NSTB)
  1097. *1343          FORMAT(6(1X,1PE12.5)/)
  1098. *              ENDIF
  1099.  
  1100.               COMJAC=DJAC*EREF*EREF
  1101.               CALL BDBSTS(BGENE,COMJAC,PKK,LRE,NSTB,REL)
  1102. *
  1103. *           les cas non prevus
  1104. *
  1105.             ELSE
  1106.               GO TO 99
  1107.             ENDIF
  1108. *
  1109.  4173    CONTINUE
  1110.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  1111. *         SEGINI XMATRI
  1112. *         IMATTT(IB)=XMATRI
  1113. *
  1114. *        Remplissage de XMATRI
  1115. *
  1116.          CALL REMPMS(REL,LRE,RE(1,1,ib))
  1117. *        SEGDES XMATRI
  1118.  3173 CONTINUE
  1119. *
  1120.       SEGSUP WRK1,WRK2,WRK3,WRK6
  1121.       IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  1122.         SEGDES MINTE2
  1123.         SEGSUP WRK4
  1124.       ENDIF
  1125. *
  1126.       IF(I195.NE.0) THEN
  1127.          INTERR(1)=I195
  1128.          CALL ERREUR(195)
  1129.          GOTO 9990
  1130.       ELSE IF(I259.NE.0) THEN
  1131.          INTERR(1)=I259
  1132.          CALL ERREUR(259)
  1133.          GOTO 9990
  1134.       ELSE IF(I367.NE.0) THEN
  1135.          INTERR(1)=I367
  1136.          CALL ERREUR(367)
  1137.          GOTO 9990
  1138.       ENDIF
  1139. *
  1140.       SEGDES xMATRI
  1141. *
  1142.       GOTO 610
  1143. *_______________________________________________________________________
  1144. *
  1145. *     JOINTS EN FORMULATION MILIEUX POREUX - SUITE
  1146. *_______________________________________________________________________
  1147. *
  1148.  185  CONTINUE
  1149. *
  1150. *     Pour ces elements NBBB : Nb de noeuds
  1151. *                       NBNO : Nb de fonctions de forme
  1152. *
  1153.       NBNO=IPORE
  1154.       NBBB=NBNN
  1155.       NSTN=IDECAP
  1156.       LPP=(NBNO-NBBB)*3/2
  1157.       LRN=IDECAP*LPP
  1158.       NSTBE=2
  1159.       IF(IFOUR.EQ.2) NSTBE=3
  1160.       NSTB=NSTBE*IDECAP
  1161.       NSTPKE=3
  1162.       NSTPK=NSTPKE*IDECAP
  1163.  
  1164. *     PRINT *,' NBNO=', NBNO
  1165. *     PRINT *,' NBBB=', NBBB
  1166. *     PRINT *,' NSTN=', NSTN
  1167. *     PRINT *,' LPP =', LPP
  1168. *     PRINT *,' LRN =', LRN
  1169. *     PRINT *,' LRE =', LRE
  1170. *     PRINT *,' NSTBE=', NSTBE
  1171. *     PRINT *,' NSTB =', NSTB
  1172. *     PRINT *,' NSTPKE =', NSTPKE
  1173. *     PRINT *,' NSTPK =', NSTPK
  1174. *
  1175.       SEGINI WRK1,WRK2,WRK3,WRK5,WRK6
  1176.       I195=0
  1177.       I259=0
  1178.       I367=0
  1179.       DO 3185  IB=1,NBELEM
  1180. *
  1181. *        On cherche les coordonnees des noeuds de l'element IB
  1182. *
  1183.          CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1184. *
  1185. *       calcul des axes locaux
  1186. *
  1187.       CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1188. *
  1189.          CALL ZERO (REL,LRE,LRE)
  1190. *
  1191.       CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE)
  1192. *
  1193. *       boucle sur les points de Gauss
  1194. *
  1195.       ISDJC=0
  1196.       DO 4185  IGAU=1,NBPGAU
  1197. *
  1198.       LHOO=NSTB
  1199.       CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOO,NSTN,XE,XEL,
  1200.      .     SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,3)
  1201.       IF(DJAC.LT.0.) ISDJC=ISDJC+1
  1202.       IF(DJAC.EQ.0.) I259 =IB
  1203.       DJAC=ABS(DJAC)*POIGAU(IGAU)
  1204. *
  1205.       EREF=1.D0
  1206.       MPTVAL=IVAMAT
  1207. *
  1208. *           le cas isotrope (le seul)
  1209. *
  1210.       IF(MATE.EQ.1) THEN
  1211.  
  1212.         FAC = DJAC*EREF*EREF
  1213.         IE=0
  1214.         DO 2184 IPR=1,IDECAP
  1215.           IPR1 = (IPR-1) * NSTPKE
  1216.           DO 2185 JPR=1,IDECAP
  1217.             JPR1 = (JPR-1) * NSTPKE
  1218.             DO 2186 I=1,NSTPKE
  1219.               II = I + IPR1
  1220.               JJ = I + JPR1
  1221.               IE=IE+1
  1222.               MELVAL=IVAL(IE)
  1223.               IGMN=MIN(IGAU,VELCHE(/1))
  1224.               IBMN=MIN(IB  ,VELCHE(/2))
  1225.               PKK(II,JJ)=VELCHE(IGMN,IBMN)*FAC
  1226. 2186        CONTINUE
  1227. 2185    CONTINUE
  1228. 2184    CONTINUE
  1229.  
  1230. *
  1231.         DO 8984 IPR=1,IDECAP
  1232.           IPR1 = (IPR-1)*NSTPKE
  1233.           IPR2 = 2*IPR
  1234.           IPPDEC=(IPR-1)*LPP
  1235.           IRRDEC=(IPR-1)*NBBB
  1236.           DO 8985 JPR=1,IDECAP
  1237.             JPR1 = (JPR-1)*NSTPKE
  1238.             JPR2 = 2*JPR
  1239.             JPPDEC=(JPR-1)*LPP
  1240.             JRRDEC=(JPR-1)*NBBB
  1241. *
  1242.             COMJAT=PKK(IPR1+1,JPR1+1)
  1243.             COMJNH=PKK(IPR1+2,JPR1+2)
  1244.             COMJNB=PKK(IPR1+3,JPR1+3)
  1245. *
  1246.             IF(IFOUR.LE.0) THEN
  1247.               DO 4284 I=1,LPP
  1248.                 II =I+IPPDEC
  1249.                 IR =I+IRRDEC
  1250.                 DO 4285 J=1,LPP
  1251.                   JJ =J+JPPDEC
  1252.                   JR =J+JRRDEC
  1253.                   REL(IR,JR)=REL(IR,JR)
  1254.      .       +COMJAT*BGENE(IPR,II)*BGENE(JPR,JJ)
  1255.      .                   *XNTT(I)*XNTT(J)
  1256.      .       +COMJNH*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTH(I,J)
  1257.      .       +COMJNB*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTB(I,J)
  1258.  4285     CONTINUE
  1259.  4284     CONTINUE
  1260. *
  1261.             ELSE
  1262.               DO 4384 I=1,LPP
  1263.                 II =I+IPPDEC
  1264.                 DO 4385 J=1,LPP
  1265.                   JJ =J+JPPDEC
  1266.                   REL(IR,JR)=REL(IR,JR)
  1267.      .            +COMJAT*XNTT(I)*XNTT(J)*
  1268.      .             (BGENE(IPR2-1,II)*BGENE(JPR2-1,JJ)
  1269.      .              + BGENE(IPR2,II)*BGENE(JPR2,JJ))
  1270.      .           +COMJNH*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTH(I,J)
  1271.      .           +COMJNB*XGENE(IPR,II)*XGENE(JPR,JJ)*XNTB(I,J)
  1272.  4385         CONTINUE
  1273.  4384         CONTINUE
  1274.           ENDIF
  1275.  
  1276.  8985   CONTINUE
  1277.  8984   CONTINUE
  1278. *
  1279.       ELSE
  1280.         GO TO 9990
  1281.       ENDIF
  1282.  
  1283. *
  1284.  4185    CONTINUE
  1285.          IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB
  1286. *
  1287. *        Remplissage de XMATRI
  1288. *
  1289.          CALL REMPMS(REL,LRE,RE(1,1,ib))
  1290.  
  1291.  3185 CONTINUE
  1292. *
  1293.       SEGSUP WRK1,WRK2,WRK3,WRK5,WRK6
  1294. *
  1295.       IF(I195.NE.0) THEN
  1296.          INTERR(1)=I195
  1297.          CALL ERREUR(195)
  1298.          GOTO 9990
  1299.       ELSE IF(I259.NE.0) THEN
  1300.          INTERR(1)=I259
  1301.          CALL ERREUR(259)
  1302.          GOTO 9990
  1303.       ELSE IF(I367.NE.0) THEN
  1304.          INTERR(1)=I367
  1305.          CALL ERREUR(367)
  1306.          GOTO 9990
  1307.       ENDIF
  1308. *
  1309.       GOTO 610
  1310. *
  1311. *     Desactivation des segment propres a la geometrie ISOUS
  1312. *
  1313.   610 CONTINUE
  1314. *
  1315.       SEGDES MELEME
  1316.       SEGDES MINTE
  1317. *
  1318.       IF (ISUP.EQ.1) THEN
  1319.         CALL DTMVAL(IVAMAT,3)
  1320.       ELSE
  1321.         CALL DTMVAL(IVAMAT,1)
  1322.       ENDIF
  1323. *
  1324.       NOMID=MOMATR
  1325.       SEGSUP NOMID
  1326.       IF (lsupdp) THEN
  1327.         NOMID=MODEPL
  1328.         SEGSUP,NOMID
  1329.       ENDIF
  1330.       IF (lsupfo) THEN
  1331.         NOMID=MOFORC
  1332.         SEGSUP,NOMID
  1333.       ENDIF
  1334. *
  1335.       COERIG(ISORI)  = 1.D0
  1336.       IRIGEL(1,ISORI)=IPMAIL
  1337.       IRIGEL(2,ISORI)=0
  1338.       IRIGEL(3,ISORI)=IPDES
  1339.       IRIGEL(4,ISORI)=xMATRI
  1340.       IRIGEL(5,ISORI)=NHRM
  1341.       IRIGEL(6,ISORI)=0
  1342.       IRIGEL(7,ISORI)=LASYM
  1343.       segact xmatri*mod 
  1344.       xmatri.symre=lasym
  1345.       SEGDES xMATRI
  1346.       IRIGEL(8,ISORI)=0
  1347.  
  1348. *  Fin de la boucle de PARTITIONNEMENT du segment XMATRI
  1349.  5000 CONTINUE
  1350.  
  1351.       SEGDES IMODEL
  1352.  
  1353.   500 CONTINUE
  1354.       IF(ISORI.NE.NRIGEL) GO TO 9999
  1355.       SEGDES MRIGID
  1356. C      SEGSUP MMODEL
  1357.       IRET = 1
  1358.       IPRIGI = MRIGID
  1359.       RETURN
  1360. *
  1361. *     Erreur dans une sous zone desactivation et retour
  1362. *
  1363.  9990 CONTINUE
  1364.       IF (ISUP.EQ.1) THEN
  1365.         CALL DTMVAL(IVAMAT,3)
  1366.       ELSE
  1367.         CALL DTMVAL(IVAMAT,1)
  1368.       ENDIF
  1369. *
  1370.  9992 CONTINUE
  1371.  
  1372.       SEGSUP DESCR
  1373.       SEGSUP xMATRI
  1374.       SEGDES MELEME
  1375.       SEGDES MINTE
  1376. *
  1377.       NOMID=MOMATR
  1378.       SEGSUP NOMID
  1379.  9999 CONTINUE
  1380.       SEGSUP MRIGID
  1381.       IRET = 0
  1382.       IPRIGI = 0
  1383.  
  1384.  9991 CONTINUE
  1385.  
  1386.       MMODEL = IPMODL
  1387.       DO isous = 1, NSOUS
  1388.         IMODEL = KMODEL(isous)
  1389.         SEGDES,IMODEL
  1390.       ENDDO
  1391. C      SEGSUP,MMODEL
  1392.  
  1393.       RETURN
  1394.       END
  1395.  
  1396.  
  1397.  
  1398.  
  1399.  
  1400.  
  1401.  

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