Télécharger la plaquette Cast3M


fr en

Télécharger rela1.eso

Retour à la liste

Numérotation des lignes : 

rela1

  1. C RELA1     SOURCE    MB234859  26/09/03    21:15:08     12636          
  2.       SUBROUTINE RELA1
  3. C
  4. C     RELAZIONE LINEARE TRA DDL
  5. C
  6. C     ICONR  NUMERO OGGETTI MELEME DELLA RELAZIONE
  7. C     ITYESR(NBSOUS) LISTA TIPI ELEMENTI CONTENUTI IN MELEME-N
  8. C     MUNESR(NBSOUS)  NUMERO ELEMENTI IN OGNI SOTTOSTRUTTURA DI MELEME-N
  9. C     LNODSR(NNR) LISTA NODI MELEME-N
  10. C     IPOR1(N)  PUNTATORE SU MWREL1
  11. C     IPOR2(N)  PUNTATORE SU MWREL2
  12. C     IPOR3(N)  PUNTATORE SU MWREL3
  13. C     COEFR(N)  COEFFICIENTI RELAZIONE
  14. C     INCREL(N) IDENTIFICATORE NOME DDL
  15. C
  16.       IMPLICIT INTEGER(I-N)
  17.       IMPLICIT REAL*8(A-H,O-Z)
  18.  
  19. -INC SMELEME
  20. -INC SMCOORD
  21. -INC SMRIGID
  22. -INC TMTRAV
  23. -INC SMCHPOI
  24. -INC SMLMOTS
  25. -INC SMLREEL
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCREEL
  29. -INC CCGEOME
  30. -INC CCHAMP
  31.  
  32.       SEGMENT /MWGGM1/(IPORE1(0))
  33.       SEGMENT /MWGGM2/(IPORE2(0))
  34.       SEGMENT /MWGGM3/(IPORE3(0))
  35.       SEGMENT /MWGGM4/(INCREL(0))
  36.       SEGMENT /MWGGM5/(COEFR(0)*D)
  37.       SEGMENT /MWREL1/(ITYESR(0)),IEL11.MWREL1,
  38.      +     IEL21.MWREL1,iel01.MWREL1
  39.       SEGMENT /MWREL2/(MUNESR(0)),IEL12.MWREL2,
  40.      +     IEL22.MWREL2,iel02.MWREL2
  41.       SEGMENT /MWREL3/(LNODSR(0)),IEL13.MWREL3,
  42.      +     IEL23.MWREL3
  43.       CHARACTER*4 MOTPV(3)
  44.       CHARACTER*4 MOCORI(7)
  45.       CHARACTER*4 MOTPM(2)
  46.       CHARACTER*4 MOTDDL(3)
  47.       CHARACTER*4 MOTBLO(3)
  48.       CHARACTER*4 MODEPL(8)
  49.       CHARACTER*4 MOROTA(5)
  50.       CHARACTER*4 MODUAL(1),MONORM(1)
  51.       CHARACTER*8 MILMOT
  52.       DIMENSION XNOR(3)
  53.       DATA MOTPV /'MINI','MAXI','FROT'/
  54.       DATA MOTPM /'+   ','-   '/
  55.       DATA MOCORI/'CORI','ENSE','ACCR','GLIS','BARY','MILI','TUYA'/
  56. C
  57.       DATA MOTBLO/'DEPL','ROTA','DIRE'/
  58.       DATA MODEPL/'UX  ','UY  ','UZ  ','UR  ','UZ  ','UT  ',
  59.      &     'ALFA','BETA'/
  60.       DATA MOROTA/'RX  ','RY  ','RZ  ','RT  ','RS  '/
  61.       DATA MODUAL/'DUAL'/
  62.       DATA MONORM/'NORM'/
  63. c
  64.       iel01=0
  65.       iel02=0
  66. C
  67. C     EST-CE UNE RELATION DE CORPS RIGIDE ?
  68. C
  69.  
  70.       CALL LIRMOT(MOCORI,7,ICORI,0)
  71. *      write(ioimp,*) 'rela1, icori=',icori
  72.       IF (ICORI.NE.0) THEN
  73. c     BARY
  74.          if (icori.eq.5) then
  75.             call relaba
  76. c     MILI
  77.          ELSE IF (ICORI.EQ.6) THEN
  78.             CALL LIROBJ('LISTMOTS',IP0,0,IRETOU)
  79.             IF (IERR.NE.0) RETURN
  80.             CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  81.             IF (IERR.NE.0) RETURN
  82.             CALL RELAMI(IP0,IPT1,IPRIG)
  83.             CALL ECROBJ('RIGIDITE',IPRIG)
  84. c     TUYA
  85.          elseif(icori.eq.7) then
  86.             call reltuy
  87. c     CORI, ENSE, ACCRO, GLIS
  88.          else
  89.             CALL RELA2(MOCORI(ICORI))
  90.          endif
  91.          RETURN
  92.       ENDIF
  93. C
  94. C     EST-CE UNE CONDITION UNILATERALE ?
  95. C
  96.       NILATE=0
  97.       CALL LIRMOT (MOTPV,3,IPO,0)
  98.       IF(IPO.EQ.1) NILATE=-1
  99.       IF(IPO.EQ.2) NILATE=+1
  100.       IF(IPO.EQ.3) NILATE=+2
  101. C
  102. C     Lecture eventuelle mot cle NORM
  103.       IRENOR=0
  104.       CALL LIRMOT(MONORM,1,IRENOR,0)
  105. C
  106. C     INITIALISATIONS
  107. C
  108.       COEX=1.D0
  109. C
  110. C     Deformations planes ou contraintes planes ou defo. plane gene :
  111.       IF (IFOUR.EQ.-1.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN
  112.          LDEPL=2
  113.          IADEPL=0
  114.          LROTA=1
  115.          IAROTA=2
  116. C     Axisymetrique :
  117.       ELSE IF (IFOUR.EQ.0) THEN
  118.          LDEPL=2
  119.          IADEPL=3
  120.          LROTA=1
  121.          IAROTA=3
  122. C     Fourier :
  123.       ELSE IF (IFOUR.EQ.1) THEN
  124.          LDEPL=3
  125.          IADEPL=3
  126.          LROTA=2
  127.          IAROTA=2
  128. C     Tridimensionnel :
  129.       ELSE IF (IFOUR.EQ.2) THEN
  130.          LDEPL=3
  131.          IADEPL=0
  132.          LROTA=3
  133.          IAROTA=0
  134. C     Massif 1D (IDIM=1) :
  135.       ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN
  136.          LDEPL=1
  137.          IADEPL=0
  138.          IF (IFOUR.GE.12) IADEPL=3
  139.          LROTA=0
  140.          IAROTA=0
  141. C     Autres cas :
  142.       ELSE
  143.          LDEPL=0
  144.          IADEPL=0
  145.          LROTA=0
  146.          IAROTA=0
  147.       ENDIF
  148. *
  149. *     LECTURE EVENTUELLE D'UN MODELE
  150. *
  151.       CALL LIROBJ('MMODEL',IPOMOD,0,IREMOD)
  152.       IF(IREMOD.NE.0) THEN
  153.          CALL RELMOD(IPOMOD)
  154.          RETURN
  155.       ENDIF
  156. *
  157. *     LECTURE EVENTUELLE D'UN CHPOINT
  158. *
  159.       CALL LIROBJ('CHPOINT',IPOCHP,0,IRECHP)
  160.       IF (IRECHP.EQ.1) GOTO 500
  161. C
  162.       SEGINI MWGGM1
  163.       SEGINI MWGGM2
  164.       SEGINI MWGGM3
  165.       SEGINI MWGGM4
  166.       SEGINI MWGGM5
  167. C
  168.  200  CONTINUE
  169. C
  170. C     LETTURA OPERANDI
  171. C
  172.       SEGINI MWREL1
  173.       SEGINI MWREL2
  174.       SEGINI MWREL3
  175.       IPORE1(**)=MWREL1
  176.       IPORE2(**)=MWREL2
  177.       IPORE3(**)=MWREL3
  178.  1320 CONTINUE
  179.       CALL LIRREE(XR,0,IRETOF)
  180.       IF(IRETOF.EQ.0) THEN
  181.         CALL LIRENT(IR,0,IRETOI)
  182.         XR=1.D0
  183.         IF(IRETOI.NE.0) XR=DBLE(IR)
  184.       ENDIF
  185.       CALL LIRMOT(NOMDD,LNOMDD,IPO,0)
  186.       IF(IPO.NE.0) THEN
  187.         COEFR(**)=XR*COEX
  188.         INCREL(**)=IPO
  189.       ELSE
  190. *
  191. *     ON REGARDE SI IL Y A DEPL OU ROTA SUIVI DE DIRECTION
  192. *
  193. C     En DIMENSION 1, le mot-cle 'DIRE' est interdit.
  194.          IF (IDIM.EQ.1) THEN
  195.             INTERR(1)=IDIM
  196.             MOTERR(1:4)=MOTBLO(3)
  197.             CALL ERREUR(971)
  198.             GOTO 559
  199.          ENDIF
  200.          CALL LIRMOT(MOTBLO,2,IMOT,1)
  201.          IF (IMOT.EQ.0) GOTO 559
  202.          IF (IMOT.EQ.1) THEN
  203.             IBDDL=LDEPL
  204.             DO 4481 IA=1,IBDDL
  205.                MOTDDL(IA)=MODEPL(IADEPL+IA)
  206.  4481       CONTINUE
  207.          ENDIF
  208.          IF (IMOT.EQ.2) THEN
  209.             IBDDL=LROTA
  210.             DO 4482 IA=1,IBDDL
  211.                MOTDDL(IA)=MOROTA(IAROTA+IA)
  212.  4482       CONTINUE
  213.          ENDIF
  214.          CALL LIRMOT(MOTBLO(3),1,IMOT,1)
  215.          IF(IMOT.EQ.0) GO TO 559
  216.          CALL LIROBJ('POINT',KPOINT,1,IRETOU)
  217.          IF(IRETOU.EQ.0) GO TO 559
  218.          YL=0.D0
  219.          XPREC=SQRT(XPETIT/XZPREC)
  220.          segact mcoord
  221.          DO 4484 IA=1,IDIM
  222.             XVAL    =XCOOR((KPOINT-1)*(IDIM+1)+IA)
  223.             XNOR(IA)=XVAL
  224.             IF(ABS(XVAL) .GT. XPREC)THEN
  225.               YL=YL+XVAL**2
  226.             ENDIF
  227.  4484    CONTINUE
  228.          IF (YL .EQ. 0.D0) THEN
  229.             CALL ERREUR(239)
  230.             GOTO 559
  231.          ENDIF
  232.          YL=1.D0/SQRT(YL)
  233.          DO 4485 IA=1,IDIM
  234.             XNOR(IA)=XNOR(IA)*YL
  235.  4485    CONTINUE
  236. *
  237. *     ON LIT LE MAILLAGE
  238. *
  239.          CALL LIROBJ('POINT   ',MILPOI,0,IRETOU)
  240.          IF(IRETOU.NE.0) THEN
  241.             MILMOT='POINT   '
  242.          ELSE
  243.             CALL LIROBJ('MAILLAGE',MILPOI,1,IRETOU)
  244.             IF(IRETOU.EQ.0)  GO TO 559
  245.             MILMOT='MAILLAGE'
  246.          ENDIF
  247. *
  248. *     PUIS ON REMET DES OBJETS DANS LA PILE
  249. *
  250.          DO 4477 IB=1,IBDDL
  251.             IF(IB.GE.2.AND.COEX.EQ. 1.D0) CALL ECRCHA(MOTPM(1))
  252.             IF(IB.GE.2.AND.COEX.EQ.-1.D0) CALL ECRCHA(MOTPM(2))
  253.             CALL ECROBJ(MILMOT,MILPOI)
  254.             CALL ECRCHA(MOTDDL(IB))
  255.             IF(IBDDL.EQ.1) THEN
  256.                XRCOO=XR
  257.             ELSE
  258.                XRCOO=XR*XNOR(IB)
  259.             ENDIF
  260.             CALL ECRREE(XRCOO)
  261.  4477    CONTINUE
  262.          GO TO 1320
  263.       ENDIF
  264. C
  265.       CALL LIROBJ('POINT   ',KPOINT,0,IRETOU)
  266.       IF(IRETOU.NE.0) GO TO 110
  267.       CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU)
  268.       IF(IRETOU.EQ.0)  GO TO 559
  269.       MELEME=KOBJET
  270.       SEGACT MELEME
  271.       NBSOUS=LISOUS(/1)
  272.       IF(NBSOUS.EQ.0) GO TO 120
  273. C     OBJET COMPLEXE
  274.       NNR=1
  275.       DO 130 IS=1,NBSOUS
  276.          IPT1=LISOUS(IS)
  277.          SEGACT IPT1
  278.          ITYESR(**)=IPT1.ITYPEL
  279.          NBNN      =IPT1.NUM(/1)
  280.          NBELEM    =IPT1.NUM(/2)
  281.          MUNESR(**)=IPT1.NUM(/2)
  282.          IF(NNR.EQ.1)LNODSR(**)=IPT1.NUM(1,1)
  283.          DO 140  I1=1,NBELEM
  284.             DO 1401  I2=1,NBNN
  285.                DO 150  I3=1,NNR
  286.                   IF(LNODSR(I3).EQ.IPT1.NUM(I2,I1)) GO TO 1401
  287.  150           CONTINUE
  288.                NNR=NNR+1
  289.                LNODSR(**)=IPT1.NUM(I2,I1)
  290.  1401       CONTINUE
  291.  140     CONTINUE
  292.          SEGDES IPT1
  293.  130  CONTINUE
  294.       SEGDES MWREL1
  295.       SEGDES MWREL2
  296.       SEGDES MWREL3
  297.       SEGDES MELEME
  298.       GO TO 160
  299. C
  300. C     OBJET SIMPLE
  301. C
  302.  120  CONTINUE
  303.       ITYESR(**)=ITYPEL
  304.       NBNN      =NUM(/1)
  305.       NBELEM    =NUM(/2)
  306.       MUNESR(**)=NUM(/2)
  307.       NNR=0
  308.       DO 170  I1=1,NBELEM
  309.          DO 1701  I2=1,NBNN
  310.             DO 180  I3=1,NNR
  311.                IF(LNODSR(I3).EQ.NUM(I2,I1)) GO TO 1701
  312.  180        CONTINUE
  313.             NNR=NNR+1
  314.             LNODSR(**)=NUM(I2,I1)
  315.  1701    CONTINUE
  316.  170  CONTINUE
  317.       SEGDES MWREL1
  318.       SEGDES MWREL2
  319.       SEGDES MWREL3
  320.       SEGDES MELEME
  321.       GO TO 160
  322.  110  CONTINUE
  323. C
  324. C     OBJET POINT
  325. C
  326.       ITYESR(**)=1
  327.       MUNESR(**)=1
  328.       LNODSR(**)=KPOINT
  329.       SEGDES MWREL1
  330.       SEGDES MWREL2
  331.       SEGDES MWREL3
  332. C     FINE OPERANDO RELAZIONE
  333.  160  CONTINUE
  334.       ICONR=IPORE1(/1)+1
  335. C     LETTURA + O -
  336.       IF(ICONR.EQ.2) CALL LIRMOT(MOTPM,2,IPO,0)
  337.       IF (IPO.EQ.0)  GO TO 300
  338. *     LIRMOT(MOTPM,2,IPO,1) goto 559
  339.       IF(ICONR.GT.2) CALL LIRMOT(MOTPM,2,IPO,0)
  340.       IF (IPO.EQ.0)  GO TO 300
  341.       COEX=1.D0
  342.       IF (IPO.EQ.2) COEX=-1.D0
  343. C     SI CERCA DI LEGGERE UN NUOVO OPERANDO DELLA RELAZIONE
  344.       GO TO  200
  345. C
  346. C     VERIFICA CONGRUENZA OPERANDI
  347. C
  348.  300  CONTINUE
  349. C
  350.       ICONR=IPORE1(/1)
  351.       IEL11=IPORE1(1)
  352.       IEL12=IPORE2(1)
  353.       IEL13=IPORE3(1)
  354. *     on autorise maintenant le point unique qui va etre applique sur
  355. *     chaque relation
  356. *     old      ityes=0
  357. *     old      munes=0
  358.       nsor=0
  359.       nnr=0
  360.       do 305 io=1,iconr
  361. *     write(ioimp,*) 'io=',io
  362.          iel11=ipore1(io)
  363.          iel12=ipore2(io)
  364.          iel13=ipore3(io)
  365.          segact iel11,iel12,iel13
  366.          nsor1=iel11.ityesr(/1)
  367.          nsor=max(nsor,nsor1)
  368.          NNR1=IEL13.LNODSR(/1)
  369.          nnr=max(nnr,nnr1)
  370. *     old         do 315 i1=1,nsor1
  371.          if (io.eq.1) then
  372.             segini,iel01=iel11
  373.             segini,iel02=iel12
  374.          else
  375.             do 315 i1=1,nsor1
  376.                ityes=iel01.ityesr(i1)
  377.                iel01.ityesr(i1)=max(iel11.ityesr(i1),ityes)
  378. *     write(ioimp,*) 'i1=',i1,' ityes=',ityes
  379. *     $              ,' iel11.ityesr(i1)',iel11.ityesr(i1)
  380.                munes=iel02.munesr(i1)
  381.                iel02.munesr(i1)=max(iel12.munesr(i1),munes)
  382. *     write(ioimp,*) 'i1=',i1,' munes=',munes
  383. *     $              ,' iel12.munesr(i1)',iel12.munesr(i1)
  384.  315        continue
  385.          endif
  386.  305  continue
  387.       DO 310 IO=1,ICONR
  388. *     write(ioimp,*) 'io=',io
  389.          IEL21=IPORE1(IO)
  390.          IEL22=IPORE2(IO)
  391.          IEL23=IPORE3(IO)
  392.          NSOR2=IEL21.ITYESR(/1)
  393. *     write(ioimp,*) 'nsor,nsor2=',nsor,nsor2
  394.          IF(NSOR.NE.NSOR2.and.nsor2.ne.1) GO TO 556
  395.          DO 320 I2=1,NSOR2
  396.             ityes=iel01.ityesr(i2)
  397.             munes=iel02.munesr(i2)
  398. *     write(ioimp,*) 'i2=',i2,' ityes=',ityes,'munes=',munes
  399. *     write(ioimp,*) 'iel21.ityesr=',iel21.ityesr(i2)
  400. *     $           ,' IEL22.MUNESR(I2)=',IEL22.MUNESR(I2)
  401.             if (iel21.ityesr(i2).eq.1.and.iel22.munesr(i2).eq.1.and.
  402.      >           nsor2.eq.1) goto 320
  403.             IF(IEL21.ITYESR(I2).NE.ityes) GO TO 556
  404.             IF(IEL22.MUNESR(I2).NE.munes) GO TO 556
  405.  320     CONTINUE
  406.          NNR2=IEL23.LNODSR(/1)
  407. *     write(ioimp,*) 'nnr,nnr2=',nnr,nnr2
  408.          IF(NNR2.NE.NNR.and.nnr2.ne.1) GO TO 556
  409.  310  CONTINUE
  410.       if (iel01.ne.0) SEGSUP,IEL01
  411.       if (iel02.ne.0) SEGSUP,IEL02
  412.       NBRELA=NNR
  413. C     OPERANDI  CONGRUENTI  FINE VERIFICA
  414. C
  415. C     CARICAMENTO MATRICE  RIGIDEZZA
  416. C
  417.       NRIGEL=1
  418.       SEGINI MRIGID
  419.       ICHOLE=0
  420.       IMGEO1=0
  421.       IMGEO2=0
  422.       ISUPEQ=0
  423.       IFORIG=IFOUR
  424.       COERIG(1)=1.D0
  425.       MTYMAT='RIGIDITE'
  426.       KRIGI=MRIGID
  427. C
  428. C     ON INITIALISE LE SEGMENT MELEME ASSOCIE AUX BLOCAGES
  429. C
  430.       SEGACT MCOORD*MOD
  431.       NBPTSO=nbpts
  432. *      write(ioimp,*) 'nbptso=',nbptso
  433.       NBPTS=NBPTSO+NBRELA
  434.       SEGADJ MCOORD
  435.       NBSOUS=0
  436.       NBREF=0
  437.       NBNN=1+ICONR
  438.       NBELEM=NBRELA
  439.       SEGINI IPT2
  440.       IRIGEL(1,1)=IPT2
  441.       IPT2.ITYPEL=22
  442. C     SI CREANO UNO PUNTO PER OGNI RELAZIONE
  443.       DO 400 I4=1,NBRELA
  444.          IPT2.ICOLOR(I4)=IDCOUL
  445.          IPTS=NBPTSO+I4
  446.          JPTS=(IPTS-1)*(IDIM+1)
  447.          DO 410 iidim=1,idim+1
  448.             XCOOR(JPTS+iidim)=0.D0
  449.  410     CONTINUE
  450. C     PUNTI ASSOCIATI AI MOLTIPLICATORI
  451.          IPT2.NUM(1,I4)=IPTS
  452.  400  CONTINUE
  453. C     CARICAMENTO NODI PSEUDO-ELEMENTI
  454.       DO 420 I8=1,ICONR
  455.          MWREL3=IPORE3(I8)
  456.          I10=I8+1
  457.          SEGACT MWREL3
  458.          LNOMAX=LNODSR(/1)
  459.          DO 430 I9=1,NBRELA
  460.             NN=LNODSR(min(I9,lnomax))
  461.             NPN=(NN-1)*(IDIM+1)
  462.             IPT2.NUM(I10,I9)=NN
  463.             NPL1=(IPT2.NUM(1,I9)-1)*(IDIM+1)
  464.             do 4301 iidim=1,idim+1
  465.                XCOOR(NPL1+iidim)=XCOOR(NPL1+iidim)+XCOOR(NPN+iidim)
  466.  4301       continue
  467.  430     CONTINUE
  468.          SEGDES MWREL3
  469.  420  CONTINUE
  470. C
  471. C     COORDINATE DEI BARICENTRI ASSOCIATI ALLE RELAZIONI
  472. C
  473. *      write(ioimp,*) 'iconr,nbrela,idim=',iconr,nbrela,idim
  474.       DO 425 I4=1,NBRELA
  475.          NPL1=(IPT2.NUM(1,I4)-1)*(IDIM+1)
  476.          do 4251 iidim=1,idim+1
  477.             XCOOR(NPL1+iidim)=XCOOR(NPL1+iidim)/ICONR
  478.  4251    continue
  479.  425  CONTINUE
  480.  
  481.       SEGDES IPT2
  482.       IRIGEL(2,1)=0
  483.       IRIGEL(5,1)=NIFOUR
  484.       IRIGEL(6,1)=NILATE
  485.       NLIGRP=ICONR+1
  486.       NLIGRD=NLIGRP
  487.       SEGINI DESCR
  488.       IRIGEL(3,1)=DESCR
  489.       LISINC(1)='LX'
  490.       LISDUA(1)='FLX'
  491.       NOELEP(1)=1
  492.       NOELED(1)=1
  493.       DO 700 I1=1,ICONR
  494.          I2=I1+1
  495.          I3=ABS(INCREL(I1))
  496.          LISINC(I2)=NOMDD(I3)
  497.          LISDUA(I2)=NOMDU(I3)
  498.          NOELEP(I2)=I2
  499.          NOELED(I2)=I2
  500.  700  CONTINUE
  501.       SEGDES DESCR
  502.       NELRIG=NBRELA
  503.       rigrel=0
  504.       SEGINI xMATRI
  505.       IRIGEL(4,1)=xMATRI
  506. *     LVAL=(NBNN*NBNN+NBNN)/2
  507.       NLIGRP=NBNN
  508.       NLIGRD=NBNN
  509. *     SEGINI XMATRI
  510.       DO 740 I6=1,NELRIG
  511. *     IMATTT(I1)=XMATRI
  512. *     740   CONTINUE
  513.          RE(1,1,i6)= 0.D0
  514.          I3=3
  515.          DO 760 I1=2,NBNN
  516.             I4=I1-1
  517.             I2=1
  518.             RE(I1,I2,i6)=COEFR(I4)
  519.             RE(I2,I1,i6)=COEFR(I4)
  520.  760     CONTINUE
  521.  740  continue
  522.       SEGDES XMATRI
  523. *     SEGDES IMATRI
  524.       call relasi(mrigid)
  525.       SEGDES MRIGID
  526.       CALL ECROBJ('RIGIDITE',KRIGI)
  527.  559  CONTINUE
  528.       ICONR=IPORE1(/1)
  529.       DO 558 I1=1,ICONR
  530.          MWREL1=IPORE1(I1)
  531.          MWREL2=IPORE2(I1)
  532.          MWREL3=IPORE3(I1)
  533.          SEGSUP MWREL1
  534.          SEGSUP MWREL2
  535.          SEGSUP MWREL3
  536.  558  CONTINUE
  537.       SEGSUP MWGGM1
  538.       SEGSUP MWGGM2
  539.       SEGSUP MWGGM3
  540.       SEGSUP MWGGM4
  541.       SEGSUP MWGGM5
  542.       RETURN
  543.  556  CONTINUE
  544. C     SEGDES IEL11,IEL12,IEL13,IEL21,IEL22,IEL23
  545.       CALL ERREUR(324)
  546.       GO TO 559
  547. *
  548. *  on arrive en 500 avec la syntaxe chp1 DUAL chp2
  549. *
  550.  500  CONTINUE
  551.       CALL ACTOBJ('CHPOINT',IPOCHP,1)
  552. C
  553. C     Lecture eventuelle mot cle DUAL
  554.       IRECHD=0
  555.       CALL LIRMOT(MODUAL,1,IRECHD,0)
  556.       IF (IRECHD.EQ.1) THEN
  557.         CALL LIROBJ('CHPOINT',IPOCHD,1,IRECHD)
  558.         IF (IERR.NE.0) RETURN
  559. C
  560. C       Deuxieme tentative de lecture du mot cle NORM, bloque par DUAL
  561.         IF (IRENOR.EQ.0) CALL LIRMOT(MONORM,1,IRENOR,0)
  562.       ENDIF
  563. C
  564. C     ipochp est le chpoin de la relation
  565. C     ipochd est le chpoin du vecteur de controle
  566.       mchpoi=ipochp
  567.       segact mchpoi
  568. *  reserver la place pour le lx
  569.       nligrp=1
  570.       nbnp=0
  571.       nbnd=0
  572.       do ims=1,ipchp(/1)
  573.         msoupo=ipchp(ims)
  574.         segact msoupo
  575.         mpoval=ipoval
  576.         segact mpoval
  577.         nligrp=nligrp+vpocha(/1)*vpocha(/2)
  578.         meleme=igeoc
  579.         segact meleme
  580.         nbnp=nbnp+num(/2)
  581.       enddo
  582. *     write(6,*) ' nligrp nbelep ',nligrp,nbelep
  583.       IF (IRECHD.EQ.1) THEN
  584.         mchpoi=ipochd
  585.         segact mchpoi
  586.         xnormd = 0.d0
  587. *  reserver la place pour le flx
  588.         nligrd=1
  589.         do ims=1,ipchp(/1)
  590.           msoupo=ipchp(ims)
  591.           segact msoupo
  592.           mpoval=ipoval
  593.           segact mpoval
  594.           nligrd=nligrd+vpocha(/1)*vpocha(/2)
  595.           if (IRENOR.EQ.0) THEN
  596.             do j=1,vpocha(/2)
  597.               do i=1,vpocha(/1)
  598.                 xnormd=max(xnormd,abs(vpocha(i,j)))
  599.               enddo
  600.             enddo
  601.             if (xnormd.lt.1d-5.or.xnormd.gt.1E5) then
  602.               reaerr(1)=xnormd
  603.               call erreur(1118)
  604.               return
  605.             endif
  606.           endif
  607.           meleme=igeoc
  608.           segact meleme
  609.           nbnd=nbnd+num(/2)
  610.         enddo
  611.         nligrp=nligrp+nligrd-1
  612.       ENDIF
  613. C
  614. C     Creation rigidite
  615.       NLIGRD=NLIGRP
  616.       SEGINI,DESCR
  617. C
  618.       NBSOUS=0
  619.       NBREF=0
  620.       NBELEM=1
  621.       NBNN=NBNP+NBND+1
  622.       SEGINI,MELEME
  623.       ITYPEL=22
  624. C
  625.       NELRIG=1
  626.       RIGREL=0
  627.       SEGINI,XMATRI
  628.       SYMRE=0
  629.       IF (IRECHD.EQ.1) SYMRE=2
  630. C
  631.       NRIGEL=1
  632.       SEGINI,MRIGID
  633.       ICHOLE=0
  634.       IMGEO1=0
  635.       IMGEO2=0
  636.       ISUPEQ=0
  637.       IFORIG=IFOUR
  638.       COERIG(1)=1.D0
  639.       MTYMAT='RIGIDITE'
  640.       irigel(1,1)=meleme
  641.       irigel(3,1)=descr
  642.       irigel(4,1)=xmatri
  643.       IRIGEL(5,1)=NIFOUR
  644.       irigel(6,1)=nilate
  645.       irigel(7,1)=symre
  646. *  remplissage maillage descripteur et valeur
  647. *  le premier noeud sera fait a la fin. C'est le multiplicateur de lagrange
  648.       iel=1
  649.       ire=1
  650.       nbnp=1
  651.       mchpoi=ipochp
  652.       do ims=1,ipchp(/1)
  653.         msoupo=ipchp(ims)
  654.         ipt1=igeoc
  655.         mpoval=ipoval
  656.         do i=1,ipt1.num(/2)
  657.           iel=iel+1
  658.           num(iel,1)=ipt1.num(1,i)
  659.           nbnp=nbnp+1
  660.           do j=1,nocomp(/2)
  661.             call place(nomdd,lnomdd,ipo,nocomp(j))
  662.             if (ipo.eq.0) then
  663.               if (irechd.EQ.0) THEN
  664.                 call place(nomdu,lnomdu,ipo,nocomp(j))
  665.               endif
  666.               if (ipo.eq.0) then
  667.                 moterr(1:4)=nocomp(j)
  668.                 moterr(5:11)='PRIMALE'
  669.                 call erreur(1117)
  670.               endif
  671.             endif
  672.             if (nocomp(j).eq.'LX  ') call erreur(1125)
  673.             IF (IERR.NE.0) RETURN
  674. C
  675.             ire=ire+1
  676.             lisinc(ire)=nomdd(ipo)
  677.             lisdua(ire)=nomdu(ipo)
  678.             noelep(ire)=iel
  679.             noeled(ire)=iel
  680.             re(1,ire,1)=vpocha(i,j)
  681.             if (irechd.ne.1) re(ire,1,1)=re(1,ire,1)
  682. *           write(6,*) ' 1 ire inc dua ',lisinc(ire),lisdua(ire)
  683.           enddo
  684.         enddo
  685.       enddo
  686.  
  687.       IF (IRECHD.EQ.1) THEN
  688.         nbnd=1
  689.         mchpoi=ipochd
  690.         do ims=1,ipchp(/1)
  691.           msoupo=ipchp(ims)
  692.           ipt1=igeoc
  693.           mpoval=ipoval
  694.           do i=1,ipt1.num(/2)
  695.             iel=iel+1
  696.             num(iel,1)=ipt1.num(1,i)
  697.             nbnd=nbnd+1
  698.             do j=1,nocomp(/2)
  699.               call place(nomdu,lnomdu,ipo,nocomp(j))
  700.               if(ipo.eq.0) then
  701.                 moterr(1:4)=nocomp(j)
  702.                 moterr(5:10)='DUALE'
  703.                 call erreur(1117)
  704.               endif
  705.               if (lisinc(ire).eq.'LX  ') call erreur(1125)
  706.               IF (IERR.NE.0) RETURN
  707. C
  708.               ire=ire+1
  709.               lisinc(ire)=nomdd(ipo)
  710.               lisdua(ire)=nocomp(j)
  711.               noelep(ire)=iel
  712.               noeled(ire)=iel
  713.               re(ire,1,1)= -vpocha(i,j)
  714. *             write(6,*) ' 2 ire inc dua ',lisinc(ire),lisdua(ire)
  715.             enddo
  716.           enddo
  717.         enddo
  718.       ENDIF
  719. C
  720. C     multiplicateur de lagrange
  721.       segact MCOORD*MOD
  722.       nbpts=nbpts+1
  723.       segadj mcoord
  724.       xl=0.d0
  725.       yl=0.d0
  726.       zl=0.d0
  727.       dl=0.d0
  728.       do i=2,num(/1)
  729.        ip=num(i,1)-1
  730.        xl=xl+xcoor(ip*(idim+1)+1)
  731.        if (idim.gt.1) yl=yl+xcoor(ip*(idim+1)+2)
  732.        if (idim.gt.2) zl=zl+xcoor(ip*(idim+1)+3)
  733.        dl=dl+xcoor(ip*(idim+1)+idim+1)
  734.       enddo
  735.       nbnnr=num(/1)-1
  736.       xl=xl/nbnnr
  737.       yl=yl/nbnnr
  738.       zl=zl/nbnnr
  739.       dl=dl/nbnnr
  740.       xcoor((nbpts-1)*(idim+1)+1)=xl
  741.       if (idim.gt.1) xcoor((nbpts-1)*(idim+1)+2)=yl
  742.       if (idim.gt.2) xcoor((nbpts-1)*(idim+1)+3)=zl
  743.       xcoor((nbpts-1)*(idim+1)+idim+1)=dl
  744.       lisinc(1)= 'LX'
  745.       lisdua(1)='FLX'
  746.       num(1,1)=nbpts
  747.       noelep(1)=1
  748.       noeled(1)=1
  749.       re(1,1,1)=0.D0
  750. C
  751. C     Mot-cle NORM
  752.       IF (IRENOR.EQ.1) THEN
  753.         JGN=LOCHPO
  754.         JGM=NLIGRP-1
  755.         SEGINI,MLMOT1
  756.         NBPRIM=0
  757.         DO 14 IB0=2,NLIGRP
  758.           DO 13 IB1=1,JGM
  759.             IF (LISINC(IB0).EQ.MLMOT1.MOTS(IB1)) GOTO 13
  760.             NBPRIM=NBPRIM+1
  761.             MLMOT1.MOTS(nbprim)=LISINC(IB0)
  762.             GOTO 14
  763.   13      CONTINUE
  764.   14    CONTINUE
  765.         JG=NBPRIM
  766.         SEGINI,MLREE1
  767.         DO IB1=2,NLIGRP
  768.           DO ICP=1,NBPRIM
  769.             IF (LISINC(IB1).EQ.MLMOT1.MOTS(ICP)) GOTO 15
  770.           ENDDO
  771.           CALL ERREUR(5)
  772.   15      CONTINUE
  773.           MLREE1.PROG(ICP)=MLREE1.PROG(ICP)+RE(1,IB1,1)
  774.         ENDDO
  775.         ZNOR=0.D0
  776.         DO IB2=1,NBPRIM
  777.           YNOR=MLREE1.PROG(IB2)
  778.           ZNOR=ZNOR+(YNOR*YNOR)
  779.         ENDDO
  780.         IF (ZNOR.EQ.0.D0) ZNOR=1.D0
  781.         ZNOR=1.D0/SQRT(ZNOR)
  782.         DO IB3=2,NLIGRP
  783.           RE(1,IB3,1)=RE(1,IB3,1)*ZNOR
  784.           IF (IRECHD.EQ.0) RE(IB3,1,1)=RE(IB3,1,1)*ZNOR
  785.         ENDDO
  786.         SEGSUP,MLREE1,MLMOT1
  787. C
  788.         IF (IRECHD.EQ.1) THEN
  789.           JGN=LOCHPO
  790.           JGM=NLIGRD-1
  791.           SEGINI,MLMOT2
  792.           NBDUAL=0
  793.           DO 17 IB0=2,NLIGRD
  794.             DO 16 IB1=1,JGM
  795.               IF (LISDUA(IB0).EQ.MLMOT2.MOTS(IB1)) GOTO 16
  796.               NBDUAL=NBDUAL+1
  797.               MLMOT2.MOTS(nbdual)=LISDUA(IB0)
  798.               GOTO 17
  799.   16        CONTINUE
  800.   17      CONTINUE
  801.           JG=NBDUAL
  802.           SEGINI,MLREE2
  803.           DO IB1=2,NLIGRD
  804.             DO ICP=1,NBDUAL
  805.               IF (LISDUA(IB1).EQ.MLMOT2.MOTS(ICP)) GOTO 18
  806.             ENDDO
  807.             CALL ERREUR(5)
  808.   18        CONTINUE
  809.             MLREE2.PROG(ICP)=MLREE2.PROG(ICP)+RE(IB1,1,1)
  810.           ENDDO
  811.           ZNOR=0.D0
  812.           DO IB2=1,nbdual
  813.             YNOR=MLREE2.PROG(IB2)
  814.             ZNOR=ZNOR+(YNOR*YNOR)
  815.           ENDDO
  816.           IF (ZNOR.EQ.0.D0) ZNOR=1.D0
  817.           ZNOR=1.D0/SQRT(ZNOR)
  818.           DO IB3=2,NLIGRD
  819.             RE(IB3,1,1)=RE(IB3,1,1)*ZNOR
  820.           ENDDO
  821.           SEGSUP,MLREE2,MLMOT2
  822.         ENDIF
  823.       ENDIF
  824. *  un resultat
  825.       call ecrobj('RIGIDITE',mrigid)
  826.       END
  827.  
  828.  
  829.  

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