C TUMU SOURCE PV090527 26/04/30 21:16:45 12529 SUBROUTINE TUMU(LCH1,IRIG1,DFLO,IRIG2) ************************************************************************ * NOM : TUMU ************************************************************************ * DESCRIPTION : Realise le produit tU*M*U ou M est une matrice carree * et U est une matrice rectangle dont les colonnes sont * donnees par un objet LISTCHPO * * Les multiplicateurs de Lagrange sont ignores * * Les inconnues duales de M doivent avoir ete definies * dans la liste NOMDU de l'include CCHAMP pour savoir * comment effectuer la multiplication a gauche * *********************************************************** * * * U[N;L] * | * | * V * +-----------+ +---+---+---+---+ * | | | C | C | | C | * | M | | H | H | . | H | * | | | P | P | . | P | * tU[L;N] | [N;N] | | # | # | . | # | * | | | | 1 | 2 | | L | * | +-----------+ +---+---+---+---+ * V * +-----------+ +-----------+ +---------------+ * | CHPOINT#1 | | CHP#1 * M | | | * +-----------+ +-----------+ | | * | CHPOINT#2 | | CHP#2 * M | | TUMU | * +-----------+ +-----------+ | | * | ... | | ... | | [L;L] | * +-----------+ +-----------+ | | * | CHPOINT#L | | CHP#L * M | | | * +-----------+ +-----------+ +---------------+ * * * avec : L = nombre de champs * N = nombre d'inconnues * (triplet noeud/composante/harmonique) * ************************************************************************ * APPELE PAR : pod.eso ************************************************************************ * ENTREES :: LCH1 = POINTEUR VERS UN OBJET LISTCHPO * IRIG1 = POINTEUR VERS UN OBJET RIGIDITE * DFLO = COEFFICIENT MULTIPLICATEUR * SORTIES :: IRIG2 = POINTEUR VERS UN OBJET RIGIDITE ************************************************************************ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMCHPOI -INC SMLCHPO -INC SMELEME -INC SMRIGID -INC CCHAMP -INC CCREEL * * ICPR(I) EST LE NUMERO LOCAL (DANS LE SUPPORT GEOMETRIQUE DU * LISTCHPO LCH1) DU I-EME NOEUD GLOBAL (DANS LA TABLE XCOORD) SEGMENT/TICPR/(ICPR(NOMAX)) * * TINCO(I) = [MCOMP(I) ; IHARM(I)] * TINCO DEFINIT LES NBINC INCONNUES DE LA MATRICE IRIG1 * ******************************** * MCOMP(I) EST LE NOM DE COMPOSANTE DE LA I-EME INCONNUE * IHARM(I) EST LE NUMERO D'HARMONIQUE DE LA I-EME INCONNUE SEGMENT,TCOMP CHARACTER*(LOCOMP) MCOMP(0) ENDSEGMENT SEGMENT,IHARM(0) * * IMAI(I) EST LE MAILLAGE ASSOCIE AU I-EME SOUPO DU PREMIER CHPOINT * IGLO(I,J) EST L'INCONNUE DE TINCO QUI CORRESPOND A L'INCONNUE * I DU SOUPO J DU PREMIER CHPOINT * IPOV(I,J) EST LE MPOVAL DU SOUPO DU J-EME CHPOINT QUI CORRESPOND * AU I-EME SOUPO DU PREMIER CHPOINT (IPOV=0 SI LE SOUPO * EST ASSOCIE AUX MULT. DE LAGRANGE) * LA COMPOSANTE I DU SOUPO J DU PREMIER CHPOINT EST EN POSITION * IINC(I,J,K) DANS LE SOUPO CORRESPONDANT DU K-EME CHPOINT SEGMENT TRAV1 INTEGER IMAI(NBSOU) INTEGER IGLO(NXMAX,NBSOU) INTEGER IPOV(NBSOU,NBCHP) INTEGER IINC(NXMAX,NBSOU,NBCHP) ENDSEGMENT * * VPO(I,J,K) EST LA VALEUR DU K-EME CHPOINT PRISE POUR LE I-EME * NOEUD LOCAL ET POUR LA J-EME INCONNUE DE TINCO * CC(I) ET DD(I) SONT LES VALEURS EXTRAITES DE VPO POUR FORMER LES * VECTEURS ELEMENTAIRES (ASSOCIES A L'ELEMENT IEL) * QUI SERONT MULTIPLIES PAR LE TABLEAU RE(*,*,IEL) * IPOSP(I) DONNE LA POSITION DANS TINCO DE LA I-EME INCONNUE PRIMALE * D'UNE MATRICE ELEMENTAIRE * IPOSD(I) DONNE LA POSITION DANS TINCO DE LA I-EME INCONNUE DUALE * D'UNE MATRICE ELEMENTAIRE SEGMENT TRAV2 REAL*8 VPO(NBPOI,NBINC,NBCHP) REAL*8 CC(NLIGMA),DD(NLIGMA) INTEGER IPOSP(NLIGMA) INTEGER IPOSD(NLIGMA) ENDSEGMENT * CHARACTER*(LOCOMP) MOCOMP * * * * +---------------------------------------------------------------+ * | | * | T R A V A I L P R E L I M I N A I R E | * | | * +---------------------------------------------------------------+ * MLCHPO=LCH1 SEGACT,MLCHPO * * ========================================================== * CONSTRUCTION DE LA LISTE DES INCONNUES DE LA MATRICE IRIG1 * (= COUPLE NOM_DE_COMPOSANTE_PRIMALE + NUMERO_HARMONIQUE) * ========================================================== * MRIGID=IRIG1 SEGACT,MRIGID NRIGEL=IRIGEL(/2) * SEGINI,TCOMP,IHARM NBINC=0 * * BOUCLE SUR LES RIGIDITES ELEMENTAIRES NLIGMA=0 DO IRI=1,NRIGEL DESCR=IRIGEL(3,IRI) NUHAR=IRIGEL(5,IRI) SEGACT,DESCR NLIGRE=LISINC(/2) IF (NLIGRE.GT.NLIGMA) NLIGMA=NLIGRE DO 1 I1=1,NLIGRE MOCOMP=LISINC(I1) DO 2 I2=1,NBINC IF (MOCOMP.NE.MCOMP(I2)) GOTO 2 IF (NUHAR .EQ.IHARM(I2)) GOTO 1 2 CONTINUE MCOMP(**)=MOCOMP IHARM(**)=NUHAR NBINC=NBINC+1 1 CONTINUE ENDDO * * * =================================== * NOMBRE DE CHPOINTS DANS LE LISTCHPO * =================================== * NBCHP=ICHPOI(/1) IF (NBCHP.EQ.0) THEN MOTERR(1:8)='LISTCHPO' INTERR(1)=LCH1 CALL ERREUR(356) RETURN ENDIF * * * =============================================================== * CORRESPONDANCE ENTRE LES INCONNUES DU PREMIER CHPOINT ET CELLES * DES CHPOINTS SUIVANTS + CORRESPONDANCE AVEC LA LISTE GLOBALE * TINCO DES INCONNUES DE LA MATRICE (REMPLISSAGE DE TRAV1) * CORRESPONDANCE ENTRE LES NUMEROTATIONS LOCALE ET GLOBALE * (REMPLISSAGE DE ICPR) * =============================================================== * MCHPO1=ICHPOI(1) SEGACT,MCHPO1 * NBSOU=MCHPO1.IPCHP(/1) IF (NBSOU.EQ.0) THEN MOTERR(1:8)='CHPOINT' CALL ERREUR(1027) RETURN ENDIF * NXMAX=3 SEGINI,TRAV1 * segact mcoord*mod NOMAX=nbpts SEGINI,TICPR NBPOI=0 * * * ************************************** * BOUCLE 1 SUR LES SOUPOS DU 1ER CHPOINT * ************************************** DO 10 IS1=1,NBSOU MSOUP1=MCHPO1.IPCHP(IS1) SEGACT,MSOUP1 * * ON IGNORE LES MULTIPLICATEURS DE LAGRANGE NX1=MSOUP1.NOCOMP(/2) MOCOMP=MSOUP1.NOCOMP(1) IF (MOCOMP.EQ.'LX'.OR.MOCOMP.EQ.'FLX') THEN SEGDES,MSOUP1 GOTO 10 ENDIF * IF (NX1.GT.NXMAX) THEN NXMAX=NX1 SEGADJ,TRAV1 ENDIF * * CORRESPONDANCE ENTRE LES INCONNUES DU SOUPO ET CELLES DE TINCO DO IX1=1,NX1 MOCOMP=MSOUP1.NOCOMP(IX1) NOHA =MSOUP1.NOHARM(IX1) DO 11 IX2=1,NBINC IF (MOCOMP.NE.MCOMP(IX2)) GOTO 11 IF (NOHA .NE.IHARM(IX2)) GOTO 11 IGLO(IX1,IS1)=IX2 11 CONTINUE IINC(IX1,IS1,1)=IX1 ENDDO * * ON VERIFIE QUE LE MAILLAGE N'EST PAS VIDE IGEO1=MSOUP1.IGEOC IF (IGEO1.LE.0) THEN MOTERR(1:8)='CHPOINT' CALL ERREUR(1027) RETURN ENDIF IMAI(IS1)=IGEO1 IPT1=IGEO1 SEGACT,IPT1 NNO1=IPT1.NUM(/2) IF (NNO1.EQ.0) GOTO 10 * * CONSTRUCTION DE LA TABLE ICPR * (NUMEROTATION GLOBALE <=> LOCALE) DO 12 I1=1,NNO1 K1=IPT1.NUM(1,I1) IF (ICPR(K1).NE.0) GOTO 12 NBPOI=NBPOI+1 ICPR(K1)=NBPOI 12 CONTINUE SEGDES,IPT1 * * POINTEUR DIRECT VERS LE SEGMENT MPOVAL IPOV(IS1,1)=MSOUP1.IPOVAL * * * ******************************** * BOUCLE 2 SUR LES AUTRES CHPOINTS * ******************************** DO 14 ICH=2,NBCHP MCHPO2=ICHPOI(ICH) SEGACT,MCHPO2 NS2=MCHPO2.IPCHP(/1) * * ********************************************** * ON VA CHERCHER LE SOUPO CORRESPONDANT A MSOUP1 * => BOUCLE 3 SUR LES SOUPOS DE MCHPO2 * ********************************************** DO 13 IS2=1,NS2 MSOUP2=MCHPO2.IPCHP(IS2) SEGACT,MSOUP1,MSOUP2 * * MEME MAILLAGE ? IGEO2=MSOUP2.IGEOC IF (IGEO1.NE.IGEO2) THEN SEGDES,MSOUP2 GOTO 13 ENDIF * * MEME NOMBRE DE COMPOSANTES ? NX2=MSOUP2.NOCOMP(/2) MOCOMP=MSOUP1.NOCOMP(1) IF (MOCOMP.EQ.'LX'.OR.MOCOMP.EQ.'FLX'.OR.NX1.NE.NX2) THEN SEGDES,MSOUP2 GOTO 13 ENDIF IF (NX2.GT.NXMAX) THEN NXMAX=NX2 SEGADJ,TRAV1 ENDIF * * MEMES LISTES DE COMPOSANTES ? * => ON FAIT LA CORRESPONDANCE ENTRE LES COMPOSANTES DES * 2 SOUPOS DO 15 IX1=1,NX1 MOCOMP=MSOUP1.NOCOMP(IX1) DO 16 IX2=1,NX2 IF (MOCOMP.NE.MSOUP2.NOCOMP(IX2)) GOTO 16 IF (MSOUP1.NOHARM(IX1).EQ.MSOUP2.NOHARM(IX2)) THEN IINC(IX1,IS1,ICH)=IX2 GOTO 15 ENDIF 16 CONTINUE GOTO 19 15 CONTINUE * * POINTEUR DIRECT VERS LE SEGMENT MPOVAL IPOV(IS1,ICH)=MSOUP2.IPOVAL * * (CHPOINT SUIVANT) SEGDES,MSOUP2,MCHPO2 GOTO 14 * 13 CONTINUE * * MESSAGE D'ERREUR * **************** 19 CONTINUE WRITE(MOTERR(1:16),FMT='(2I8)') MCHPO1,MCHPO2 CALL ERREUR(135) RETURN * 14 CONTINUE SEGDES,MSOUP1 * 10 CONTINUE SEGDES,MCHPO1,MLCHPO * * * ======================================================== * STOCKAGE DES VALEURS DU LISTCHPO DANS UN TABLEAU ORDONNE * SELON LA GEOMETRIE LOCALE (ICPR) ET SELON LES INCONNUES * DE LA MATRICE D'ENTREE (TINCO) => REMPLISSAGE DE TRAV2 * ======================================================== * SEGINI,TRAV2 DO ICH=1,NBCHP DO 20 ISOU=1,NBSOU IPO1=IPOV(ISOU,ICH) IF (IPO1.EQ.0) GOTO 20 MPOVAL=IPO1 MELEME=IMAI(ISOU) SEGACT,MELEME,MPOVAL NNO=VPOCHA(/1) NIX=VPOCHA(/2) DO 21 IX=1,NIX IX1=IGLO(IX,ISOU) IF (IX1.EQ.0) GOTO 21 IIX=IINC(IX,ISOU,ICH) DO INO=1,NNO N1=ICPR(NUM(1,INO)) VPO(N1,IX1,ICH)=VPOCHA(INO,IIX) ENDDO 21 CONTINUE SEGDES,MELEME,MPOVAL 20 CONTINUE ENDDO SEGSUP,TRAV1 * * * * +---------------------------------------------------------------+ * | | * | C R E A T I O N D U S U P E R - E L E M E N T | * | | * +---------------------------------------------------------------+ * NBSOUS=0 NBELEM=1 NBNN=NBCHP NBREF=0 SEGINI,IPT1 IPT1.ITYPEL=28 segact mcoord*mod NBPT1=nbpts NBPTS=NBPT1+NBNN SEGADJ,MCOORD DO K=1,NBNN K1=(NBPT1+K-1)*(IDIM+1) XCOOR(K1+1)=K XCOOR(K1+2)=0 IF (IDIM.EQ.3) XCOOR(K1+3)=0 IPT1.NUM(K,1)=NBPT1+K ENDDO SEGDES,IPT1 * * * * +---------------------------------------------------------------+ * | | * | D E S C R I P T E U R D E L A M A T R I C E | * | | * +---------------------------------------------------------------+ * NLIGRP=NBCHP NLIGRD=NBCHP SEGINI,DES1 DO K=1,NBCHP DES1.LISINC(K)='ALFA' DES1.LISDUA(K)='FALF' DES1.NOELEP(K)=K DES1.NOELED(K)=K ENDDO SEGDES,DES1 * * * * +---------------------------------------------------------------+ * | | * | R E M P L I S S A G E D U C O N T E N U | * | | * +---------------------------------------------------------------+ * NELRIG=1 rigrel=0 SEGINI,XMATR1 * ISYM1=IRIGEL(7,1) * * (BOUCLE 1 SUR LES SOUS-MATRICES) DO 30 IRI=1,NRIGEL MELEME=IRIGEL(1,IRI) SEGACT,MELEME IF (ITYPEL.EQ.22) GOTO 30 * DESCR=IRIGEL(3,IRI) SEGACT,DESCR NINCP=LISINC(/2) NINCD=LISDUA(/2) * ON VERIFIE QUE LA MATRICE EST CARREE IF (NINCP.NE.NINCD) THEN CALL ERREUR(21) RETURN ENDIF * * ON VERIFIE LA CORRESPONDANCE ENTRE LA LIGNE ET LA COLONNE D'UN * * NOEUD DONNE * DO K=1,NINCP * IF (NOELED(K).NE.NOELEP(K)) THEN * CALL ERREUR(21) * RETURN * ENDIF * ENDDO * NUHAR=IRIGEL(5,IRI) * ISYME=IRIGEL(7,IRI) IF (ISYME.NE.ISYM1) ISYM1=2 * * * ========================================================= * CORRESPONDANCE ENTRE LES INCONNUES PRIMALES DE LA MATRICE * ELEMENTAIRE ET LES INCONNUES "GLOBALES" DU SEGMENT TINCO * (SELON LESQUELLES SONT ORDONNES TOUS LES VPOCHA) * ========================================================= * DO 31 IN1=1,NINCP MOCOMP=LISINC(IN1) DO 32 IN2=1,NBINC IF (MOCOMP.NE.MCOMP(IN2)) GOTO 32 IF (NUHAR .EQ.IHARM(IN2)) THEN IPOSP(IN1)=IN2 GOTO 31 ENDIF 32 CONTINUE 31 CONTINUE * * * ======================================================== * CORRESPONDANCE ENTRE LES INCONNUES DUALES DE LA MATRICE * ELEMENTAIRE ET LES INCONNUES "GLOBALES" DU SEGMENT TINCO * L'ASSOCIATION EST FAITE GRACE AUX LISTES DEFINIES DANS * L'INCLUDE CCHAMP * ======================================================== * DO 33 IN1=1,NINCD MOCOMP=LISDUA(IN1) * ON CHERCHE LA PRIMALE ASSOCIEE A LA DUALE NUMERO IN1 DO IC1=1,LNOMDU IF (MOCOMP.EQ.NOMDU(IC1)) GOTO 34 ENDDO * ERREUR : COMPOSANTE NON REFERENCEE DANS CCHAMP MOTERR=MOCOMP CALL ERREUR(108) RETURN 34 CONTINUE MOCOMP=NOMDD(IC1) DO 35 IN2=1,NBINC IF (MOCOMP.NE.MCOMP(IN2)) GOTO 35 IF (NUHAR .EQ.IHARM(IN2)) THEN IPOSD(IN1)=IN2 GOTO 33 ENDIF 35 CONTINUE CALL ERREUR(21) RETURN 33 CONTINUE * * * ================================================ * CALCUL DES PRODUITS MATRICE/VECTEUR ELEMENTAIRES * ================================================ * NNO=NUM(/1) NEL=NUM(/2) COER=COERIG(IRI) XMATRI=IRIGEL(4,IRI) SEGACT,XMATRI * * ******************************************** * BOUCLE 2 SUR LES ELEMENTS DE LA SOUS-MATRICE * ******************************************** DO 36 IEL=1,NEL * * ON VERIFIE QUE LA RIGIDITE ELEMENTAIRE POSSEDE UN SUPPORT * GEOMETRIQUE COMPATIBLE AVEC LE LISTCHPO DO INO=1,NNO IF (ICPR(NUM(INO,IEL)).NE.0) GOTO 37 ENDDO GOTO 36 37 CONTINUE * * ************************************ * BOUCLE 3 SUR LES CHPOINTS "A DROITE" * ************************************ DO ICH=1,NBCHP * * FABRICATION DU VECTEUR ELEMENTAIRE POUR LA * MULTIPLICATION A DROITE DO INP=1,NINCP I1=ICPR(NUM(NOELEP(INP),IEL)) IF (I1.EQ.0) THEN CC(INP)=0.D0 ELSE I2=IPOSP(INP) CC(INP)=VPO(I1,I2,ICH) ENDIF ENDDO * * REMPLISSAGE DE LA DIAGONALE VA=XTMXMU(CC,RE(1,1,IEL),NINCP,ISYME)*COER XMATR1.RE(ICH,ICH,1)=XMATR1.RE(ICH,ICH,1)+VA * * * REMPLISSAGE DU TRIANGLE SUPERIEUR * ************************************ * BOUCLE 4 SUR LES CHPOINTS "A GAUCHE" * ************************************ DO JCH=ICH+1,NBCHP * * FABRICATION DU VECTEUR ELEMENTAIRE POUR LA * MULTIPLICATION A GAUCHE DO IND=1,NINCD I1=ICPR(NUM(NOELED(IND),IEL)) IF (I1.EQ.0) THEN DD(IND)=0.D0 ELSE I2=IPOSD(IND) DD(IND)=VPO(I1,I2,JCH) ENDIF ENDDO * * * ********************************************* * BOUCLES 5 ET 6 SUR LES LIGNES ET LES COLONNES * DE LA MATRICE ELEMENTAIRE * ********************************************* VA=0.D0 * * MATRICE ELEMENTAIRE SYMETRIQUE * => ON MULTIPLIE PLUTOT PAR LA MATRICE TRANSPOSEE CAR * C'EST PLUS PERFORMANT EN FORTRAN IF (ISYME.EQ.0) THEN DO IND=1,NINCD IF (ABS(DD(IND)).GT.XPETIT) THEN VB=0.D0 DO INP=1,NINCP VB=VB+RE(INP,IND,IEL)*CC(INP) ENDDO VA=VA+DD(IND)*VB ENDIF ENDDO * * MATRICE ELEMENTAIRE DIAGONALE ELSEIF (ISYME.EQ.3) THEN DO IN=1,NINCP VA=VA+DD(IN)*RE(IN,IN,IEL)*CC(IN) ENDDO * * MATRICE ELEMENTAIRE ANTI-SYMETRIQUE OU QUELCONQUE ELSE DO IND=1,NINCD IF (ABS(DD(IND)).GT.XPETIT) THEN VB=0.D0 DO INP=1,NINCP VB=VB+RE(IND,INP,IEL)*CC(INP) ENDDO VA=VA+DD(IND)*VB ENDIF ENDDO ENDIF * XMATR1.RE(ICH,JCH,1)=XMATR1.RE(ICH,JCH,1)+COER*VA * ENDDO ENDDO 36 CONTINUE * SEGDES,XMATRI,DESCR,MELEME * 30 CONTINUE * SEGDES,MRIGID SEGSUP,TRAV2,TCOMP,IHARM,TICPR * * * ================================= * REMPLISSAGE DU TRIANGLE INFERIEUR * ================================= DO ICH=1,NBCHP DO JCH=ICH+1,NBCHP XMATR1.RE(JCH,ICH,1)=XMATR1.RE(ICH,JCH,1) ENDDO ENDDO * * * * * +---------------------------------------------------------------+ * | | * | C H A P E A U D U M R I G I D | * | | * +---------------------------------------------------------------+ * NRIGEL=1 SEGINI,MRIGID IRIG2=MRIGID MTYMAT='RIGIDITE' COERIG(1)=DFLO IRIGEL(1,1)=IPT1 IRIGEL(2,1)=0 IRIGEL(3,1)=DES1 IRIGEL(4,1)=XMATR1 IRIGEL(5,1)=0 IRIGEL(6,1)=0 IRIGEL(7,1)=ISYM1 xmatr1.symre=isym1 SEGDES,XMATR1 ICHOLE=0 IMGEO1=0 IMGEO2=0 IFORIG=IFOUR ISUPEQ=0 JRCOND=0 JRDEPP=0 JRDEPD=0 JRELIM=0 JRGARD=0 JRTOT=0 IMLAG=0 IPROFO=0 IVECRI=0 SEGDES,MRIGID * * RETURN * END * *