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

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