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    PV090527  26/04/30    21:15:47     12529          
  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.             CALL ZERO(REL,LRE,LRE)
  1089. c
  1090. c   recuperation de l'epaisseur
  1091. c
  1092.             IF (IFOUR.EQ.-2.AND.IPCHE2.NE.0) THEN
  1093.               MPTVAL=IVACAR
  1094.               MELVAL=IVAL(1)
  1095.               IF (MELVAL.NE.0) THEN
  1096.                 IBMN=MIN(IB,VELCHE(/2))
  1097.                 DIM3=VELCHE(1,IBMN)
  1098.               ELSE
  1099.                 DIM3 = 1.D0
  1100.               ENDIF
  1101.             ENDIF
  1102. c
  1103. c     on cherche les contraintes on les met dans work...
  1104. c
  1105.             JC = 0
  1106.             MPTVAL=IVASTR
  1107.             DO 5044 IGAU=1,NBPGAU
  1108.               DO 5045 ICOMP=1,NSTRS
  1109.                 MELVAL=IVAL(ICOMP)
  1110.                 IGMN=MIN(IGAU,VELCHE(/1))
  1111.                 IBMN=MIN(IB  ,VELCHE(/2))
  1112.                 JC=JC+1
  1113.                 WORK(JC)=VELCHE(IGMN,IBMN)
  1114.  5045         CONTINUE
  1115.  5044       CONTINUE
  1116. c
  1117. c     appel a coque2 ksigma...
  1118. c
  1119.             AN=NHRM
  1120.             CALL CQ2KSG(XE,1.D0,DIM3,IFOUR,AN,NBPGAU,WORK(1),WORK(19),
  1121.      1                  WORK(22),QSIGAU,POIGAU,WORK(25),WORK(30),
  1122.      2                  WORK(35),WORK(42),WORK(49),WORK(113),WORK(177),
  1123.      3                  WORK(241),WORK(305),LRE,REL)
  1124. c
  1125. c     remplissage de xmatri
  1126. c
  1127.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1128. c
  1129.  3044     CONTINUE
  1130. c
  1131. C 9044     CONTINUE
  1132.           SEGSUP,MWRK1,MWRK3,MWRK4
  1133.           GOTO 510
  1134. c
  1135. c_______________________________________________________________________
  1136. c
  1137. c  elements barre et cercle (et TIMO)
  1138. c_______________________________________________________________________
  1139.  46       CONTINUE
  1140. C Cas particulier :
  1141.           IF (MELE.EQ.95.AND.IFOUR.NE.0.AND.IFOUR.NE.1) GOTO 99
  1142. C
  1143.           NBBB = NBNN
  1144.           SEGINI MWRK1,MWRK3
  1145. c
  1146.           DO 3046 IB=1,NBELEM
  1147. c
  1148. c  on cherche les coordonnees des noeuds de l elementib
  1149. c
  1150.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1151. c
  1152. c  mise a zero de la raideur geometrique
  1153. c
  1154.             CALL ZERO(REL,LRE,LRE)
  1155. c
  1156. c  on cherche l'effort
  1157. c
  1158.             MPTVAL=IVASTR
  1159.             MELVAL=IVAL(1)
  1160.             NBPTEL=VELCHE(/1)
  1161.             IBMN=MIN(IB,VELCHE(/2))
  1162. c
  1163.             IF (NBPTEL.EQ.1) THEN
  1164.               EFFORT=VELCHE(1,IBMN)
  1165.             ELSE IF (NBPTEL.EQ.2) THEN
  1166.               EFF1=VELCHE(1,IBMN)
  1167.               EFF2=VELCHE(2,IBMN)
  1168.               EFFORT=0.5D0*(EFF1+EFF2)
  1169.             ENDIF
  1170. c
  1171. c  on calcule la rigidite geometrique
  1172. c
  1173.             IF (MELE.EQ.46.or.MELE.eq.84)
  1174.      &         CALL BARKSG(REL,LRE,EFFORT,XE,KERRE)
  1175.             IF (MELE.EQ.95) CALL CERKSG(REL,LRE,EFFORT,XE,KERRE)
  1176.             IF (KERRE.NE.0) THEN
  1177.               INTERR(1)=ISOUS
  1178.               INTERR(2)=IB
  1179.               CALL ERREUR(128)
  1180.               GO TO 9046
  1181.             ENDIF
  1182. c
  1183. c  remplissage de xmatri
  1184. c
  1185. c           cas particulier TIMO : on saute les ddls de rotation
  1186.             IF (MELE.EQ.84) THEN
  1187.               NCOMPU=NCOMP/2
  1188.               ii=0
  1189.               iii=0
  1190.               DO 841 INOEUD=1,NBNNS
  1191.               DO 842 ICOMP=1,NCOMP
  1192.                 ii=ii+1
  1193.                 if(ii.gt.NCOMPU) goto 842
  1194.                 iii=iii+1
  1195.                 jj=0
  1196.                 jjj=0
  1197.                 DO 843 JNOEUD=1,NBNNS
  1198.                 DO 844 JCOMP=1,NCOMP
  1199.                   jj=jj+1
  1200.                   if(jj.gt.ii) goto 842
  1201.                   if(jj.gt.NCOMPU) goto 844
  1202.                   jjj=jjj+1
  1203.                   RE(ii,jj,ib)=REL(iii,jjj)
  1204.                   RE(jj,ii,ib)=REL(iii,jjj)
  1205.  844            CONTINUE
  1206.  843            CONTINUE
  1207.  842          CONTINUE
  1208.  841          CONTINUE
  1209.             ELSE
  1210.               CALL REMPMT(REL,LRE,RE(1,1,ib))
  1211.             ENDIF
  1212. C
  1213.  3046     CONTINUE
  1214.  
  1215.  9046     CONTINUE
  1216.           SEGSUP,MWRK1,MWRK3
  1217.           GOTO 510
  1218.  
  1219. c_______________________________________________________________________
  1220. c
  1221. c     element coq4
  1222. c_______________________________________________________________________
  1223.  49       CONTINUE
  1224.           NBBB=NBNN
  1225.           NBNO=NBNN
  1226.           SEGINI,MWRK1,MWRK2,MWRK4
  1227.  
  1228.           DO 3049 IB=1,NBELEM
  1229. c
  1230. c     on cherche les coordonnees des noeuds de  l element ib
  1231. c
  1232.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1233.  
  1234.             CALL ZERO(REL,LRE,LRE)
  1235.             CALL CQ4LOC(XE,XEL,BPSS,IRRT,0)
  1236. C
  1237. C   attention  :  rien de prevu en cas d'excentrement
  1238. C
  1239.             CALL BCOQ4(5,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XZERO,0,IRRT,0)
  1240. c
  1241.             MPTVAL=IVASTR
  1242.             DO 5049 ICOMP=1,NSTRS
  1243.               MELVAL=IVAL(ICOMP)
  1244.               IGMN=MIN(5,VELCHE(/1))
  1245.               IBMN=MIN(IB  ,VELCHE(/2))
  1246.               XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  1247.  5049       CONTINUE
  1248. c
  1249.             CALL CQ4KSG(DJAC,XSTRS,SHPWRK, REL)
  1250.             CALL TRANSK(REL,BPSS,LRE,4,0)
  1251. c
  1252. c     remplissage de xmatri
  1253. c
  1254.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1255. C
  1256.  3049     CONTINUE
  1257.  
  1258. C 9049     CONTINUE
  1259.           SEGSUP,MWRK1,MWRK2,MWRK4
  1260.           GOTO 510
  1261.  
  1262. c_______________________________________________________________________
  1263. c
  1264. c     element cof3
  1265. c_______________________________________________________________________
  1266. c
  1267.  51       CONTINUE
  1268. c
  1269.           NBBB=NBNN
  1270.           SEGINI,MWRK1,MWRK3,MWRK4
  1271. c
  1272.           CALL ERREUR(19)
  1273.           GOTO 9051
  1274.  
  1275. C          DO 3051 IB=1,NBELEM
  1276. Cc
  1277. Cc     on cherche les coordonnees des noeuds de  l element ib
  1278. Cc
  1279. C            CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1280. C
  1281. C            CALL ZERO(REL,LRE,LRE)
  1282. C
  1283. C            MPTVAL=IVACAR
  1284. C            MELVAL=IVAL(1)
  1285. C            IF (MELVAL.NE.0) THEN
  1286. C              IBMN=MIN(IB  ,VELCHE(/2))
  1287. C              EPAI=VELCHE(1,IBMN)
  1288. C            ELSE
  1289. C              EPAI=XZERO
  1290. C            ENDIF
  1291. Cc
  1292. Cc     on cherche les contraintes on les met dans work...
  1293. Cc
  1294. C            JC=0
  1295. C            MPTVAL=IVASTR
  1296. C            DO 5051 IGAU=1,NBPGAU
  1297. C              DO 5051 ICOMP=1,NSTRS
  1298. C                MELVAL=IVAL(ICOMP)
  1299. C                IGMN=MIN(IGAU,VELCHE(/1))
  1300. C                IBMN=MIN(IB  ,VELCHE(/2))
  1301. C                JC=JC+1
  1302. C                WORK(JC)=VELCHE(IGMN,IBMN)
  1303. C 5051       CONTINUE
  1304. Cc
  1305. Cc     appel a coque2 ksigma...
  1306. Cc
  1307. C            AN=NHRM
  1308. CC             call cq3ksg(xe,epai,an,nbpgau,work(1),work(19),work(22),
  1309. CC     1                   work(25),work(30),work(35),work(42),work(49),
  1310. CC     2                   work(113),work(177),work(241),work(305),rel)
  1311. Cc
  1312. Cc     remplissage de xmatri
  1313. Cc
  1314. C            CALL REMPMT(REL,LRE,RE(1,1,ib))
  1315. Cc
  1316. C 3051     CONTINUE
  1317.  
  1318.  9051     CONTINUE
  1319.           SEGSUP MWRK1,MWRK3,MWRK4
  1320.           GOTO 510
  1321.  
  1322. c_______________________________________________________________________
  1323. c
  1324. c     element shb8
  1325. c_______________________________________________________________________
  1326.  260      CONTINUE
  1327.           NBBB=NBNN
  1328.           SEGINI,MWRK1,MWRK7
  1329. C           write(6,*) ' nbnn nbpgau nstrs lre' , NBNN,nbpgau,nstrs,lre
  1330.  
  1331.           DO 3260 IB=1,NBELEM
  1332. c
  1333. c     on cherche les coordonnees des noeuds de  l element ib
  1334. c
  1335.             CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1336.  
  1337.             MPTVAL=IVASTR
  1338.             IE=0
  1339.             do 3268 igau=1,nbpgau
  1340.               DO 3269 ICOMP=1,NSTRS
  1341.                 iE=IE+1
  1342.                 MELVAL=IVAL(ICOMP)
  1343.                 IGMN=MIN(IGAU,VELCHE(/1))
  1344.                 IBMN=MIN(IB  ,VELCHE(/2))
  1345.                 work1(ie)=VELCHE(IGMN,IBMN)
  1346. C                 write(6,*)' xstrs(icomp)',icomp,XSTRS(ICOMP)
  1347.  3269         CONTINUE
  1348.  3268       CONTINUE
  1349.             propel(1)=0.
  1350.             call shb8 (9,xe,D,propel,work1,rel,out)
  1351.  
  1352. C      remplissage de xmatri
  1353.             CALL REMPMT(REL,LRE,RE(1,1,ib))
  1354. C
  1355.  3260     CONTINUE
  1356. c
  1357. C 9260     CONTINUE
  1358.           SEGSUP,MWRK1,MWRK7
  1359.           GOTO 510
  1360.  
  1361. c_______________________________________________________________________
  1362. c
  1363. c     desactivation des segments propres a la zone geometrique isous
  1364. c_______________________________________________________________________
  1365. c
  1366.  510      CONTINUE
  1367.  597      CONTINUE
  1368.           IF (ISUP1L.EQ.1 .OR. nblprt.GT.1) THEN
  1369.             CALL DTMVAL(IVASTR,3)
  1370.           ELSE
  1371.             CALL DTMVAL(IVASTR,1)
  1372.           ENDIF
  1373.           IF (ISUP2L.EQ.1 .OR. nblprt.GT.1) THEN
  1374.             CALL DTMVAL(IVACAR,3)
  1375.           ELSE
  1376.             CALL DTMVAL(IVACAR,1)
  1377.           ENDIF
  1378.           xmatri = ipmatr
  1379.           SEGDES,xmatri
  1380.  
  1381. C- Sortie prematuree en cas d'erreur
  1382.           IF (IERR.NE.0) GOTO 598
  1383.  
  1384. C- Stockage de la matrice
  1385.           nrigel=irigel(/2) +1
  1386.           segadj,mrigid
  1387. C           jrige = NRIGE0 + irige
  1388.           jrige=nrigel
  1389.           COERIG(jrige)   = 1.
  1390.           IRIGEL(1,jrige) = ipmail
  1391.           IRIGEL(2,jrige) = 0
  1392.           IRIGEL(3,jrige) = ipdesc
  1393.           IRIGEL(4,jrige) = ipmatr
  1394.           IRIGEL(5,jrige) = NIFOUR
  1395.           IRIGEL(6,jrige) = 0
  1396.           IRIGEL(7,jrige) = 0
  1397.           IRIGEL(8,jrige) = 0
  1398.  
  1399.         ENDDO
  1400. C- Fin de la boucle sur les partitions
  1401. C
  1402.  598    CONTINUE
  1403.         IF (MOSTRS.NE.0) THEN
  1404.           nomid = MOSTRS
  1405.           IF (lsupco) SEGSUP,nomid
  1406.           notype = MOTYPS
  1407.           SEGSUP,notype
  1408.         ENDIF
  1409.         IF (MOCARA.NE.0) THEN
  1410.           NOMID = MOCARA
  1411.           SEGSUP,NOMID
  1412.           notype = MOTYPC
  1413.           SEGSUP,notype
  1414.         ENDIF
  1415. C
  1416.         NOMID=MODEPL
  1417.         IF (lsupde) SEGSUP,NOMID
  1418.         NOMID = MOFORC
  1419.         IF (lsupfo) SEGSUP,NOMID
  1420.  
  1421.  599    CONTINUE
  1422.  
  1423.         IF (IERR.NE.0) GOTO 999
  1424. C
  1425.  500  CONTINUE
  1426. C* Fin de la boucle sur les modeles elementaires
  1427.  
  1428.  999  CONTINUE
  1429.       IF (IERR.NE.0) THEN
  1430.         ktrace = -1
  1431. **  En situation d'erreur, on laisse le menage faire son travail
  1432. **      CALL DERIGI(MRIGID,ktrace,msorse)
  1433. **      SEGSUP,MRIGID
  1434.         IPRIGG = 0
  1435.       ELSE
  1436.         if(irigel(/2).eq.0) then
  1437.           call erreur (86)
  1438.           return
  1439.         endif
  1440.         SEGDES,MRIGID
  1441.         IPRIGG = MRIGID
  1442.       ENDIF
  1443.  
  1444. c      RETURN
  1445.       END
  1446.  
  1447.  
  1448.  
  1449.  
  1450.  
  1451.  

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