Télécharger la plaquette Cast3M


fr en

Télécharger ksigmp.eso

Retour à la liste

Numérotation des lignes : 

ksigmp

  1. C KSIGMP    SOURCE    JK148537  26/06/23    21:15:04     12579          
  2.  
  3.       SUBROUTINE KSIGMP(IPMODL,IPCHE1,IPCHE2,IFLAM, IPRIGG)
  4.  
  5. c_______________________________________________________________________
  6. c
  7. c
  8. c     construction de la matrice de raideur geometrique a partir d'un
  9. c     mchaml de contraintes
  10. c
  11. c     entr{es:
  12. c     ________
  13. c
  14. c       ipmodl  pointeur sur un mmodel
  15. c       ipche1  pointeur sur un mchaml de contraintes
  16. c       ipche2  pointeur sur un mchaml de caracteristiques
  17. c       iflam   flag de flambage
  18. c
  19. c     sorties:
  20. c     ________
  21. c
  22. c       iprigg  pointeur sur un objet rigidite
  23. c               = 0 en cas d'erreur
  24. c_______________________________________________________________________
  25. c
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8(A-H,O-Z)
  28.  
  29. -INC PPARAM
  30. -INC CCOPTIO
  31. -INC CCHAMP
  32. -INC CCREEL
  33.  
  34. -INC SMCHAML
  35. -INC SMCOORD
  36. -INC SMELEME
  37. -INC SMINTE
  38. -INC SMMODEL
  39. -INC SMRIGID
  40.  
  41. -INC TMPTVAL
  42.  
  43.       INTEGER oooval
  44.  
  45.       SEGMENT NOTYPE
  46.         CHARACTER*16 TYPE(NBTYPE)
  47.       ENDSEGMENT
  48.  
  49.       SEGMENT MWRK1
  50.         REAL*8 REL(LRE,LRE) ,XE(3,NBBB) ,XSTRS(NSTRS)
  51.       ENDSEGMENT
  52. C
  53.       SEGMENT MWRK2
  54.         REAL*8 SHPWRK(6,NBNO) ,BGENE(NSTRS,LRE)
  55.       ENDSEGMENT
  56. C
  57.       SEGMENT MWRK3
  58.         REAL*8 WORK(LW)
  59.       ENDSEGMENT
  60. C
  61.       SEGMENT MWRK4
  62.         REAL*8 BPSS(3,3) ,XEL(3,NBBB)
  63.       ENDSEGMENT
  64. C
  65.       SEGMENT MWRK5
  66.         REAL*8 GEOM(20), tabw(6,9), tabrot(4,9), XX(3), YY(3)
  67.       ENDSEGMENT
  68. C
  69. C   segment pour shb8
  70.       SEGMENT MWRK7
  71.         REAL*8 PROPEL(1),out(1),d(1), work1(30)
  72.       ENDSEGMENT
  73. C
  74.       character*6 msorse
  75.       CHARACTER*8 CMATE
  76.       CHARACTER*(NCONCH) CONM
  77.       PARAMETER ( NINF=3 )
  78.       INTEGER INFOS(NINF)
  79.       LOGICAL lsupfo,lsupde,lsupco,BDPGE
  80.       INTEGER ISUP1,ISUP2
  81.  
  82.       ISUP1 = 0
  83.       ISUP2 = 0
  84.  
  85.       IPRIGG = 0
  86.  
  87.       IDIMP1 = IDIM+1
  88. C
  89. C      verification du lieu support du mchaml de contraintes
  90. C
  91.       CALL QUESUP(IPMODL,IPCHE1,3,0,ISUP1,IRET1C)
  92.       IF (ISUP1.GT.1) RETURN
  93. C
  94. C      verification du lieu support du mchaml de caracteristiques
  95. C
  96.       IF (IPCHE2.NE.0) THEN
  97.         CALL QUESUP(IPMODL,IPCHE2,3,0,ISUP2,iret2c)
  98.         IF (ISUP2.GT.1) RETURN
  99.       ENDIF
  100. c
  101. c_______________________________________________________________________
  102. c
  103. c     initialisation du chapeau de l objet rigidite
  104. c_______________________________________________________________________
  105. c
  106.       NRIGEL = 0
  107.       SEGINI,MRIGID
  108.       IFORIG = IFOUR
  109.       ICHOLE = 0
  110.       IMGEO1 = 0
  111.       IMGEO2 = 0
  112.       ISUPEQ = 0
  113.       IF (IFLAM.NE.0) THEN
  114.         MTYMAT = 'MASSE   '
  115.       ELSE
  116.         MTYMAT = 'RIGIDITE'
  117.       ENDIF
  118. c
  119. c_______________________________________________________________________
  120. c
  121. c     activation du modele
  122. c_______________________________________________________________________
  123. c
  124.       MMODEL = IPMODL
  125.       SEGACT,MMODEL
  126.       NSOUS = KMODEL(/1)
  127. c
  128. c      boucle sur les modeles elementaires
  129. c
  130.       DO 500 ISOUS = 1,NSOUS
  131. c
  132. c     traitement du modele
  133. c
  134.         IMODEL = KMODEL(ISOUS)
  135.         IF (FORMOD(1).EQ.'CHARGEMENT') GOTO 500
  136.         SEGACT,IMODEL
  137. C
  138. C     INITIALISATIONS
  139. C
  140.         IPMINT = 0
  141.         IPMIN1 = 0
  142.  
  143.         MOSTRS = 0
  144.         MOCARA = 0
  145.         MOTYPS = 0
  146.         MOTYPC = 0
  147.  
  148.         MODEPL = 0
  149.         MOFORC = 0
  150.         lsupde = .false.
  151.         lsupfo = .false.
  152.         lsupco = .false.
  153.  
  154.         IDESCR = 0
  155.  
  156. C- Recuperation d'informations sur le maillage elementaire
  157.         IPT1 = IMAMOD
  158.         SEGACT,IPT1
  159.         NBNOE1 = IPT1.NUM(/1)
  160.         NBELE1 = IPT1.NUM(/2)
  161.  
  162. C- Quelques informations sur le modele
  163.         IIPDPG = imodel.IPDPGE
  164.         IIPDPG = IPTPOI(IIPDPG)
  165.  
  166.         CONM  = CONMOD
  167.         CMATE = CMATEE
  168. C       MATE  = IMATEE
  169. c*      INAT  = INATUU
  170.  
  171.         IRTD = 1
  172.         CALL IDENT(IPT1,CONM,IPCHE1,IPCHE2, INFOS,IRTD)
  173.         IF (IRTD.EQ.0) GOTO 599
  174. C
  175. C- Recuperation d'informations sur l'element fini
  176.         MELE  = NEFMOD
  177. c  pour l'el. timo on utilise l'el. barr
  178. c         IF (MELE .EQ. 84) MELE = 46
  179. c  bp: comme il n y a plus elquoi, ce n'est pas ici que ca intervient...
  180.  
  181. c   coque integree ou pas ?
  182.         NPINT = INFMOD(1)
  183.         IF (NPINT.NE.0) THEN
  184.           CALL ERREUR(615)
  185.           GOTO 599
  186.         ENDIF
  187.  
  188. C       LHOOK = INFELE(10)
  189. c*      LHOO2 = LHOOK*LHOOK
  190.         NSTRS = INFELE(16)
  191.         MFR   = INFELE(13)
  192.         LW    = INFELE(7)
  193. C        NDDL  = INFELE(15)
  194.         LRE   = INFELE(9)
  195. C        IPORE = INFELE(8)
  196.         NHRM  = NIFOUR
  197.  
  198.         IPPORE = 0
  199.         IF (MFR.EQ.33) IPPORE = NBNOE1
  200.  
  201. c_______________________________________________________________________
  202. C   segments d'integration     *
  203. c_______________________________________________________________________
  204. C  minte : 1er  segment d'integration, il existe pour tous les e.f.
  205. C  minte1: 2eme segment d'integration, uniquement pour certains e.f.
  206. C          en particulier pour coq6 et coq8
  207. C  nbpg:nb de points de gauss = nbpgau du segment minte
  208. C  iele:no d'element geometrique associe a l'e.f. mele
  209. C  nbff:nb de fonctions de forme = nbno du segment minte
  210.         NBPGAU = INFELE(  6)
  211. C        IELE   = INFELE( 14)
  212. c*        ICARA  = INFELE(  5)
  213.         IPMINT = INFMOD(5)
  214. c*        IPMINT  = INFELE(11)
  215.         IPMIN1 = INFMOD(3)
  216.         MINTE = IPMINT
  217.         IF (IPMINT.NE.0) SEGACT,MINTE
  218.  
  219. c_______________________________________________________________________
  220. c
  221. C  initialisation du segment descr, segment descripteur des *
  222. C  des inconnues relatives a la matrice de rigidite         *
  223. c_______________________________________________________________________
  224.         if (lnomid(1).ne.0) then
  225.           MODEPL = lnomid(1)
  226.         else
  227.           lsupde = .true.
  228.           CALL IDPRIM(IMODEL,MFR,MODEPL,NDEPL,NDUM)
  229.         endif
  230.         nomid = MODEPL
  231.         segact,nomid
  232.         ndepl = lesobl(/2)
  233. c*        ndum  = lesfac(/2)
  234.  
  235.         if (lnomid(2).ne.0) then
  236.           moforc = lnomid(2)
  237.         else
  238.           lsupfo=.true.
  239.           CALL IDDUAL(IMODEL,MFR,MOFORC,NFORC,NDUM)
  240.         endif
  241.         nomid = MOFORC
  242.         segact,nomid
  243.         nforc = lesobl(/2)
  244. c*      ndum  = lesfac(/2)
  245.  
  246.         IF (NDEPL.EQ.0.OR.NFORC.EQ.0.OR.NDEPL.NE.NFORC) CALL ERREUR(5)
  247.  
  248. C CB215821 : On copie RIGI2 pour les deformations planes generalisees
  249.         IF (IIPDPG.GT.0) THEN
  250.           IF (IFOUR.EQ.-3) THEN
  251.             BDPGE=.TRUE.
  252.             IREF=(IIPDPG-1)*(IDIM+1)
  253. C            XDPGE=XCOOR(IREF+1)
  254. C            YDPGE=XCOOR(IREF+2)
  255.           ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  256.      &             IFOUR.EQ.10 .OR. IFOUR.EQ.11 .OR. IFOUR.EQ.14) THEN
  257.             BDPGE=.TRUE.
  258. C            XDPGE=XZero
  259. C            YDPGE=XZero
  260.           ELSE
  261.             CALL ERREUR(21)
  262.             RETURN
  263.           ENDIF
  264.         ELSE
  265.           BDPGE=.FALSE.
  266. C          XDPGE=XZero
  267. C          YDPGE=XZero
  268.         ENDIF
  269.  
  270.         NCOMP=NDEPL
  271.         NBNNS=NBNOE1
  272.         IF (MFR.EQ.33)                NCOMP = NDEPL-1
  273.         IF (MFR.EQ.19 .OR. MFR.EQ.21) NBNNS = NBNOE1/2
  274. C
  275. C   REMPLISSAGE DU SEGMENT DESCRIPTEUR
  276.         IF(BDPGE)THEN
  277. C         CB215821 : En 'DPGE' on traite sans le Pt support
  278.           NCOMP = NCOMP - 3
  279.           LRE   = LRE   - 3
  280.         ENDIF
  281.         NLIGRP = LRE
  282.         NLIGRD = LRE
  283.         SEGINI,DESCR
  284.  
  285.         IDDL = 1
  286.         DO 1004 INOEUD=1,NBNNS
  287.           DO 1005 ICOMP=1,NCOMP
  288.             NOMID=MODEPL
  289.             LISINC(IDDL)=LESOBL(ICOMP)
  290.             NOMID=MOFORC
  291.             LISDUA(IDDL)=LESOBL(ICOMP)
  292.             NOELEP(IDDL)=INOEUD
  293.             NOELED(IDDL)=INOEUD
  294.             IDDL=IDDL+1
  295.  1005     CONTINUE
  296.  1004   CONTINUE
  297.  
  298. C           cas des milieux poreux
  299. C
  300. C         if (mfr.eq.33) then
  301. C           ipos = nspos(iele)
  302. C           do 1104 inoeud=1,nbsom(iele)
  303. C             nomid=modepl
  304. C             lisinc(iddl)=lesobl(ndepl)
  305. C             nomid=moforc
  306. C             lisdua(iddl)=lesobl(ndepl)
  307. C             i = ibsom(ipos+inoeud-1)
  308. C             noelep(iddl)=i
  309. C             noeled(iddl)=i
  310. C             iddl=iddl+1
  311. C 1104      continue
  312. C         endif
  313.  
  314. C           cas des element raccord
  315. C
  316.         IF (MFR.EQ.19.OR.MFR.EQ.21) THEN
  317.           CALL IDPRIM(IMODEL,MFR+1000,MODPL1,NDEPL,NDUM)
  318.           CALL IDDUAL(IMODEL,MFR+1000,MOFRC1,NFORC,NDUM)
  319.           DO 1106 INOEUD=NBNNS+1,NBNOE1
  320.             DO 1107 ICOMP=1,NDEPL
  321.               NOMID=MODPL1
  322.               LISINC(IDDL)=LESOBL(ICOMP)
  323.               NOMID=MOFRC1
  324.               LISDUA(IDDL)=LESOBL(ICOMP)
  325.               NOELEP(IDDL)=INOEUD
  326.               NOELED(IDDL)=INOEUD
  327.               IDDL=IDDL+1
  328. 1107        continue
  329. 1106      continue
  330.           NOMID=MODPL1
  331.           SEGSUP,NOMID
  332.           NOMID=MOFRC1
  333.           SEGSUP,NOMID
  334.         ENDIF
  335.  
  336.         SEGDES,DESCR
  337.         IDESCR = DESCR
  338. c_______________________________________________________________________
  339. c
  340. C  composantes de contraintes necessaires                      *
  341. c_______________________________________________________________________
  342.         if (lnomid(4).ne.0) then
  343.           MOSTRS = lnomid(4)
  344.         else
  345.           lsupco=.true.
  346.           CALL IDCONT(IMODEL,IFOUR,MOSTRS,NSTR,NFAC)
  347.         endif
  348.         nomid = MOSTRS
  349.         segact,nomid
  350.         nstr=lesobl(/2)
  351. c*        nfac=lesfac(/2)
  352. c*        write(6,*) 'mostrts',mostrs,nstr,nfac
  353.         nbtype = 1
  354.         SEGINI,notype
  355.         TYPE(1)='REAL*8'
  356.         MOTYPS = notype
  357.  
  358.         ifai = 1
  359.         if (mele.eq.260.and.IRET1C.eq.5) ifai = 0
  360.         ISUP1L = 0
  361.         IF (ISUP1.EQ.1.AND.ifai.eq.1) ISUP1L = 1
  362.  
  363. c____________________________________________________________________
  364. c
  365. C  traitement des champs de caracteristiques                   *
  366. c____________________________________________________________________
  367.         NBROBL = 0
  368.         NBRFAC = 0
  369.         IVECT  = 0
  370.         notype = 0
  371.         nomid  = 0
  372. C
  373. C       v1x v1y dans le cas de la coque dst orthotrope
  374. C
  375.         IF (MFR.EQ.9) THEN
  376.           IF (MELE.EQ.93.AND.CMATE.NE.'ISOTROPE')THEN
  377.             NBROBL=2
  378.             SEGINI NOMID
  379.             LESOBL(1)='V1X '
  380.             LESOBL(2)='V1Y '
  381. C
  382.             NBTYPE=1
  383.             SEGINI NOTYPE
  384.             TYPE(1)='REAL*8'
  385.           ENDIF
  386. C
  387. C  epaisseur dans le cas massif et coq2 en contraintes planes
  388. C
  389.         ELSE IF ( (MFR.EQ.1.OR.MFR.EQ.3.OR.MFR.EQ.31) .AND.
  390.      +            IFOUR.EQ.-2 .AND. IPCHE2.NE.0) THEN
  391.           NBRFAC=1
  392.           SEGINI NOMID
  393.           LESFAC(1)='DIM3'
  394. C
  395.           NBTYPE=1
  396.           SEGINI NOTYPE
  397.           TYPE(1)='REAL*8'
  398. C
  399. C  epaisseur et excentrement dans le cas des coques epaisses
  400. C
  401.         ELSE IF (MFR.EQ.5 .OR. (MFR.EQ.3.AND.IFOUR.NE.-2)) THEN
  402.           NBROBL=1
  403.           NBRFAC=1
  404.           SEGINI NOMID
  405.           LESOBL(1)='EPAI'
  406.           LESFAC(1)='EXCE'
  407.  
  408.           NBTYPE=1
  409.           SEGINI NOTYPE
  410.           TYPE(1)='REAL*8'
  411. C
  412. C  caracteristiques pour les poutres
  413. C
  414.         ELSE IF (MFR.EQ.7 ) THEN
  415.           IF (IFOUR.EQ.-2.OR.IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN
  416.             NBROBL=2
  417.             NBRFAC=1
  418.             SEGINI NOMID
  419.             LESOBL(1)='SECT'
  420.             LESOBL(2)='INRZ'
  421.             LESFAC(1)='SECY'
  422. C
  423.             NBTYPE=1
  424.             SEGINI NOTYPE
  425.             TYPE(1)='REAL*8'
  426.           ELSE
  427.             NBROBL=4
  428.             NBRFAC=5
  429.             IVECT =1
  430.             SEGINI NOMID
  431.             LESOBL(1)='TORS'
  432.             LESOBL(2)='INRY'
  433.             LESOBL(3)='INRZ'
  434.             LESOBL(4)='SECT'
  435.             LESFAC(1)='SECY'
  436.             LESFAC(2)='SECZ'
  437.             LESFAC(3)='VX  '
  438.             LESFAC(4)='VY  '
  439.             LESFAC(5)='VZ  '
  440. C
  441.             NBTYPE=9
  442.             SEGINI NOTYPE
  443.             TYPE(1)='REAL*8'
  444.             TYPE(2)='REAL*8'
  445.             TYPE(3)='REAL*8'
  446.             TYPE(4)='REAL*8'
  447.             TYPE(5)='REAL*8'
  448.             TYPE(6)='REAL*8'
  449.             TYPE(7)='REAL*8'
  450.             TYPE(8)='REAL*8'
  451.             TYPE(9)='REAL*8'
  452.           ENDIF
  453. C
  454. C  caracteristiques pour les tuyaux
  455. C
  456.         ELSE IF (MFR.EQ.13) THEN
  457.           NBROBL = 2
  458.           NBRFAC = 5
  459.           IVECT  = 1
  460.           SEGINI NOMID
  461.           LESOBL(1)='EPAI'
  462.           LESOBL(2)='RAYO'
  463.           LESFAC(1)='RACO'
  464.           LESFAC(2)='CISA'
  465.           LESFAC(3)='VX  '
  466.           LESFAC(4)='VY  '
  467.           LESFAC(5)='VZ  '
  468. C
  469.           NBTYPE = 7
  470.           SEGINI NOTYPE
  471.           TYPE(1)='REAL*8'
  472.           TYPE(2)='REAL*8'
  473.           TYPE(3)='REAL*8'
  474.           TYPE(4)='REAL*8'
  475.           TYPE(5)='REAL*8'
  476.           TYPE(6)='REAL*8'
  477.           TYPE(7)='REAL*8'
  478.         ENDIF
  479. C
  480.         MOCARA = NOMID
  481.         MOTYPC = NOTYPE
  482.         NCARA = NBROBL
  483.         NCARF = NBRFAC
  484.         NCARR = NCARA+NCARF
  485. **      write(6,*) 'mfr ncarr en 498',mfr,ncarr,ncara,ncarf
  486.  
  487.         IF (MOCARA.NE.0 .AND. IPCHE2.EQ.0) THEN
  488.           MOTERR(1:8)  = 'CARACTER'
  489.           MOTERR(9:12) = NOMTP(MELE)
  490.           MOTERR(13:20)= 'KSIGMA'
  491.           CALL ERREUR(145)
  492.           GOTO 598
  493.         ENDIF
  494.  
  495.         ifai = 1
  496.         IF (mele.EQ.260) ifai = 0
  497.         ISUP2L = 0
  498.         IF (ISUP2.EQ.1.AND.ifai.eq.1) ISUP2L = 1
  499.  
  500. C- Partionnement si necessaire de la matrice de capacite
  501. C- determinant ainsi le nombre d'objets elementaires de MRIGID
  502.         LTRK = oooval(1,4)
  503.         IF (LTRK.EQ.0) LTRK = oooval(1,1)
  504.       LTRK=MAX(LTRK,2**24)
  505. C  Ajout a la taille en mots de la matrice des infos du segment
  506.         LSEG = LRE*LRE*NBELE1 + 16
  507.         NBLPRT = (LSEG-1)/LTRK + 1
  508.         NBLMAX = (NBELE1-1)/NBLPRT + 1
  509.         NBLPRT = (NBELE1-1)/NBLMAX + 1
  510. C         write(ioimp,*) ' capa1 : nblprt nblmax = ',nblprt,nblmax,nbele1
  511.  
  512. C Ajout de la matrice a la matrice globale
  513. C ========================================
  514. C        NRIGE0 =  IRIGEL(/2)
  515. C         NRIGEL = NRIGE0 + NBLPRT
  516. C         SEGADJ,MRIGID
  517.  
  518.         descr  = IDESCR
  519.         meleme = IPT1
  520.         NBNN   = NBNOE1
  521.         nbelem = NBELE1
  522.         nbsous = 0
  523.         nbref  = 0
  524. C
  525. C ***********************************************************************
  526. C                           P H A S E   2
  527. C
  528. C  Boucle sur les PARTITIONS elementaires de la matrice
  529. C
  530. C ***********************************************************************
  531.         DO irige = 1, NBLPRT
  532.  
  533.           IF (NBLPRT.GT.1) THEN
  534. C- Partitionnement du maillage support de la matrice elementaire
  535. C- (IPT1 peut etre desactive suite a l'appel a KOMCHA !)
  536.             SEGACT,IPT1
  537.             ielem  = (irige-1)*NBLMAX
  538.             nbelem = MIN(NBLMAX,NBELE1-ielem)
  539. C           write(ioimp,*) ' creation segment ',nbnn,nbelem
  540.             SEGINI,meleme
  541.             itypel = IPT1.itypel
  542.             DO ielt = 1, nbelem
  543.               jelt = ielt + ielem
  544.               DO inoe = 1, NBNN
  545.                 num(inoe,ielt) = IPT1.NUM(inoe,jelt)
  546.               ENDDO
  547.               icolor(ielt) = IPT1.ICOLOR(jelt)
  548.             ENDDO
  549. C- Recopie du descripteur
  550.             des1 = IDESCR
  551.             SEGINI,descr=des1
  552.             SEGDES,descr
  553.           ENDIF
  554.  
  555.           ipmail = meleme
  556.           ipdesc = descr
  557.  
  558. C- Initialisation de la matrice de rigidite elementaire (xmatri)
  559.           NELRIG = nbelem
  560.           RIGREL=0
  561.           SEGINI,xmatri
  562.           ipmatr = xmatri
  563.  
  564. C- Recuperation des valeurs des contraintes et proprietes geometriques
  565.           IVASTR = 0
  566.           IVACAR = 0
  567.           IVECTL = IVECT
  568.           NCARR1 = NCARR
  569. C
  570.           CALL KOMCHA(IPCHE1,ipmail,CONM,MOSTRS,MOTYPS,1,INFOS,3,IVASTR)
  571.           IF (IERR.NE.0) GOTO 597
  572.           IF (ISUP1L.EQ.1) THEN
  573.             CALL VALCHE(IVASTR,NSTR,IPMINT,IPPORE,MOSTRS,MELE)
  574.             IF (IERR.NE.0) THEN
  575.               ISUP1L = 0
  576.               GOTO 597
  577.             ENDIF
  578.           ENDIF
  579.  
  580.           IF (MOCARA.NE.0)  THEN
  581.             CALL KOMCHA(IPCHE2,ipmail,CONM,MOCARA,MOTYPC,1,
  582.      &                  INFOS,3,IVACAR)
  583.             IF (IERR.NE.0) GOTO 597
  584.             IF (ISUP2L.EQ.1) THEN
  585.               CALL VALCHE(IVACAR,NCARR,IPMINT,IPPORE,MOCARA,MELE)
  586.               IF (IERR.NE.0) THEN
  587.                 ISUP2L = 0
  588.                 GOTO 597
  589.               ENDIF
  590.             ENDIF
  591.             IF (IVECT.EQ.1) THEN
  592. ***           MPTVAL = IVACAR
  593. *****         NCARR1 = NCARR - 3
  594.  
  595. ***           IF (IVAL(NCARR1).EQ.0) IVECTL = 2
  596.             ENDIF
  597.           ENDIF
  598. **        write(6,*) ' dans ksigmp ivect ivectl ',ivect,ivectl
  599.  
  600. c_______________________________________________________________________
  601. c
  602. c     numero des etiquettes      :
  603. c     etiquettes de 1 a 98 pour traitement specifique a l element
  604. c     dans la zone specifique a chaque element commencant par :
  605. c     5  continue
  606. c     element 5   etiquettes 1005 2005 3005 4005   ...
  607. c     44 continue
  608. c     element 44  etiquettes 1044 2044 3044 4044   ...
  609. c_______________________________________________________________________
  610.           IF (MELE.LE.100) THEN
  611.       GOTO (99,99,99, 4,99, 4,99, 4,99, 4,99,99,99, 4, 4, 4, 4,99,99,99,
  612.      1      99,99, 4, 4, 4, 4,27,28,29,99,99,99,99,99,99,99,99,99,99,99,
  613.      2      41,29,43,44,99,46,99,99,49,99,51,99,99,99,99,41,99,99,99,99,
  614.      3      99,99,99,99,99,99,99,99, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,99,99,
  615.      4      99,99,99,29,99,99,99,99,99,99,99,99,28,99,46,99,99,99,99,99
  616.      5     ),MELE
  617.           ELSE IF (MELE.LE.200) THEN
  618.       GOTO (99,99,99,99,99,99,99,99,99,99, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
  619.      1       4, 4,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  620.      2      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  621.      3      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  622.      4      99,99, 4, 4,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99
  623.      5     ),MELE-100
  624.           ELSE IF (MELE.LE.300) THEN
  625.       GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  626.      1      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  627.      2      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  628.      &      260,
  629.      3      99,99,99,99,99,99,99,99,99,99,99,99,99, 4, 4,99,99,99,99,99,
  630.      4      99,99, 4, 4,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99
  631.      5     ),MELE-200
  632.           ENDIF
  633. c
  634.  99       CONTINUE
  635. C           MOTERR(1:4)  = NOMTP(MELE)
  636. C           MOTERR(5:12) = 'KSIGMP  '
  637. C          CALL ERREUR(86)
  638.           GOTO 510
  639. c
  640. c_______________________________________________________________________
  641. c
  642. c     secteur de calcul pour les elements massifs
  643. c_______________________________________________________________________
  644.  4        CONTINUE
  645.           NBNO = NBNN
  646.           NBBB = NBNN
  647.           SEGINI,MWRK1,MWRK2
  648. c   recuperation de l'epaisseur
  649.           DIM3 = 1.D0
  650.           MEPDI3 = 0
  651. c*         IF (IFOUR.EQ.-2.AND.IPCHE2.NE.0) THEN
  652.           IF (IVACAR.NE.0) THEN
  653.             MPTVAL = IVACAR
  654.             MEPDI3 = IVAL(1)
  655.           ENDIF
  656.  
  657.           DO 3004 IB=1,NBELEM
  658. c
  659. c     on cherche  les coordonnees des noeuds de l element ib
  660. c
  661.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  662.             CALL ZERO(REL,LRE,LRE)
  663. c
  664. c     boucle sur les points de gauss
  665. c
  666.             ISDJC = 0
  667.             DO 4004 IGAU=1,NBPGAU
  668. c
  669. c   recuperation de l'epaisseur
  670.               IF (MEPDI3.NE.0) THEN
  671.                 MELVAL = MEPDI3
  672.                 IGMN=MIN(IGAU,VELCHE(/1))
  673.                 IBMN=MIN(  IB,VELCHE(/2))
  674.                 DIM3=VELCHE(IGMN,IBMN)
  675.               ENDIF
  676. c
  677.               DO 100 IA=1,NBNN
  678.                 DO 101 IO=1,IDIMP1
  679.                   SHPWRK(IO,IA)=SHPTOT(IO,IA,IGAU)
  680.  101            CONTINUE
  681.  100          CONTINUE
  682.               CALL DEVOLU(XE,SHPWRK,MFR,NBNN,IFOUR,NIFOUR,IDIM,DIM3,
  683.      &                    RR,DJAC)
  684. c
  685. c     verification du signe du jacobien
  686. c
  687.               IF (DJAC.LT.0.) ISDJC=ISDJC+1
  688.               DJAC = ABS(DJAC)
  689.               IF (DJAC.LT.XPETIT) THEN
  690.                 INTERR(1) = IB
  691.                 CALL ERREUR(259)
  692.                 GOTO 9004
  693.               ENDIF
  694.               DJAC = DJAC * POIGAU(IGAU)
  695. c
  696. c     on recupere les contraintes
  697. c
  698.               MPTVAL=IVASTR
  699.               DO 5004 ICOMP=1,NSTR
  700.                 MELVAL=IVAL(ICOMP)
  701.                 IGMN = MIN(IGAU,VELCHE(/1))
  702.                 IBMN = MIN(IB  ,VELCHE(/2))
  703.                 XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  704.  5004         CONTINUE
  705. c
  706.               IF (IFOUR.EQ.1) THEN
  707.                 IF (NIFOUR.EQ.0) THEN
  708.                   CALL THSIG1(SHPWRK,DJAC,XSTRS,NBNN,LRE,REL,RR)
  709.                 ELSE
  710.                   CALL THSIG2(SHPWRK,DJAC,XSTRS,NBNN,LRE,REL,NIFOUR,RR)
  711.                 ENDIF
  712.               ELSE IF (IFOUR.EQ.0) THEN
  713.                 CALL THSIG3(SHPWRK,DJAC,XSTRS,NBNN,LRE,REL,RR)
  714.               ELSE
  715.                 CALL THSIGH(SHPWRK,DJAC,XSTRS,NBNN,IDIM,LRE,REL)
  716.               ENDIF
  717. c
  718.  4004       CONTINUE
  719.  
  720.             IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  721.               INTERR(1) = IB
  722.               CALL ERREUR(195)
  723.               GOTO 9004
  724.             ENDIF
  725. c
  726. c     remplissage de xmatri
  727. c
  728.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  729. c
  730.  3004     CONTINUE
  731.  
  732.  9004     CONTINUE
  733.           SEGSUP MWRK1,MWRK2
  734.           GOTO 510
  735. c
  736. c_______________________________________________________________________
  737. c
  738. ccccccccccccccccccc element coq3
  739. c_______________________________________________________________________
  740.  27       CONTINUE
  741.           NBBB = NBNN
  742.           SEGINI,MWRK1,MWRK3,MWRK4
  743.  
  744.           DO 3027 IB = 1, NBELEM
  745. c
  746. c     on cherche les coordonnees des noeuds de l element ib
  747. c
  748.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  749. c
  750.             CALL ZERO(REL,LRE,LRE)
  751. c
  752. c     on cherche les contraintes
  753. c
  754.             MPTVAL=IVASTR
  755.             DO 5027 ICOMP=1,NSTR
  756.               MELVAL=IVAL(ICOMP)
  757.               IBMN=MIN(IB  ,VELCHE(/2))
  758.               XSTRS(ICOMP)=VELCHE(1,IBMN)
  759.  5027       CONTINUE
  760. c
  761. ccccccc   on calcule k(sigma)
  762. c
  763.             CALL COQ3KS(REL,XSTRS,XE,1.D0,WORK)
  764. c
  765. c     remplissage de xmatri
  766. c
  767.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  768.  
  769.  3027     CONTINUE
  770.  
  771. C 9027     CONTINUE
  772.           SEGSUP,MWRK1,MWRK3,MWRK4
  773.           GOTO 510
  774. c
  775. c_______________________________________________________________________
  776. c
  777. c     element dkt , dst
  778. c_______________________________________________________________________
  779.  28       CONTINUE
  780.           DIM3 = 1.D0
  781.           NBNO = NBNN
  782.           IDI2=IDIM-1
  783.           NBBB=NBNN
  784.           SEGINI MWRK1,MWRK2,MWRK4,MWRK5
  785.           XX(1)=.5D0
  786.           XX(2)=.0D0
  787.           XX(3)=.5D0
  788.           YY(1)=.0D0
  789.           YY(2)=.5D0
  790.           YY(3)=.5D0
  791. C Pour la recuperation de l'epaisseur des elements DKT
  792.           IEPDKT = 0
  793.           IF (MFR.EQ.3 .AND. IFOUR.NE.-2) IEPDKT = IVACAR
  794. c*of 2011/06/22 : Quid de l'epaisseur pour les DST ????? EPAI = 0 ici !!
  795. C Pour la recuperation des axes d'orthotropie des elements DST
  796.           IAODST = 0
  797.           IF (MELE.EQ.93.AND.CMATE.NE.'ISOTROPE') IAODST = IVACAR
  798.  
  799.           DO 3028  IB=1,NBELEM
  800.  
  801. c     on cherche les coordonnees des noeuds de l element ib
  802. c
  803.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  804.  
  805.             CALL ZERO(REL,LRE,LRE)
  806.             CALL VPAST(XE,BPSS)
  807. c     bpss    stocke la matrice de passage
  808.             CALL VCORLC (XE,XEL,BPSS)
  809. c
  810. c     boucle sur les points de gauss
  811. c
  812.             DO 4028 IGAU=1,NBPGAU
  813. c
  814. c   recuperation de l'epaisseur (element DKT)
  815.               IF (IEPDKT.NE.0) THEN
  816.                 MPTVAL=IEPDKT
  817.                 MELVAL=IVAL(1)
  818.                 IGMN=MIN(IGAU,VELCHE(/1))
  819.                 IBMN=MIN(IB,VELCHE(/2))
  820.                 EPAI=VELCHE(IGMN,IBMN)
  821.               ELSE
  822.                 EPAI = XZERO
  823.               ENDIF
  824. c
  825.               call DKTSHP(IGAU,XEL,tabw,DJAC)
  826.               call GEOCST(XEL,GEOM)
  827.               call BBGFDK(XX(IGAU),YY(IGAU),GEOM,tabrot)
  828.  
  829.               DO 6028 IC=1,NBNN
  830.                 DO 60281 ID=1,6
  831.                   SHPWRK(ID,IC)=SHPTOT(ID,IC,IGAU)
  832. 60281           CONTINUE
  833.  6028         CONTINUE
  834.  
  835.               CALL DEVOLU(XEL,SHPWRK,MFR,NBNN,IFOUR,NIFOUR,IDI2,DIM3,
  836.      &                    RR,DJAC)
  837.               DJAC=DJAC*POIGAU(IGAU)
  838. c
  839. c     on cherche les contraintes
  840. c
  841.               MPTVAL=IVASTR
  842.               DO 5028 ICOMP=1,NSTRS
  843.                 MELVAL=IVAL(ICOMP)
  844.                 IGMN=MIN(IGAU,VELCHE(/1))
  845.                 IBMN=MIN(IB  ,VELCHE(/2))
  846.                 XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  847. C                 write(6,*)' xstrs(icomp)',icomp,XSTRS(ICOMP)
  848.  5028         CONTINUE
  849.  
  850. C Recuperation des axes d'orthotropie (element DST)
  851.               IF (IAODST.NE.0) THEN
  852.                 MPTVAL=IAODST
  853.                 MELVAL=IVAL(1)
  854.                 IBMN=MIN(IB  ,VELCHE(/2))
  855.                 IGMN=MIN(IGAU,VELCHE(/1))
  856.                 COSA=VELCHE(IGMN,IBMN)
  857.                 MELVAL=IVAL(2)
  858.                 IBMN=MIN(IB  ,VELCHE(/2))
  859.                 IGMN=MIN(IGAU,VELCHE(/1))
  860.                 SINA=VELCHE(IGMN,IBMN)
  861.                 CC=COSA*COSA
  862.                 SS=SINA*SINA
  863.                 CS=SINA*COSA
  864. C
  865. C        chgt d'axes
  866. C
  867.                 SIG1=CC*XSTRS(1)+SS*XSTRS(2)-2.D0*CS*XSTRS(3)
  868.                 SIG2=CC*XSTRS(2)+SS*XSTRS(1)+2.D0*CS*XSTRS(3)
  869.                 SIG3=CS*(XSTRS(1)-XSTRS(2))+(CC-SS)*XSTRS(3)
  870.                 XSTRS(1)=SIG1
  871.                 XSTRS(2)=SIG2
  872.                 XSTRS(3)=SIG3
  873.               ENDIF
  874. c
  875.               CALL DKTHSH(SHPWRK,tabw,tabrot,DJAC,XSTRS,REL,EPAI)
  876.  4028       CONTINUE
  877.  
  878.             CALL TRANSK(REL,BPSS,LRE,3,1)
  879. c
  880. c     remplissage de xmatri
  881. c
  882.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  883. c
  884.  3028     CONTINUE
  885.  
  886. C 9028     CONTINUE
  887.           SEGSUP,MWRK1,MWRK2,MWRK4,MWRK5
  888.           GOTO 510
  889.  
  890. c_______________________________________________________________________
  891. c
  892. c  element poutre
  893. c_______________________________________________________________________
  894.  29       CONTINUE
  895.           NBBB = NBNN
  896.           SEGINI,MWRK1,MWRK3
  897.  
  898.           DO 3029 IB=1,NBELEM
  899. c
  900. c  on cherche les coordonnees des noeuds de l elementib
  901. c
  902.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  903. c
  904. c  il faudrait aussi modifier le vecteur local de la poutre
  905. c
  906. c  mise a zero de la raideur geometrique
  907. c
  908.             CALL ZERO(REL,LRE,LRE)
  909. c
  910. c     rangement des caracteristiques dans work
  911. c
  912.             MPTVAL=IVACAR
  913.             DO 6029 IC=1,NCARR
  914.               WORK(IC)=XZERO
  915.               IF (IVAL(IC).NE.0) THEN
  916.                 MELVAL=IVAL(IC)
  917.                 IBMN=MIN(IB,VELCHE(/2))
  918.                 DO 4029 IGAU=1,NBNN
  919.                   IGMN=MIN(IGAU,VELCHE(/1))
  920.                   IF (IGMN.GT.0.AND.IBMN.GT.0) THEN
  921.                     WORK(IC)=WORK(IC)+VELCHE(IGMN,IBMN)
  922.                   ENDIF
  923.  4029           CONTINUE
  924.                 WORK(IC)=WORK(IC)/NBNN
  925.               ENDIF
  926.  6029       CONTINUE
  927. c
  928. c  cas des tuyaux  -  on calcule les caracteristiques de la poutre
  929. c                     equivalente
  930. c
  931.             IF (MELE.EQ.42) THEN
  932. **            do ic=5,ncarr-1
  933. **            work(ic)=work(ic+1)
  934. **            enddo
  935. **            write (6,*) 'work isous ib',isous,ib
  936. **            write(6,*) (work(ic),ic=1,ncarr)
  937.  
  938.  
  939.  
  940.               CISA=WORK(4)
  941.               VX=WORK(5)
  942.               VY=WORK(6)
  943.               VZ=WORK(7)
  944. **            write(6,*) 'ksigmp vx vy vz',vx,vy,vz,isous,nbelem
  945.               CALL TUYCAR(WORK,CISA,VX,VY,VZ,KERRE,0)
  946.               IF (KERRE.EQ.77) THEN
  947.                 CALL ERREUR(77)
  948.                 GOTO 9029
  949.               ENDIF
  950.             ENDIF
  951. c
  952. c  on cherche les contraintes - on les met dans work
  953. c
  954.             IE = 9
  955.             MPTVAL=IVASTR
  956.             DO 7029 ID=1,2
  957.               ID2=ID
  958.               IF (NBPGAU.EQ.1.AND.ID.EQ.2) ID2=1
  959.               DO 70291 ICOMP=1,NSTR
  960.                 IE = IE+1
  961.                 MELVAL=IVAL(ICOMP)
  962.                 IGMN=MIN(ID2 ,VELCHE(/1))
  963.                 IBMN=MIN(IB  ,VELCHE(/2))
  964.                 WORK(IE)=VELCHE(IGMN,IBMN)
  965. 70291         CONTINUE
  966.  7029       CONTINUE
  967. c
  968. c  on calcule la rigidite geometrique
  969. c
  970.             IF (IFOUR.EQ.-2.OR.IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN
  971.               CALL POUKS2(REL,LRE,WORK(10),WORK,XE,WORK(22),KERRE)
  972.             ELSE
  973.               CALL POUKSG(REL,LRE,WORK(10),WORK,XE,WORK(22),KERRE)
  974.             ENDIF
  975.  
  976.             IF (KERRE.NE.0) THEN
  977.               INTERR(1)=ISOUS
  978.               INTERR(2)=IB
  979.               CALL ERREUR(128)
  980.               GOTO 9029
  981.             ENDIF
  982. c
  983. c  remplissage de xmatri
  984. c
  985.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  986. C
  987.  3029     CONTINUE
  988. c
  989.  9029     CONTINUE
  990.           SEGSUP,MWRK1,MWRK3
  991. c
  992.           GOTO 510
  993. c_______________________________________________________________________
  994. c
  995. c  elements coq8 et coq6
  996. c_______________________________________________________________________
  997.  41       CONTINUE
  998.           NBBB=NBNN
  999.           LRI =NBNN*5
  1000.           SEGINI,MWRK1,MWRK3
  1001. c
  1002.           MINTE1 = IPMIN1
  1003.           SEGACT,MINTE1
  1004. c
  1005.           DO 3041 IB=1,NBELEM
  1006. c
  1007. c  on cherche les coordonnees des noeuds de l elementib
  1008. c
  1009.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1010.  
  1011.             CALL ZERO(REL,LRE,LRE)
  1012. c
  1013. c  on cherche les caracteristiques de l element ib
  1014. c
  1015.             MPTVAL = IVACAR
  1016.             MELVAL = IVAL(1)
  1017.             IF (MELVAL.NE.0) THEN
  1018.               IBMN = MIN(IB  ,VELCHE(/2))
  1019.               DO IGAU = 1, NBNN
  1020.                 IGMN = MIN(IGAU,VELCHE(/1))
  1021.                 WORK(IGAU) = VELCHE(IGMN,IBMN)
  1022.               ENDDO
  1023.             ELSE
  1024.               DO IGAU = 1, NBNN
  1025.                 WORK(IGAU)=XZERO
  1026.               ENDDO
  1027.             ENDIF
  1028. c
  1029. c  on cherche les contraintes - on les met dans work
  1030. c
  1031.             IE = 9
  1032.             MPTVAL=IVASTR
  1033.             DO 7041 IGAU=1,NBPGAU
  1034.               DO 7042 ICOMP=1,NSTRS
  1035.                 MELVAL=IVAL(ICOMP)
  1036.                 IGMN=MIN(IGAU,VELCHE(/1))
  1037.                 IBMN=MIN(IB  ,VELCHE(/2))
  1038.                 WORK(IE)=VELCHE(IGMN,IBMN)
  1039.                 IE=IE+1
  1040.  7042         CONTINUE
  1041.  7041       CONTINUE
  1042. c
  1043. c  on calcule la rigidite geometrique
  1044. c
  1045.             CALL COQ8KS(REL,XE,SHPTOT,MINTE1.SHPTOT,
  1046.      &                  NBPGAU,POIGAU,DZEGAU,
  1047.      &                  WORK(1),WORK(9),NBNN,LRE,LRI,WORK(51))
  1048. c
  1049. c  remplissage de xmatri
  1050. c
  1051.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1052. C
  1053.  3041     CONTINUE
  1054.  
  1055. C 9041     CONTINUE
  1056.           SEGSUP,MWRK1,MWRK3
  1057.           GO TO 510
  1058. c_______________________________________________________________________
  1059. c
  1060. c     tuyau fissure
  1061. c_______________________________________________________________________
  1062.  43       CONTINUE
  1063. c     ksigma n a pas de sens evident pour cet element
  1064. c     on cree une matrice nulle
  1065. c          DO 3043 IB=1,NBELEM
  1066. c            do 4043 ic=1,lval
  1067. c              re(ic,ic,ib)=XZERO
  1068. c 4043       continue
  1069. c 3043     CONTINUE
  1070.           GOTO 510
  1071. c
  1072. c_______________________________________________________________________
  1073. c
  1074. c     element coq2
  1075. c_______________________________________________________________________
  1076. c
  1077.  44       CONTINUE
  1078.           DIM3=1.D0
  1079.           NBBB=NBNN
  1080.           SEGINI MWRK1,MWRK3,MWRK4
  1081. c
  1082.           DO 3044 IB=1,NBELEM
  1083. c
  1084. c     on cherche les coordonnees des noeuds de  l element ib
  1085. c
  1086.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1087.  
  1088.        IF (IFOUR.EQ.1.AND.IDIM.EQ.3) THEN
  1089. c jk148537 assume 1->r 3->z 
  1090.          do in=1,nbnn
  1091.            xe(2,in) = xe(3,in)
  1092.          enddo
  1093.        ENDIF
  1094.  
  1095.             CALL ZERO(REL,LRE,LRE)
  1096. c
  1097. c   recuperation de l'epaisseur
  1098. c
  1099.             IF (IFOUR.EQ.-2.AND.IPCHE2.NE.0) THEN
  1100.               MPTVAL=IVACAR
  1101.               MELVAL=IVAL(1)
  1102.               IF (MELVAL.NE.0) THEN
  1103.                 IBMN=MIN(IB,VELCHE(/2))
  1104.                 DIM3=VELCHE(1,IBMN)
  1105.               ELSE
  1106.                 DIM3 = 1.D0
  1107.               ENDIF
  1108.             ENDIF
  1109. c
  1110. c     on cherche les contraintes on les met dans work...
  1111. c
  1112.             JC = 0
  1113.             MPTVAL=IVASTR
  1114.             DO 5044 IGAU=1,NBPGAU
  1115.               DO 5045 ICOMP=1,NSTRS
  1116.                 MELVAL=IVAL(ICOMP)
  1117.                 IGMN=MIN(IGAU,VELCHE(/1))
  1118.                 IBMN=MIN(IB  ,VELCHE(/2))
  1119.                 JC=JC+1
  1120.                 WORK(JC)=VELCHE(IGMN,IBMN)
  1121.  5045         CONTINUE
  1122.  5044       CONTINUE
  1123. c
  1124. c     appel a coque2 ksigma...
  1125. c
  1126.             AN=NHRM
  1127.             CALL CQ2KSG(XE,1.D0,DIM3,IFOUR,AN,NBPGAU,WORK(1),WORK(19),
  1128.      1                  WORK(22),QSIGAU,POIGAU,WORK(25),WORK(30),
  1129.      2                  WORK(35),WORK(42),WORK(49),WORK(113),WORK(177),
  1130.      3                  WORK(241),WORK(305),LRE,REL)
  1131. c
  1132. c     remplissage de xmatri
  1133. c
  1134.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1135. c
  1136.  3044     CONTINUE
  1137. c
  1138. C 9044     CONTINUE
  1139.           SEGSUP,MWRK1,MWRK3,MWRK4
  1140.           GOTO 510
  1141. c
  1142. c_______________________________________________________________________
  1143. c
  1144. c  elements barre et cercle (et TIMO)
  1145. c_______________________________________________________________________
  1146.  46       CONTINUE
  1147. C Cas particulier :
  1148.           IF (MELE.EQ.95.AND.IFOUR.NE.0.AND.IFOUR.NE.1) GOTO 99
  1149. C
  1150.           NBBB = NBNN
  1151.           SEGINI MWRK1,MWRK3
  1152. c
  1153.           DO 3046 IB=1,NBELEM
  1154. c
  1155. c  on cherche les coordonnees des noeuds de l elementib
  1156. c
  1157.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1158. c
  1159. c  mise a zero de la raideur geometrique
  1160. c
  1161.             CALL ZERO(REL,LRE,LRE)
  1162. c
  1163. c  on cherche l'effort
  1164. c
  1165.             MPTVAL=IVASTR
  1166.             MELVAL=IVAL(1)
  1167.             NBPTEL=VELCHE(/1)
  1168.             IBMN=MIN(IB,VELCHE(/2))
  1169. c
  1170.             IF (NBPTEL.EQ.1) THEN
  1171.               EFFORT=VELCHE(1,IBMN)
  1172.             ELSE IF (NBPTEL.EQ.2) THEN
  1173.               EFF1=VELCHE(1,IBMN)
  1174.               EFF2=VELCHE(2,IBMN)
  1175.               EFFORT=0.5D0*(EFF1+EFF2)
  1176.             ENDIF
  1177. c
  1178. c  on calcule la rigidite geometrique
  1179. c
  1180.             IF (MELE.EQ.46.or.MELE.eq.84)
  1181.      &         CALL BARKSG(REL,LRE,EFFORT,XE,KERRE)
  1182.             IF (MELE.EQ.95) CALL CERKSG(REL,LRE,EFFORT,XE,KERRE)
  1183.             IF (KERRE.NE.0) THEN
  1184.               INTERR(1)=ISOUS
  1185.               INTERR(2)=IB
  1186.               CALL ERREUR(128)
  1187.               GO TO 9046
  1188.             ENDIF
  1189. c
  1190. c  remplissage de xmatri
  1191. c
  1192. c           cas particulier TIMO : on saute les ddls de rotation
  1193.             IF (MELE.EQ.84) THEN
  1194.               NCOMPU=NCOMP/2
  1195.               ii=0
  1196.               iii=0
  1197.               DO 841 INOEUD=1,NBNNS
  1198.               DO 842 ICOMP=1,NCOMP
  1199.                 ii=ii+1
  1200.                 if(ii.gt.NCOMPU) goto 842
  1201.                 iii=iii+1
  1202.                 jj=0
  1203.                 jjj=0
  1204.                 DO 843 JNOEUD=1,NBNNS
  1205.                 DO 844 JCOMP=1,NCOMP
  1206.                   jj=jj+1
  1207.                   if(jj.gt.ii) goto 842
  1208.                   if(jj.gt.NCOMPU) goto 844
  1209.                   jjj=jjj+1
  1210.                   RE(ii,jj,ib)=REL(iii,jjj)
  1211.                   RE(jj,ii,ib)=REL(iii,jjj)
  1212.  844            CONTINUE
  1213.  843            CONTINUE
  1214.  842          CONTINUE
  1215.  841          CONTINUE
  1216.             ELSE
  1217.               CALL REMPMT(REL,LRE,RE(1,1,ib))
  1218.             ENDIF
  1219. C
  1220.  3046     CONTINUE
  1221.  
  1222.  9046     CONTINUE
  1223.           SEGSUP,MWRK1,MWRK3
  1224.           GOTO 510
  1225.  
  1226. c_______________________________________________________________________
  1227. c
  1228. c     element coq4
  1229. c_______________________________________________________________________
  1230.  49       CONTINUE
  1231.           NBBB=NBNN
  1232.           NBNO=NBNN
  1233.           SEGINI,MWRK1,MWRK2,MWRK4
  1234.  
  1235.           DO 3049 IB=1,NBELEM
  1236. c
  1237. c     on cherche les coordonnees des noeuds de  l element ib
  1238. c
  1239.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1240.  
  1241.             CALL ZERO(REL,LRE,LRE)
  1242.             CALL CQ4LOC(XE,XEL,BPSS,IRRT,0)
  1243. C
  1244. C   attention  :  rien de prevu en cas d'excentrement
  1245. C
  1246.             CALL BCOQ4(5,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XZERO,0,IRRT,0)
  1247. c
  1248.             MPTVAL=IVASTR
  1249.             DO 5049 ICOMP=1,NSTRS
  1250.               MELVAL=IVAL(ICOMP)
  1251.               IGMN=MIN(5,VELCHE(/1))
  1252.               IBMN=MIN(IB  ,VELCHE(/2))
  1253.               XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  1254.  5049       CONTINUE
  1255. c
  1256.             CALL CQ4KSG(DJAC,XSTRS,SHPWRK, REL)
  1257.             CALL TRANSK(REL,BPSS,LRE,4,0)
  1258. c
  1259. c     remplissage de xmatri
  1260. c
  1261.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1262. C
  1263.  3049     CONTINUE
  1264.  
  1265. C 9049     CONTINUE
  1266.           SEGSUP,MWRK1,MWRK2,MWRK4
  1267.           GOTO 510
  1268.  
  1269. c_______________________________________________________________________
  1270. c
  1271. c     element cof3
  1272. c_______________________________________________________________________
  1273. c
  1274.  51       CONTINUE
  1275. c
  1276.           NBBB=NBNN
  1277.           SEGINI,MWRK1,MWRK3,MWRK4
  1278. c
  1279.           CALL ERREUR(19)
  1280.           GOTO 9051
  1281.  
  1282. C          DO 3051 IB=1,NBELEM
  1283. Cc
  1284. Cc     on cherche les coordonnees des noeuds de  l element ib
  1285. Cc
  1286. C            CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1287. C
  1288. C            CALL ZERO(REL,LRE,LRE)
  1289. C
  1290. C            MPTVAL=IVACAR
  1291. C            MELVAL=IVAL(1)
  1292. C            IF (MELVAL.NE.0) THEN
  1293. C              IBMN=MIN(IB  ,VELCHE(/2))
  1294. C              EPAI=VELCHE(1,IBMN)
  1295. C            ELSE
  1296. C              EPAI=XZERO
  1297. C            ENDIF
  1298. Cc
  1299. Cc     on cherche les contraintes on les met dans work...
  1300. Cc
  1301. C            JC=0
  1302. C            MPTVAL=IVASTR
  1303. C            DO 5051 IGAU=1,NBPGAU
  1304. C              DO 5051 ICOMP=1,NSTRS
  1305. C                MELVAL=IVAL(ICOMP)
  1306. C                IGMN=MIN(IGAU,VELCHE(/1))
  1307. C                IBMN=MIN(IB  ,VELCHE(/2))
  1308. C                JC=JC+1
  1309. C                WORK(JC)=VELCHE(IGMN,IBMN)
  1310. C 5051       CONTINUE
  1311. Cc
  1312. Cc     appel a coque2 ksigma...
  1313. Cc
  1314. C            AN=NHRM
  1315. CC             call cq3ksg(xe,epai,an,nbpgau,work(1),work(19),work(22),
  1316. CC     1                   work(25),work(30),work(35),work(42),work(49),
  1317. CC     2                   work(113),work(177),work(241),work(305),rel)
  1318. Cc
  1319. Cc     remplissage de xmatri
  1320. Cc
  1321. C            CALL REMPMT(REL,LRE,RE(1,1,ib))
  1322. Cc
  1323. C 3051     CONTINUE
  1324.  
  1325.  9051     CONTINUE
  1326.           SEGSUP MWRK1,MWRK3,MWRK4
  1327.           GOTO 510
  1328.  
  1329. c_______________________________________________________________________
  1330. c
  1331. c     element shb8
  1332. c_______________________________________________________________________
  1333.  260      CONTINUE
  1334.           NBBB=NBNN
  1335.           SEGINI,MWRK1,MWRK7
  1336. C           write(6,*) ' nbnn nbpgau nstrs lre' , NBNN,nbpgau,nstrs,lre
  1337.  
  1338.           DO 3260 IB=1,NBELEM
  1339. c
  1340. c     on cherche les coordonnees des noeuds de  l element ib
  1341. c
  1342.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1343.  
  1344.             MPTVAL=IVASTR
  1345.             IE=0
  1346.             do 3268 igau=1,nbpgau
  1347.               DO 3269 ICOMP=1,NSTRS
  1348.                 iE=IE+1
  1349.                 MELVAL=IVAL(ICOMP)
  1350.                 IGMN=MIN(IGAU,VELCHE(/1))
  1351.                 IBMN=MIN(IB  ,VELCHE(/2))
  1352.                 work1(ie)=VELCHE(IGMN,IBMN)
  1353. C                 write(6,*)' xstrs(icomp)',icomp,XSTRS(ICOMP)
  1354.  3269         CONTINUE
  1355.  3268       CONTINUE
  1356.             propel(1)=0.
  1357.             call shb8 (9,xe,D,propel,work1,rel,out)
  1358.  
  1359. C      remplissage de xmatri
  1360.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1361. C
  1362.  3260     CONTINUE
  1363. c
  1364. C 9260     CONTINUE
  1365.           SEGSUP,MWRK1,MWRK7
  1366.           GOTO 510
  1367.  
  1368. c_______________________________________________________________________
  1369. c
  1370. c     desactivation des segments propres a la zone geometrique isous
  1371. c_______________________________________________________________________
  1372. c
  1373.  510      CONTINUE
  1374.  597      CONTINUE
  1375.           IF (ISUP1L.EQ.1 .OR. nblprt.GT.1) THEN
  1376.             CALL DTMVAL(IVASTR,3)
  1377.           ELSE
  1378.             CALL DTMVAL(IVASTR,1)
  1379.           ENDIF
  1380.           IF (ISUP2L.EQ.1 .OR. nblprt.GT.1) THEN
  1381.             CALL DTMVAL(IVACAR,3)
  1382.           ELSE
  1383.             CALL DTMVAL(IVACAR,1)
  1384.           ENDIF
  1385.           xmatri = ipmatr
  1386.           SEGDES,xmatri
  1387.  
  1388. C- Sortie prematuree en cas d'erreur
  1389.           IF (IERR.NE.0) GOTO 598
  1390.  
  1391. C- Stockage de la matrice
  1392.           nrigel=irigel(/2) +1
  1393.           segadj,mrigid
  1394. C           jrige = NRIGE0 + irige
  1395.           jrige=nrigel
  1396.           COERIG(jrige)   = 1.
  1397.           IRIGEL(1,jrige) = ipmail
  1398.           IRIGEL(2,jrige) = 0
  1399.           IRIGEL(3,jrige) = ipdesc
  1400.           IRIGEL(4,jrige) = ipmatr
  1401.           IRIGEL(5,jrige) = NIFOUR
  1402.           IRIGEL(6,jrige) = 0
  1403.           IRIGEL(7,jrige) = 0
  1404.           IRIGEL(8,jrige) = 0
  1405.  
  1406.         ENDDO
  1407. C- Fin de la boucle sur les partitions
  1408. C
  1409.  598    CONTINUE
  1410.         IF (MOSTRS.NE.0) THEN
  1411.           nomid = MOSTRS
  1412.           IF (lsupco) SEGSUP,nomid
  1413.           notype = MOTYPS
  1414.           SEGSUP,notype
  1415.         ENDIF
  1416.         IF (MOCARA.NE.0) THEN
  1417.           NOMID = MOCARA
  1418.           SEGSUP,NOMID
  1419.           notype = MOTYPC
  1420.           SEGSUP,notype
  1421.         ENDIF
  1422. C
  1423.         NOMID=MODEPL
  1424.         IF (lsupde) SEGSUP,NOMID
  1425.         NOMID = MOFORC
  1426.         IF (lsupfo) SEGSUP,NOMID
  1427.  
  1428.  599    CONTINUE
  1429.  
  1430.         IF (IERR.NE.0) GOTO 999
  1431. C
  1432.  500  CONTINUE
  1433. C* Fin de la boucle sur les modeles elementaires
  1434.  
  1435.  999  CONTINUE
  1436.       IF (IERR.NE.0) THEN
  1437.         ktrace = -1
  1438. **  En situation d'erreur, on laisse le menage faire son travail
  1439. **      CALL DERIGI(MRIGID,ktrace,msorse)
  1440. **      SEGSUP,MRIGID
  1441.         IPRIGG = 0
  1442.       ELSE
  1443.         if(irigel(/2).eq.0) then
  1444.           call erreur (86)
  1445.           return
  1446.         endif
  1447.         SEGDES,MRIGID
  1448.         IPRIGG = MRIGID
  1449.       ENDIF
  1450.  
  1451. c      RETURN
  1452.       END
  1453.  
  1454.  
  1455.  
  1456.  
  1457.  
  1458.  
  1459.  

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