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

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