C PROJRG SOURCE BP208322 10/12/01 21:15:14 6798 SUBROUTINE PROJRG(MRIGID,MTAB1,POS,RI1,RI2) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * *********************************************************** * PROJECTION D'UNE MATRICE SUR UNE BASE DE MODES * * _______________________________________________________ * * * * CREATION : Nicolas BENECH (11 Avril 1995) * * MODIFICATIONS : * * -18/11/2010 Benoit PRABEL : AMELIORATION PERFORMANCE * * _______________________________________________________ * * * * MODULE(S) APPELANT(S) : PJBA * * * * MODULE(S) APPELE(S) : ACCTAB, YTMX * * mucpri, corrsp, xty * * _______________________________________________________ * * * * EN ENTREE : * * MRIGID : Matrice a projeter * * MTAB1 : Base de modes, reels ou complexes * * 'REEL' : indique que l'on utilise le produit * * scalaire reel (pas de conjugaison) * * * * EN SORTIE : * * RI1 : Matrice projetee (partie reelle) * * RI2 : Matrice projetee (partie imaginaire) * * _______________________________________________________ * * * * REMARQUE : * * L'operation realisee est : * * (MTAB1)t . MRIGID. MTAB1 * * Dans le cas complexe, la transposition est accompagnee * * de la conjugaison (si REEL n'est pas mentionne). * *********************************************************** * -INC CCOPTIO -INC SMCHPOI -INC SMELEME -INC SMLCHPO -INC SMRIGID -INC SMTABLE -INC SMLMOTS pointeur IPLMOX.mlmots, IPLMOY.mlmots * * Declarations * REAL*8 XVAL, RMAX CHARACTER*8 LETYPE CHARACTER*8 TYPMOD,CHARRE LOGICAL MODCOM INTEGER I, J, NBMOD, POS, IREEL, IVALRE, IOBRE, isym REAL*8 XVALRE LOGICAL LOGRE * * par defaut, on considere le vrai produit scalaire avec le conjugué * transposé y^H.M.x , mais avec POS=1, on prend seulement y^T.M.x IREEL = -1 IF(POS.EQ.1) IREEL = 1 * * recup de MTYMAT + symetrie / antisymetrie de la matrice a projeter segact,MRIGID LETYPE = MRIGID.MTYMAT NRIGEL=IRIGEL(/2) isym=IRIGEL(7,1) if(NRIGEL.le.1) goto 09 do iel=2,NRIGEL if (isym.ne.IRIGEL(7,iel)) then isym=2 goto 09 endif enddo 09 CONTINUE segdes,MRIGID * LETYPE = ' ' CALL ACCTAB(MTAB1,'MOT',0,0.0D0,'MODES',.TRUE.,0, & 'TABLE',IVALRE,XVALRE,CHARRE,LOGRE,MTAB2) SEGACT, MTAB2 NBMOD = MTAB2.MLOTAB-2 SEGDES, MTAB2 * N1 = NBMOD SEGINI, MLCHP1 SEGINI, MLCHP2 * * Constitution du maillage support et du segment descriptif * NBNN = NBMOD NBELEM = 1 NBSOUS = 0 NBREF = 0 SEGINI, MELEME ITYPEL = 1 * NLIGRD=NBMOD NLIGRP=NBMOD nelrig=1 SEGINI,DESCR * DO 10, I=1, NBMOD IPT1 = 0 * CALL ACCTAB(MTAB2,'ENTIER',I,0.0D0,' ',.TRUE.,0, & 'TABLE',IVALRE,XVALRE,CHARRE,LOGRE,MTAB3) CALL ACCTAB(MTAB3,'MOT',0,0.0D0,'SOUSTYPE',.TRUE.,0, & 'MOT',IVALRE,XVALRE,TYPMOD,LOGRE,IOBRE) * * Le calcul est impossible : * IF (TYPMOD.EQ.'MODE_ANN') THEN IF (LETYPE.NE.'ANNULE') THEN WRITE (*,*) 'Calcul impossible : modes annules !!!' LETYPE = 'ANNULE' ENDIF GOTO 5 ENDIF * * Cas reel ou cas complexe ? * IF (TYPMOD .EQ. 'MODE_COM') THEN MODCOM=.TRUE. CALL ACCTAB(MTAB3,'MOT',0,0.0D0,'DEFORMEE_MODALE_REELLE', & .TRUE.,0,'CHPOINT',IVALRE,XVALRE,CHARRE,LOGRE,MCHPOI) MLCHP1.ICHPOI(I) = MCHPOI CALL ACCTAB(MTAB3,'MOT',0,0.0D0,'DEFORMEE_MODALE_IMAGINAIRE', & .TRUE.,0,'CHPOINT',IVALRE,XVALRE,CHARRE,LOGRE,MCHPOI) MLCHP2.ICHPOI(I) = MCHPOI ELSE MODCOM = .FALSE. CALL ACCTAB(MTAB3,'MOT',0,0.0D0,'DEFORMEE_MODALE',.TRUE.,0, & 'CHPOINT',IVALRE,XVALRE,CHARRE,LOGRE,MCHPOI) MLCHP1.ICHPOI(I) = MCHPOI ENDIF * CALL ACCTAB(MTAB3,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0, & 'POINT',IVALRE,XVALRE,CHARRE,LOGRE,IPT1) * 5 MELEME.NUM(I,1)=IPT1 * DESCR.LISINC(I) = 'ALFA' DESCR.LISDUA(I) = 'FALF' DESCR.NOELEP(I) = I DESCR.NOELED(I) = I * 10 CONTINUE SEGDES, DESCR SEGDES, MELEME * * Constitution des segments XMATRI * NLIGRD=NBMOD NLIGRP=NBMOD * IF (LETYPE .EQ. 'ANNULE') THEN SEGINI, XMATR1 IF (MODCOM) THEN SEGINI, XMATR2 SEGDES, XMATR2 ENDIF SEGDES, XMATR1 GOTO 55 ENDIF * * Cas reel * SEGINI, XMATR1 c DO 20, I=1, NBMOD c MCHPO1 = MLCHP1.ICHPOI(I) c DO 20, J=1, NBMOD c MCHPO2 = MLCHP1.ICHPOI(J) c CALL YTMX (MCHPO2, MCHPO1, MRIGID, XVAL) c XMATR1.RE(I,J,1)=XVAL c 20 CONTINUE c bp 19/11/2010 : on ameliore les perf en remplacant ci dessus par : if (isym.eq.0) then DO 20 J=1, NBMOD MCHPO2 = MLCHP1.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) DO I=J, NBMOD MCHPO1 = MLCHP1.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR1.RE(I,J,1)=XVAL XMATR1.RE(J,I,1)=XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 20 CONTINUE elseif (isym.eq.1) then DO 21 J=1, NBMOD XMATR1.RE(J,J,1)=0.D0 if(J.ge.NBMOD) goto 21 MCHPO2 = MLCHP1.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) JP1=J+1 DO I=JP1, NBMOD MCHPO1 = MLCHP1.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR1.RE(I,J,1)=XVAL XMATR1.RE(J,I,1)=-1.D0*XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 21 CONTINUE else DO 22, J=1, NBMOD MCHPO2 = MLCHP1.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) DO I=1, NBMOD MCHPO1 = MLCHP1.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR1.RE(I,J,1)=XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 22 CONTINUE endif * * Cas complexe : calcul de termes complementaires * IF (MODCOM) THEN c partie reelle = phiR_i^T.M.phiR_j +/- phiI_i^T.M. phiI_j c DO 30, I=1, NBMOD c MCHPO1 = MLCHP2.ICHPOI(I) c DO 30, J=1, NBMOD c MCHPO2 = MLCHP2.ICHPOI(J) c CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL) c XMATR1.RE(I,J,1)=XMATR1.RE(I,J,1)-IREEL*XVAL c 30 CONTINUE DO 30, J=1, NBMOD MCHPO2 = MLCHP2.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) DO I=1, NBMOD MCHPO1 = MLCHP2.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR1.RE(I,J,1)=XMATR1.RE(I,J,1)-IREEL*XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 30 CONTINUE * 2eme matrice (=partie imaginaire) c partie imaginaire = phiR_i^T.M.phiI_j -/+ phiI_i^T.M. phiR_j SEGINI,XMATR2 c DO 40, I=1, NBMOD c MCHPO1 = MLCHP1.ICHPOI(I) c DO 40, J=1, NBMOD c MCHPO2 = MLCHP2.ICHPOI(J) c CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL) c XMATR2.RE(I,J,1)=XVAL c 40 CONTINUE c DO 50, I=1, NBMOD c MCHPO1 = MLCHP2.ICHPOI(I) c DO 50, J=1, NBMOD c MCHPO2 = MLCHP1.ICHPOI(J) c CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL) c XMATR2.RE(I,J,1)=XMATR2.RE(I,J,1)+IREEL*XVAL c 50 CONTINUE DO 40, J=1, NBMOD MCHPO2 = MLCHP2.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) DO I=1, NBMOD MCHPO1 = MLCHP1.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR2.RE(I,J,1)=XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 40 CONTINUE DO 50, J=1, NBMOD MCHPO2 = MLCHP1.ICHPOI(J) CALL MUCPRI(MCHPO2,MRIGID,MCHPO3) DO I=1, NBMOD MCHPO1 = MLCHP2.ICHPOI(I) CALL CORRSP(MRIGID,MCHPO1,MCHPO3,IPLMOX,IPLMOY) CALL XTY1(MCHPO1,MCHPO3,IPLMOX,IPLMOY,XVAL) XMATR2.RE(I,J,1)=XMATR2.RE(I,J,1)+IREEL*XVAL ENDDO segsup,IPLMOX,IPLMOY segsup,MCHPO3 50 CONTINUE SEGDES,XMATR2 ENDIF * SEGDES,XMATR1 SEGDES,MLCHP1,MLCHP2 * * Creation des segments IMATRI * 55 NELRIG = 1 * SEGINI, IMATR1 * IMATR1.IMATTT(1) = XMATR1 * SEGDES, IMATR1 IF (MODCOM) THEN * SEGINI, IMATR2 * IMATR2.IMATTT(1) = XMATR2 * SEGDES, IMATR2 ENDIF * * Creation des rigidites calculees * NRIGE=7 NRIGEL=1 SEGINI, RI1 RI1.MTYMAT = LETYPE RI1.IFORIG = IFOMOD RI1.IMGEO1 = 0 RI1.IMGEO2 = 0 RI1.COERIG = 1.D0 RI1.IRIGEL(1,1) = MELEME RI1.IRIGEL(2,1) = 0 RI1.IRIGEL(3,1) = DESCR RI1.IRIGEL(4,1) = xMATR1 RI1.IRIGEL(5,1) = NIFOUR RI1.IRIGEL(6,1) = 0 c RI1.IRIGEL(7,1) = 2 RI1.IRIGEL(7,1) = isym SEGDES, RI1 IF (MODCOM) THEN RI1.IRIGEL(7,1) = 2 SEGINI, RI2 = RI1 RI2.IRIGEL(4,1) = xMATR2 SEGDES, RI2 ELSE RI2 = 0 SEGSUP, MLCHP2 ENDIF * RETURN END