Télécharger la plaquette Cast3M


fr en

Télécharger d2vtra.eso

Retour à la liste

Numérotation des lignes : 

d2vtra

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

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