C PROFCH SOURCE CB215821 26/08/24 21:17:58 12622 SUBROUTINE PROFCH(MATRAK,IMPR,KPIMP,NIMP) C************************************************************************ C C OBJET : C CE SOUS PROGRAMME PREPARE L'INVERSION DE LA MATRICE DE PRESSION C C KTYPI= 1 OU 5, METHODE DIRECTE C C 1/ IOPTI = 0 PAS DE RENUMEROTATION C IOPTI = 1 RENUMEROTATION DE LA MATRICE C DE PRESSION C 2/ CALCUL DU PROFIL ET DE LA TAILLE C DE LA DEMI-MATRICE C C MATRAK : objet de type MATRAK C Si KIZCL=0 On cree IZL C Sinon rien C IMPR 1 ou 2 niveau d'impressions C KPIMP = 0 pas de pression imposee C KPIMP = 1 pression imposée en un point (le dernier) C KPIMP = 2 pression moyenne imposée C NIMP numero de la pression imposée C C************************************************************************ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*** C-INC SMMATRAKANC C************************************************************************* C C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees C * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange * (points CENTRE ) pour chaque operateur de contrainte * KGEOC SPG pour la totalite des points CENTRE. * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse) * KLEMC Connectivites de l'ensemble des contraintes * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones SEGMENT MATRAK INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP) INTEGER LIZAFM(NBSOUS) INTEGER IKAM0 (NBSOUS) INTEGER IMEM (NBELC) INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC ENDSEGMENT SEGMENT IZAFM REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX) ENDSEGMENT POINTEUR IPMJ.IZAFM,IPMK.IZAFM C******************************************************************* -INC SMMATRK1 -INC SMCOORD -INC SMLENTI POINTEUR IZIPAD.MLENTI -INC SMELEME POINTEUR IZK.MELEME -INC SMTABLE POINTEUR MTABP.MTABLE PARAMETER (IBLSIZ=50000) SEGMENT IZTRAV INTEGER ITAB(NJA) ENDSEGMENT C C IVK = Tableau de travail PARAMETER (MAXIVK=100) INTEGER IVK(MAXIVK) C*** CHARACTER*8 TYPE CHARACTER*4 LMOT(1) DATA LMOT/'IMPR'/ DATA IOPTI/1/ C*** SEGACT MATRAK*MOD IF(KIZCL.NE.0)THEN SEGDES MATRAK RETURN ENDIF KGEO=KGEOC CALL KNBEL(KGEO,NBELC) MELEME=KLEMC C C--- CALNOE renvoie dans IZTKB un tableau contenant pour chaque noeud C du domaine la liste des {l{ments poss{dant ce noeud. Les noeuds C sont rep{r{s dans la num{rotation locale(domaniale), les {l{ments C sont pris dans l'ordre fourni par la num{rotation naturelle. C Il est possible, partant de ce tableau, de construire pour chaque C {l{ment la liste des {l{ments voisins. Ceci est tr}s utile pour C le calcul de la matrice de pression ou Akl.ne.0 si k et l sont C deux {l{ments voisins. C C Cette technique n'est pas (encore) utilis{e si KTYPI = 1 ou 5. C Les maillages trait{s dans ce cas {tant "petits", ce n'est pas C trop grave. Par contre lorsque KTYPI=2,3,4, c'est @ dire pour les C tr}s gros maillages, son emploi est obligatoire. C C Cet appel ne doit etre fait que si on stocke la matrice de C pression (KTPI.LE.5). C CALL KALNO0(MELEME,IZK,MLENTI,IZIPAD,IRET) IF(IRET.EQ.0)GO TO 90 C--- Activation de tous les MELEMEs du domaine SEGACT MELEME,IZK,MLENTI NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 DO 400 NS=1,NBSOUS IF(NBSOUS.EQ.1)IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(NS) SEGACT IPT1 400 CONTINUE C C************************************************************************ C C H O L E S K I C************************************************************************ C NJA=NBELC NJAN=NBELC NJAB=NBELC SEGINI IZL,IZTRAV KIZCL=IZL LONG=0 CALL INITI(NUAN,NBELC,0) CALL INITI(NUNA,NBELC,0) NUAN(1)=1 NUNA(1)=1 C********* ON STOKE LES POINTEURS DES MELEME KI=1 DO 40 NS=1,NBSOUS IF(NBSOUS.EQ.1)IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(NS) NEL=IPT1.NUM(/2) IPOINT=IPT1 CALL INITI(IMEL(KI),NEL,IPOINT) DO 41 I=1,NEL IMJ(KI+I-1)=I 41 CONTINUE KI=KI+NEL 40 CONTINUE C********* IF(IOPTI.EQ.1)GO TO 33 C********* ON NE RENUMEROTE PAS IF(IMPR.GE.2)WRITE(6,*)' PAS DE RENUMEROTATION DE LA PRESSION ' DO 32 I=1,NBELC NUAN(I)=I NUNA(I)=I 32 CONTINUE GO TO 34 C********* OPTIMISATION DE LA NUMEROTAION DE LA MATRICE DE PRESSION 33 CONTINUE IF(IMPR.GE.2)WRITE(6,*)' ON RENUMEROTE LA MATRICE DE PRESSION ' C********* C write(6,*)' SUB PROFCH : IPADL :' C nbp=ipadl(/1) C write(6,1001)(IPADL(II),II=1,NBP) C im2=lect(/1) C im1=izk.num(/1) C write(6,*)' IM1=',IM1 ,' Taille IZK.NUM' C do 1873 i=1,im2 C nbb=lect(i) C write(6,*)' I=->',i,' nb=',lect(i) C write(6,1001)(izk.num(ii,i),ii=1,nbb) C1873 continue KC=1 DO 30 KK=1,NBELC K=NUNA(KK) IF(K.NE.0)GO TO 310 C ON A PLUS D'ELEMENT NOUVELLEMENT NUMEROTE WRITE(6,*)' EST CE LA DISCONTINUITE ',K,KC,KK KC=KC+1 DO 311 KKA=1,NBELC IF(NUAN(KKA).EQ.0)GO TO 312 C SI NUAN(KKA)=0 C'EST QUE L'ELEMENT KKA N'A PAS ETE NUMEROTE 311 CONTINUE WRITE(6,*)' IL YA UN PROBLEME DE NUMEROTATION' CALL ARRET(0) 312 CONTINUE NUAN(KKA)=KC NUNA(KC)=KKA K=NUNA(KC) 310 CONTINUE IPT1=IMEL(K) NP=IPT1.NUM(/1) DO 1931 I=1,NP NU=IPT1.NUM(I,IMJ(K)) C NU=IPT1.NUM(I,K) NU=IZIPAD.LECT(NU) NELV=LECT(NU) DO 31 NN=1,NELV KN=IZK.NUM(NN,NU) IF(NUAN(KN).NE.0)GO TO 31 KC=KC+1 NUAN(KN)=KC NUNA(KC)=KN 31 CONTINUE 1931 CONTINUE 30 CONTINUE DO 35 I=1,NBELC ITAB(I)=IMJ(NUNA(I)) 35 CONTINUE CALL RSETI(IMJ,ITAB,NBELC) DO 36 I=1,NBELC ITAB(I)=IMEL(NUNA(I)) 36 CONTINUE CALL RSETI(IMEL,ITAB,NBELC) SEGSUP IZTRAV C ********************************************** 34 CONTINUE IF(IMPR.GE.2)THEN WRITE(6,1929) 1929 FORMAT(/10X,' SUB PROFIM : RENUMEROTAION ZONE PRESSION ') WRITE(6,*)' CORRESPONDANCE NUMER. ANC. --> NUMER. NOUVELLE ' WRITE(6,1930)(N,NUAN(N),N=1,NBELC) WRITE(6,*)' CORRESPONDANCE NUMER. NOUVELLE --> NUMER. ANC. ' WRITE(6,1930)(N,NUNA(N),N=1,NBELC) 1930 FORMAT(7(1X,'**',I5,3X,I5,1X,'**')) ENDIF C DO 10 KK=1,NBELC K=IMJ(KK) IPT1=IMEL(KK) NP=IPT1.NUM(/1) K0=KK DO 1932 I=1,NP NU=IPT1.NUM(I,K) NU=IZIPAD.LECT(NU) NELV=LECT(NU) DO 11 NN=1,NELV KN=IZK.NUM(NN,NU) KN=NUAN(KN) IF(KN.LT.K0)K0=KN 11 CONTINUE 1932 CONTINUE IF(KK.NE.NIMP.OR.KPIMP.NE.2)THEN KZA(KK)=KK-K0+1 ELSE KZA(KK)=NBELC-1 ENDIF LONG=LONG+KZA(KK) 10 CONTINUE C C La matrice est stockée par blocs de plusieurs lignes. En effet, il faut C l'équivalent de 900 opérations pour activer/désactiver un segment. Si C une ligne de la matrice comprend 50 éléments, on fait 100 opérations C utiles et 900 inutiles lors de la résolution. C La taille des blocs est fixée par le paramètre IBLSIZ. La diagonale est C stockée dans un segment à part. C C Il nous faut constituer un certain nombre de tableaux. Le premier segment C décrit la matrice. Il contient un tableau donnant les numéros des premières C lignes de chaque bloc et un tableau de pointeurs sur des segments C descripteurs de bloc. C C Il faut d'abord déterminer le nombre de blocs. C NBLK=1 IF(LONG-NBELC.GT.IBLSIZ) THEN ILBL=0 DO 12 KK=2,NBELC LA=KZA(KK)-1 IF(LA.GT.IBLSIZ) THEN NBLK=NBLK+2 ILBL=0 GOTO 12 ENDIF ILBL=ILBL+LA IF(ILBL.GT.IBLSIZ) THEN NBLK=NBLK+1 ILBL=LA GOTO 12 ENDIF 12 CONTINUE ENDIF C C On peut allouer le segment descripteur de la matrice. On va pouvoir C remplir des numéros de ligne des débuts de bloc. C SEGINI IDMAT IF(NBLK.EQ.1) THEN NLDBLK(1)=2 NLDBLK(2)=NBELC+1 ELSE IBLK=1 ILBL=0 NLDBLK(1)=2 NLDBLK(NBLK+1)=NBELC+1 DO 13 KK=2,NBELC LA=KZA(KK)-1 IF(LA.GT.IBLSIZ) THEN NLDBLK(IBLK+1)=KK NLDBLK(IBLK+2)=KK+1 IBLK=IBLK+2 ILBL=0 GOTO 13 ENDIF ILBL=ILBL+LA IF(ILBL.GT.IBLSIZ) THEN NLDBLK(IBLK+1)=KK ILBL=LA IBLK=IBLK+1 GOTO 13 ENDIF 13 CONTINUE ENDIF C C On peut maintenant remplir le segment descripteur de chaque bloc. C DO 15 IBLK=1,NBLK NLBLK=NLDBLK(IBLK+1)-NLDBLK(IBLK) SEGINI IDBLK IDESCR(IBLK)=IDBLK ILON=0 IDEBLK(1)=1 DO 16 KK=NLDBLK(IBLK),NLDBLK(IBLK+1)-1 LA=KZA(KK)-1 ILON=ILON+LA IDEBLK(KK-NLDBLK(IBLK)+2)=IDEBLK(KK-NLDBLK(IBLK)+1)+LA 16 CONTINUE SEGDES IDBLK 15 CONTINUE SEGDES IDMAT C C On stocke dans KZA1 la référence au segment IDMAT à partir duquel tout se C retrouve. C C write(6,*)' sub PROFCH on met IDMAT -> KZA1 ',IDMAT KZA1=IDMAT 14 CONTINUE C IF(IMPR.GE.2)WRITE(6,*)(KZA(M),M=1,NBELC),LONG C IF(IMPR.GE.1)WRITE(6,200)LONG SEGDES IZL SEGSUP IZK,IZIPAD,MLENTI C--- D{sactivation de tous les MELEMEs DO 51 NS=1,NBSOUS IF(NBSOUS.EQ.1)IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(NS) SEGDES IPT1 51 CONTINUE SEGDES MELEME RETURN 90 CONTINUE WRITE(6,*)' Retour anormal de PROFCH' RETURN C************************************************************************ 100 FORMAT(40X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7, *' MOTS',13X,1H*/40X,1H*,72X,1H*) 200 FORMAT(10X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7, *' MOTS',13X,1H*) 1001 FORMAT(20(1X,I5)) END