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

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