Télécharger la plaquette Cast3M


fr en

Télécharger devtra.eso

Retour à la liste

Numérotation des lignes : 

devtra

  1. C DEVTRA    SOURCE    CB215821  25/04/23    21:15:14     12247          
  2.       SUBROUTINE DEVTRA(ITBAS,ITKM,ITA,KTKAM,IPMAIL,KTRES,KTNUM,KPREF,
  3.      &                  KTPHI,KTLIAB,RIGIDE,ITCARA,LMODYN)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. *--------------------------------------------------------------------*
  7. *                                                                    *
  8. *     Operateur DYNE : algorithme de Fu - de Vogelaere               *
  9. *     ________________________________________________               *
  10. *                                                                    *
  11. *     Transpose l'information des objets de Castem2000 dans des      *
  12. *     tableaux de travail.                                           *
  13. *                                                                    *
  14. *     Parametres:                                                    *
  15. *                                                                    *
  16. * e   ITBAS   Table representant une base modale                     *
  17. * e   ITKM    Table contenant les matrices XK et XM                  *
  18. * e   ITA     Table contenant la matrice XASM                        *
  19. * es  KTKAM   Segment contenant les matrices XK, XASM et XM          *
  20. * e   IPMAIL  Maillage de reference pour les CHPOINTs resultats      *
  21. * es  KTRES   Segment de sauvegarde des resultats                    *
  22. * e   KPREF   Segment des points de reference                        *
  23. * es  KTPHI   Segment des deformees modales                          *
  24. * e   KTLIAB  Segment des liaisons sur base B                        *
  25. * e   RIGIDE  Vrai si corps rigide, faux sinon                       *
  26. *                                                                    *
  27. *     Auteur, date de creation:                                      *
  28. *                                                                    *
  29. *     Denis ROBERT-MOUGIN, le 26 mai 1989.                           *
  30. *                                                                    *
  31. *--------------------------------------------------------------------*
  32. -INC SMRIGID
  33. -INC SMCOORD
  34. -INC SMELEME
  35.  
  36. -INC PPARAM
  37. -INC CCOPTIO
  38. -INC CCREEL
  39. *
  40.       SEGMENT,MTKAM
  41.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  42.          REAL*8 XOPER(NB1,NB1,NOPER)
  43.       ENDSEGMENT
  44.       SEGMENT,MTPHI
  45.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  46.          INTEGER IAROTA(NSB)
  47.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  48.       ENDSEGMENT
  49.       SEGMENT MTLIAB
  50.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  51.          REAL*8 XPALB(NLIAB,NXPALB)
  52.          REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  53.       ENDSEGMENT
  54.       SEGMENT,MTNUM
  55.          REAL*8 XDT(NPC1),XTEMPS(NPC1)
  56.       ENDSEGMENT
  57.       SEGMENT,MPREF
  58.          INTEGER IPOREF(NPREF)
  59.       ENDSEGMENT
  60. *
  61. *
  62.       segment mtbas
  63.         integer itbmod,lsstru(np1),nsstru
  64.       endsegment
  65.  
  66. c     segment local : enregistre les positions d'indices
  67.       segment IN2IA(LVAL)
  68. c     segment local : calcule l operateur et son inverse
  69.       SEGMENT   MOP
  70.         REAL*8  XOP(NB1,NB1)
  71.         INTEGER INVOP(NB1)
  72.         REAL*8  XOPM1(NB1,NB1)
  73.       ENDSEGMENT
  74.       POINTEUR  MOP2.MOP
  75.  
  76.       LOGICAL L0,L1,RIGIDE,LMODYN,LPLEIN(3)
  77.       CHARACTER*4 CMOT,MOINC
  78.       CHARACTER*8 TYPRET,CHARRE
  79.       CHARACTER*40 MONMOT
  80. *
  81.       MTKAM  = KTKAM
  82.       MTPHI  = KTPHI
  83.       MTLIAB = KTLIAB
  84.       MPREF  = KPREF
  85.       MTNUM  = KTNUM
  86. *
  87.       NPLB  = IBASB(/1)
  88.       NSB   = INMSB(/1)
  89.       NA2   = XPHILB(/3)
  90.       IDIMB = XPHILB(/4)
  91.       NLIAB = IPALB(/1)
  92.       NA1  = XASM(/1)
  93.       NB1K = XK(/2)
  94.       NB1C = XASM(/2)
  95.       NB1M = XM(/2)
  96.       NB1  = XOPER(/1)
  97.       NOPER= XOPER(/3)
  98.       NPREF=IPOREF(/1)
  99. *
  100.       IA1 = 0
  101.       DEUXPI = 2.D0 * XPI
  102.       RIGIDE =.FALSE.
  103.       LPLEIN(1)=.FALSE.
  104.       LPLEIN(2)=.FALSE.
  105.       LPLEIN(3)=.FALSE.
  106. *
  107. *     Traitement des matrices de variables generalisees:
  108. *
  109.       IF (ITBAS.NE.0 .AND.ITKM.EQ.0.AND.(.NOT.LMODYN)) THEN
  110.          IF (IIMPI.EQ.333)
  111.      &   WRITE(IOIMP,*) 'DEVTRA: cas table BASE_DE_MODES, quel type?'
  112.          CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'SOUSTYPE',L0,IP0,
  113.      &                     'MOT',I1,X1,MONMOT,L1,IP1)
  114.          IF (IERR.NE.0) RETURN
  115. *
  116. *        Cas ou la base est unique
  117. *
  118.          IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  119.             IF (IIMPI.EQ.333)
  120.      &      WRITE(IOIMP,*) 'DEVTRA: lecture table BASE_MODALE'
  121. *
  122. *           On recupere la base de modes
  123. *
  124.             CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  125.      &                        'TABLE',I1,X1,' ',L1,IBAS)
  126.             IF (IERR.NE.0) RETURN
  127.             CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,
  128.      &                  RIGIDE,ITCARA,LMODYN,ITKM)
  129. c            IF (RIGIDE) THEN
  130. c                RIGIDE =.FALSE.
  131. c                DO 80 ILIA =1,NLIAB
  132. c                    ITYP = IPALB(ILIA,1)
  133. c                    IF (ITYP.EQ.35.OR.ITYP.EQ.36) THEN
  134. c                        RIGIDE =.TRUE.
  135. c                    ENDIF
  136. c 80             CONTINUE
  137. c            ENDIF
  138.             IF (IERR.NE.0) RETURN
  139. *
  140. *        Cas ou on a un ensemble de bases
  141. *
  142.          ELSE IF (MONMOT(1:17).EQ.'ENSEMBLE_DE_BASES') THEN
  143.             IF (IIMPI.EQ.333)
  144.      &      WRITE(IOIMP,*) 'DEVTRA: lecture table ENSEMBLE_DE_BASES'
  145. *
  146. *           On boucle sur le nombre de bases
  147. *
  148.             IT = 0
  149.             NPLSB = 0
  150.  10         CONTINUE
  151.             TYPRET = ' '
  152.             IT = IT + 1
  153.             CALL ACCTAB(ITBAS,'ENTIER',IT,X0,' ',L0,IP0,
  154.      &                          TYPRET,I1,X1,CHARRE,L1,ITTBAS)
  155.             IF (IERR.NE.0) RETURN
  156.             IF (ITTBAS.NE.0) THEN
  157.                IF (TYPRET.EQ.'TABLE   ') THEN
  158.                   CALL ACCTAB(ITTBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  159.      &                             'TABLE',I1,X1,' ',L1,IBAS)
  160.                   IF (IERR.NE.0) RETURN
  161.                   CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  162.      &                        RIGIDE,ITCARA,LMODYN,ITKM)
  163.                   IF (IERR.NE.0) RETURN
  164.                   NPLSB = MAX(NPLSB,ICOMP)
  165.                   GOTO 10
  166.                ELSE
  167.                   CALL ERREUR(491)
  168.                   RETURN
  169.                ENDIF
  170.             ENDIF
  171.  
  172. *         le segadj n'est plus necessaire
  173. *         on le fait dans dyne26
  174. *                  MP
  175. *           SEGADJ,MTPHI
  176.          ENDIF
  177. *
  178. *     ELSE IF (ITBAS.NE.0.AND.ITKM.NE.0) THEN
  179. *        WRITE(IOIMP,*)
  180. *   &    'DYNE : TBAS et TKM coexistent ---> @ implementer.'
  181. *        IERR = 2
  182. *        RETURN
  183. *
  184.       ELSE IF (LMODYN) THEN
  185.          IF (IIMPI.EQ.333)
  186.      &   WRITE(IOIMP,*) 'DEVTRA: cas table PASAPAS -> appel DYNE26'
  187.          mtbas = itbas
  188.          NPLSB = 0
  189.          do isstru = 1,nsstru
  190.             CALL DYNE26(itbas,KTKAM,KTLIAB,KTPHI,IA1,isstru,ICOMP,
  191.      &                        RIGIDE,ITCARA,LMODYN,ITKM)
  192.                   IF (IERR.NE.0) RETURN
  193.                   NPLSB = MAX(NPLSB,ICOMP)
  194.          enddo
  195. *
  196. *
  197.       ELSE IF (ITKM.NE.0) THEN
  198.          IF (IIMPI.EQ.333)
  199.      &   WRITE(IOIMP,*) 'DEVTRA: cas table RAIDEUR_ET_MASSE'
  200.          TYPRET = ' '
  201.          CALL ACCTAB(ITKM,'MOT',I0,X0,'RAIDEUR',L0,IP0,
  202.      &                   TYPRET,I1,X1,CHARRE,L1,IRIGI)
  203.          IF (IERR.NE.0) RETURN
  204.          IF (IRIGI.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  205.             IF (IIMPI.EQ.333)
  206.      &      WRITE(IOIMP,*) 'DEVTRA: lecture table RAIDEUR ok'
  207. c           nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  208.               TYPRET= ' '
  209.               CALL ACCTAB(ITKM,'MOT',I0,X0,'NATURE_RAIDEUR',L0,IP0,
  210.      &                        TYPRET,I1,X1,CHARRE,L1,IP1)
  211.               IF(TYPRET(1:3).eq.'MOT') THEN
  212.                if(iimpi.EQ.333) write(ioimp,*) 'Nature de K =',CHARRE
  213.                IF(CHARRE(1:6).eq.'PLEINE') THEN
  214.                 LPLEIN(1)=.TRUE.
  215.                 NB1K=NA1
  216.                 NB1=max(NB1,NB1K)
  217.                 SEGADJ,MTKAM
  218.                ELSEIF(CHARRE(1:8).ne.'DIAGONAL') THEN
  219.                 write(ioimp,*) 'Nature de K =',CHARRE,' non comprise !'
  220.                 call erreur(251)
  221.                ENDIF
  222.               ELSEIF(TYPRET(1:8).NE.'        ') THEN
  223.           write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  224.           write(ioimp,*) 'et pas un ',TYPRET
  225.                MOTERR(1:8)='MOT     '
  226.                call erreur(37)
  227.               ELSE
  228.                if(iimpi.EQ.333) write(ioimp,*)'DIAGO par defaut',NB1K
  229.               ENDIF
  230.             MRIGID = IRIGI
  231.             SEGACT,MRIGID
  232.             NRIGI = IRIGEL(/2)
  233.             DO 20 I=1,NRIGI
  234.                COEF = COERIG(I)
  235.                MELEME = IRIGEL(1,I)
  236.                DESCR  = IRIGEL(3,I)
  237.                XMATRI = IRIGEL(4,I)
  238.                SEGACT,DESCR,MELEME,XMATRI
  239.                NRIG = RE(/3)
  240.                LVAL = RE(/1)
  241.                DO 30 IRIG=1,NRIG
  242.                   IF(LPLEIN(1)) segini,IN2IA
  243. c                 boucle sur les lignes (ddls duals)
  244.                   DO 35 IN=1,LVAL
  245.                     INODE=NOELED(IN)
  246.                     IF(INODE.ne.NOELEP(IN)) THEN
  247.                        WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  248.      &                 'primales et duales de la matrice RAIDEUR'
  249.                        CALL ERREUR(47)
  250.                        RETURN
  251.                     ENDIF
  252.                     NNODE=NUM(INODE,IRIG)
  253. c                   position de cette inconnue dans IPOREF de MPREF
  254.                     DO 36 IA=1,NPREF
  255.                       IF (IPOREF(IA).EQ.NNODE) GOTO 39
  256.  36                 CONTINUE
  257.                     write(ioimp,*) 'DEVTRA: Incoherence entre les ',
  258.      &              'points de reference et la matrice RAIDEUR'
  259.                     CALL ERREUR(504)
  260.  39                 CONTINUE
  261. c          write(ioimp,*) 'DEVTRA:   + noeud dual trouvé en position',IA
  262.                     IF(LPLEIN(1)) THEN
  263. c                     on enregistre la position
  264.                       IN2IA(IN)=IA
  265. c                     boucle sur les ddl duals JN >= IN (depuis le coin)
  266.                       DO 37 JN=1,IN
  267.                         IB = IN2IA(JN)
  268. *                       Matrice pleine ...
  269.                         XK(IA,IB) = XK(IA,IB)
  270.      &                              + (RE(IN,JN,IRIG) * COEF)
  271. c                       attention a ne pas remplir 2 fois la diagonale...
  272.                         IF(IA.eq.IB) GOTO 37
  273.                         XK(IB,IA) = XK(IB,IA)
  274.      &                              + (RE(JN,IN,IRIG) * COEF)
  275.  37                   CONTINUE
  276.                     ELSE
  277. *                     Partie diagonale seulement ...
  278.                       XK(IA,1) = XK(IA,1) + (RE(IN,IN,IRIG) * COEF)
  279.                     ENDIF
  280.  35               CONTINUE
  281.                   IF(LPLEIN(1)) segsup,IN2IA
  282.  30               CONTINUE
  283.                SEGDES,XMATRI,MELEME,DESCR
  284.  20            CONTINUE
  285.             SEGDES,MRIGID
  286.          ELSE
  287.             CALL ERREUR(483)
  288.             RETURN
  289.          ENDIF
  290. *
  291.          IF (IIMPI.EQ.333)
  292.      &   WRITE(IOIMP,*) 'DEVTRA: lecture table MASSE'
  293.          TYPRET = ' '
  294.          CALL ACCTAB(ITKM,'MOT',I0,X0,'MASSE',L0,IP0,
  295.      &                   TYPRET,I1,X1,CHARRE,L1,IMASS)
  296.          IF (IERR.NE.0) RETURN
  297.          IF (IMASS.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  298.             IF (IIMPI.EQ.333)
  299.      &      WRITE(IOIMP,*) 'DEVTRA: lecture table MASSE ok'
  300. c           nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  301.               TYPRET= ' '
  302.               CALL ACCTAB(ITKM,'MOT',I0,X0,'NATURE_MASSE',L0,IP0,
  303.      &                        TYPRET,I1,X1,CHARRE,L1,IP1)
  304.               IF(TYPRET(1:3).eq.'MOT') THEN
  305.                if(iimpi.EQ.333) write(ioimp,*) 'Nature de M =',CHARRE
  306.                IF(CHARRE(1:6).eq.'PLEINE') THEN
  307.                 LPLEIN(3)=.TRUE.
  308.                 NB1M=NA1
  309.                 NB1=max(NB1,NB1M)
  310.                 NOPER=3
  311.                 SEGADJ,MTKAM
  312.                ELSEIF(CHARRE(1:8).ne.'DIAGONAL') THEN
  313.                 write(ioimp,*) 'Nature de M =',CHARRE,' non comprise !'
  314.                 call erreur(251)
  315.                ENDIF
  316.               ELSEIF(TYPRET(1:8).NE.'        ') THEN
  317.           write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  318.           write(ioimp,*) 'et pas un ',TYPRET
  319.                MOTERR(1:8)='MOT     '
  320.                call erreur(37)
  321.               ELSE
  322.                 if(iimpi.EQ.333) write(ioimp,*)'DIAGO par defaut',NB1M
  323.               ENDIF
  324.             MRIGID = IMASS
  325.             SEGACT,MRIGID
  326.             NMASS = IRIGEL(/2)
  327.             DO 40 I=1,NMASS
  328.                COEF = COERIG(I)
  329.                MELEME = IRIGEL(1,I)
  330.                DESCR  = IRIGEL(3,I)
  331.                XMATRI = IRIGEL(4,I)
  332.                SEGACT,DESCR,MELEME,XMATRI
  333.                NRIG = RE(/3)
  334.                LVAL = RE(/1)
  335.                DO 50 IRIG=1,NRIG
  336.                   IF(LPLEIN(3)) segini,IN2IA
  337. c                 boucle sur les lignes (ddls duals)
  338.                   DO 55 IN=1,LVAL
  339.                     INODE=NOELED(IN)
  340.                     IF(INODE.ne.NOELEP(IN)) THEN
  341.                        WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  342.      &                 'primales et duales de la matrice MASSE'
  343.                        CALL ERREUR(47)
  344.                        RETURN
  345.                     ENDIF
  346.                     NNODE=NUM(INODE,IRIG)
  347. c                   position de cette inconnue dans IPOREF de MPREF
  348.                     DO 56 IA=1,NPREF
  349.                       IF (IPOREF(IA).EQ.NNODE) GOTO 59
  350.  56                 CONTINUE
  351.                     write(ioimp,*) 'DEVTRA: Incoherence entre les ',
  352.      &              'points de reference et la matrice MASSE'
  353.                     CALL ERREUR(504)
  354.  59                 CONTINUE
  355.                     IF(LPLEIN(3)) THEN
  356. c                     on enregistre la position
  357.                       IN2IA(IN)=IA
  358. c                     boucle sur les ddl duals JN >= IN (depuis le coin)
  359.                       DO 57 JN=1,IN
  360.                         IB = IN2IA(JN)
  361. *                       Matrice pleine ...
  362.                         XM(IA,IB) = XM(IA,IB)
  363.      &                              + (RE(IN,JN,IRIG) * COEF)
  364. c                       attention a ne pas remplir 2 fois la diagonale...
  365.                         IF(IA.eq.IB) GOTO 57
  366.                         XM(IB,IA) = XM(IB,IA)
  367.      &                              + (RE(JN,IN,IRIG) * COEF)
  368.  57                   CONTINUE
  369.                     ELSE
  370. *                     Partie diagonale seulement ...
  371.                       XM(IA,1) = XM(IA,1) + (RE(IN,IN,IRIG) * COEF)
  372.                     ENDIF
  373.  55               CONTINUE
  374.                   IF(LPLEIN(3)) segsup,IN2IA
  375.  50            CONTINUE
  376.                SEGDES,XMATRI,MELEME,DESCR
  377.  40         CONTINUE
  378.             SEGDES,MRIGID
  379.          ELSE
  380.             CALL ERREUR(484)
  381.             RETURN
  382.          ENDIF
  383.  
  384. *        traitement de la base modale (necessaire si liaison B)
  385. *        --> remplissage de XPHILB
  386.          TYPRET = ' '
  387.          CALL ACCTAB(ITKM,'MOT',I0,X0,'BASE_MODALE',L0,IP0,
  388.      &               TYPRET,I1,X1,CHARRE,L1,ITBAS2)
  389.          IF(ITBAS2.eq.0) GOTO 599
  390.          TYPRET = ' '
  391.          CALL ACCTAB(ITBAS2,'MOT',IMODE,X0,'SOUSTYPE',L0,IP0,
  392.      &                      'MOT',I1,X1,MONMOT,L1,IP1)
  393.          IF (IERR.NE.0) RETURN
  394. *        Cas ou la base est unique
  395.          IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  396. *           On recupere la base de modes
  397.             CALL ACCTAB(ITBAS2,'MOT',IMODE,X0,'MODES',L0,IP0,
  398.      &                        'TABLE',I1,X1,' ',L1,IBAS2)
  399.             IF (IERR.NE.0) RETURN
  400.             CALL DYNE26(IBAS2,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,
  401.      &                  RIGIDE,ITCARA,LMODYN,ITKM)
  402. *        Cas ou la base est un ensemble de bases
  403.          ELSEIF(MONMOT(1:17).EQ.'ENSEMBLE_DE_BASES') THEN
  404.             IT = 0
  405.             NPLSB = 0
  406.  510        CONTINUE
  407.             TYPRET = ' '
  408.             IT = IT + 1
  409.             CALL ACCTAB(ITBAS2,'ENTIER',IT,X0,' ',L0,IP0,
  410.      &                          TYPRET,I1,X1,CHARRE,L1,ITTBAS)
  411.             IF (IERR.NE.0) RETURN
  412.             IF (ITTBAS.NE.0) THEN
  413.                IF (TYPRET.EQ.'TABLE   ') THEN
  414.                   CALL ACCTAB(ITTBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  415.      &                             'TABLE',I1,X1,' ',L1,IBAS2)
  416.                   IF (IERR.NE.0) RETURN
  417.                   CALL DYNE26(IBAS2,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  418.      &                        RIGIDE,ITCARA,LMODYN,ITKM)
  419.                   IF (IERR.NE.0) RETURN
  420.                   NPLSB = MAX(NPLSB,ICOMP)
  421.                   GOTO 510
  422.                ELSE
  423.                   CALL ERREUR(491)
  424.                   RETURN
  425.                ENDIF
  426.             ENDIF
  427.          ENDIF
  428.  599     CONTINUE
  429.          IF (IERR.NE.0) RETURN
  430.  
  431.       ENDIF
  432. *
  433. *     Traitement de la matrice d'amortissement
  434. *
  435.       IF (ITA.NE.0) THEN
  436.          IF (IIMPI.EQ.333)
  437.      &   WRITE(IOIMP,*) 'DEVTRA: cas table AMORTISSEMENT'
  438.          if (lmodyn) then
  439.            iamor = ita
  440.            typret='RIGIDITE'
  441.          else
  442.            TYPRET = ' '
  443.            CALL ACCTAB(ITA,'MOT',I0,X0,'AMORTISSEMENT',L0,IP0,
  444.      &                    TYPRET,I1,X1,CHARRE,L1,IAMOR)
  445.            IF (IERR.NE.0) RETURN
  446.          endif
  447.          IF (IAMOR.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  448.             IF (IIMPI.EQ.333)
  449.      &      WRITE(IOIMP,*) 'DEVTRA: lecture table AMORTISSEMENT ok'
  450. c           nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  451.             if(.not.lmodyn) then
  452.               TYPRET= ' '
  453.               CALL ACCTAB(ITA,'MOT',I0,X0,'NATURE',L0,IP0,
  454.      &                        TYPRET,I1,X1,CHARRE,L1,IP1)
  455.               IF(TYPRET(1:3).eq.'MOT') THEN
  456.                if(iimpi.EQ.333) write(ioimp,*) 'Nature de C =',CHARRE
  457.                IF(CHARRE(1:6).eq.'PLEINE') THEN
  458.                 LPLEIN(2)=.TRUE.
  459.                 NB1C=NA1
  460.                 NB1=max(NB1,NB1C)
  461.                 NOPER=max(2,NOPER)
  462.                 SEGADJ,MTKAM
  463.                ELSEIF(CHARRE(1:8).ne.'DIAGONAL') THEN
  464.                 write(ioimp,*) 'Nature de C =',CHARRE,' non comprise !'
  465.                 call erreur(251)
  466.                ENDIF
  467.               ELSEIF(TYPRET(1:8).NE.'        ') THEN
  468.           write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  469.           write(ioimp,*) 'et pas un ',TYPRET
  470.                MOTERR(1:8)='MOT     '
  471.                call erreur(37)
  472.               ENDIF
  473.             endif
  474.             MRIGID = IAMOR
  475.             SEGACT,MRIGID
  476.             NAMOR = IRIGEL(/2)
  477.             DO 60 I=1,NAMOR
  478.                COEF = COERIG(I)
  479. c          write(ioimp,*) 'DEVTRA: sous rigidite ',I,'/',NAMOR,COEF
  480.                MELEME = IRIGEL(1,I)
  481.                DESCR  = IRIGEL(3,I)
  482.                XMATRI = IRIGEL(4,I)
  483.                SEGACT,DESCR,MELEME,XMATRI
  484.                NRIG = RE(/3)
  485.                LVAL = RE(/1)
  486.                DO 70 IRIG=1,NRIG
  487. c          write(ioimp,*) 'DEVTRA: + element',IRIG,'/',NRIG
  488.                   IF(LPLEIN(2)) segini,IN2IA
  489. c                 boucle sur les lignes (ddls duals)
  490.                   DO 75 IN=1,LVAL
  491.                     INODE=NOELED(IN)
  492.                     IF(INODE.ne.NOELEP(IN)) THEN
  493.                        WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  494.      &                 'primales et duales de la matrice AMORTISSEMENT'
  495.                        CALL ERREUR(47)
  496.                        RETURN
  497.                     ENDIF
  498.                     NNODE=NUM(INODE,IRIG)
  499. c          write(ioimp,*) 'DEVTRA:   + noeud dual',IN,'/',LVAL,' #',NNODE
  500. c                   position de cette inconnue dans IPOREF de MPREF
  501.                     DO 76 IA=1,NPREF
  502.                       IF (IPOREF(IA).EQ.NNODE) GOTO 79
  503.  76                 CONTINUE
  504.                     write(ioimp,*) 'DEVTRA: Incoherence entre les ',
  505.      &              'points de reference et la matrice AMORTISSEMENT'
  506.                     CALL ERREUR(504)
  507.  79                 CONTINUE
  508. c          write(ioimp,*) 'DEVTRA:   + noeud dual trouvé en position',IA
  509.                     IF(LPLEIN(2)) THEN
  510. c                     on enregistre la position
  511.                       IN2IA(IN)=IA
  512. c                     boucle sur les ddl duals JN >= IN (depuis le coin)
  513.                       DO 77 JN=1,IN
  514.                         IB = IN2IA(JN)
  515. *                       Matrice pleine ...
  516.                         XASM(IA,IB) = XASM(IA,IB)
  517.      &                              + (RE(IN,JN,IRIG) * COEF)
  518. c                       attention a ne pas remplir 2 fois la diagonale...
  519.                         IF(IA.eq.IB) GOTO 77
  520.                         XASM(IB,IA) = XASM(IB,IA)
  521.      &                              + (RE(JN,IN,IRIG) * COEF)
  522.  77                   CONTINUE
  523.                     ELSE
  524. *                     Partie diagonale seulement ...
  525.                       XASM(IA,1) = XASM(IA,1) + (RE(IN,IN,IRIG) * COEF)
  526.                     ENDIF
  527.  75               CONTINUE
  528.                   IF(LPLEIN(2)) segsup,IN2IA
  529.  70            CONTINUE
  530.                SEGDES,XMATRI,MELEME,DESCR
  531.  60         CONTINUE
  532.             SEGDES,MRIGID
  533.          ELSE
  534.             CALL ERREUR(485)
  535.             RETURN
  536.          ENDIF
  537.       ENDIF
  538.  
  539. *     on calcule les operateurs = inverses de :
  540. *     1:      XOP = 4I + dt*M^-1*C
  541. *     2: MOP2.XOP = 6I + dt*M^-1*C
  542. *     3:            M
  543.  
  544.       IF(NOPER.GT.0) THEN
  545.          IF (IIMPI.EQ.333)
  546.      &   WRITE(IOIMP,*)'DEVTRA : calcul des operateurs'
  547.          SEGINI,MOP,MOP2
  548. c        le pas de temps est constant (seule possibilite aujourd'hui)
  549.          pdt = xdt(1)
  550.  
  551. *     ---si M pleine, on l'inverse---
  552.          IF(LPLEIN(3)) THEN
  553. c          Inversion de M + stockage dans XOPER( , ,3)
  554.            CALL IVMAT(NA1,XM,INVOP,XOPM1,DETOP,0,IRET)
  555.            IF(IRET.ne.0) RETURN
  556.            DO 610 IA=1,NA1
  557.            DO 610 IB=1,NA1
  558.               XOPER(IA,IB,3)= XOPM1(IA,IB)
  559.  610       CONTINUE
  560.          ENDIF
  561.  
  562. *     ---Calcul des operateurs---
  563.          DO 611 IA=1,NA1
  564.                   XOP(IA,IA) = 4.D0
  565.              MOP2.XOP(IA,IA) = 6.D0
  566.  611     CONTINUE
  567. c     ...C pleine, M pleine
  568.          IF(LPLEIN(2).AND.LPLEIN(3)) THEN
  569.            DO 612 IA=1,NA1
  570.            DO 612 IB=1,NA1
  571.               AUX = 0.D0
  572.               DO 613 IC=1,NA1
  573.                  AUX = AUX + XOPM1(IA,IC)*XASM(IC,IB)
  574.  613          CONTINUE
  575.               AUX = pdt*AUX
  576.                    XOP(IA,IB) =      XOP(IA,IB) + AUX
  577.               MOP2.XOP(IA,IB) = MOP2.XOP(IA,IB) + AUX
  578.  612       CONTINUE
  579. c     ...C pleine, M diago
  580.          ELSEIF(LPLEIN(2)) THEN
  581.            DO 614 IA=1,NA1
  582.            DO 614 IB=1,NA1
  583.               AUX = pdt*XASM(IA,IB)/XM(IA,1)
  584.                    XOP(IA,IB) =      XOP(IA,IB) + AUX
  585.               MOP2.XOP(IA,IB) = MOP2.XOP(IA,IB) + AUX
  586.  614       CONTINUE
  587. c     ...C diago, M pleine
  588.          ELSEIF(LPLEIN(3)) THEN
  589.            DO 615 IA=1,NA1
  590.            DO 615 IB=1,NA1
  591.               AUX = pdt*XOPM1(IA,IB)*XASM(IB,1)
  592.                    XOP(IA,IB) =      XOP(IA,IB) + AUX
  593.               MOP2.XOP(IA,IB) = MOP2.XOP(IA,IB) + AUX
  594.  615       CONTINUE
  595.          ELSE
  596. c          pour le cas tout diago, on na pas besoin d'operateur
  597.            WRITE(IOIMP,*) 'DEVTRA: option PLEINE incoherente !'
  598.            CALL ERREUR(5)
  599.            RETURN
  600.          ENDIF
  601.  
  602. *     ---Inversion des operateurs---
  603.          CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET)
  604.          IF(IRET.ne.0) RETURN
  605.          CALL IVMAT(NA1,MOP2.XOP,MOP2.INVOP,MOP2.XOPM1,DETOP2,0,
  606.      &              IRET)
  607.          IF(IRET.ne.0) RETURN
  608. *        rem : ici IVMAT, mais il existe aussi INVER, INVER1 ...
  609.          DO 69 IA=1,NA1
  610.          DO 69 IB=1,NB1
  611.             XOPER(IA,IB,1)=     XOPM1(IA,IB)
  612.             XOPER(IA,IB,2)=MOP2.XOPM1(IA,IB)
  613.  69      CONTINUE
  614.          SEGSUP,MOP,MOP2
  615.       ENDIF
  616. *
  617.       IF (IIMPI.EQ.333) THEN
  618.          WRITE(IOIMP,*)'     segment MTPHI'
  619.          WRITE(IOIMP,*)'DEVTRA : valeur de NPLB  :',IBASB(/1)
  620.          WRITE(IOIMP,*)'DEVTRA : valeur de NSB   :',XPHILB(/1)
  621.          WRITE(IOIMP,*)'DEVTRA : valeur de NPLSB :',XPHILB(/2)
  622.          WRITE(IOIMP,*)'DEVTRA : valeur de NA2   :',XPHILB(/3)
  623.          WRITE(IOIMP,*)'DEVTRA : valeur de IDIMB :',XPHILB(/4)
  624.          WRITE(IOIMP,*)'     segment MTKAM'
  625.          WRITE(IOIMP,*)'NA1,NB1K,NB1C,NB1M,NB1,NOPER='
  626.          WRITE(IOIMP,*) NA1,NB1K,NB1C,NB1M,NB1,NOPER
  627.          WRITE(IOIMP,*) 'LPLEIN=',(LPLEIN(jou),jou=1,3)
  628.          if(NB1K.gt.1) then
  629.             do iou=1,NA1
  630.                WRITE(IOIMP,*) 'XK=',(XK(iou,jou),jou=1,NB1K)
  631.             enddo
  632.          else
  633.             do iou=1,NA1
  634.                WRITE(IOIMP,*) 'XK(',iou,',1)=',XK(iou,1)
  635.             enddo
  636.          endif
  637.          if(NB1C.gt.1) then
  638.             do iou=1,NA1
  639.               WRITE(IOIMP,*) 'XASM=',(XASM(iou,jou),jou=1,NB1C)
  640.             enddo
  641.          else
  642.             do iou=1,NA1
  643.               WRITE(IOIMP,*) 'XASM(',iou,',1)=',XASM(iou,1)
  644.             enddo
  645.          endif
  646.          if(NB1M.gt.1) then
  647.             do iou=1,NA1
  648.                WRITE(IOIMP,*) 'XM=',(XM(iou,jou),jou=1,NB1M)
  649.             enddo
  650.          else
  651.             do iou=1,NA1
  652.                WRITE(IOIMP,*) 'XM(',iou,',1)=',XM(iou,1)
  653.             enddo
  654.          endif
  655.          if(NOPER.ge.1) then
  656.          do iop=1,NOPER
  657.            do iou=1,NB1
  658.              WRITE(IOIMP,*) 'XOPER(',iou,',:,',iop,')=',
  659.      &                      (XOPER(iou,jou,iop),jou=1,NB1)
  660.            enddo
  661.          enddo
  662.          endif
  663.       ENDIF
  664. *
  665.       END
  666.  
  667.  
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  
  683.  
  684.  

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