C RAFF3D SOURCE CB215821 20/11/25 13:38:23 10792 SUBROUTINE RAFF3D(IPT2,IPT3,ICPR,KARPOS,KARETE,KMILIE,MELVA2, .NACREE,KARAF,IPT4,JPLANS,JPLAN3,JPLCOM,JNOEFA,IPT7,JFARAF,KARET2, .XDEN) c entrees de raff3D: c - IPT2: maillage élémentaire à raffiner c - IPT3: maillage de SURE précédant c - ICPR: tableau de passage noeuds local/global c - KARPOS: tableau du nb d'arretes par noeuds c - KARETE: tableau des d'arretes KARETE(i,j)=k c - KMILIE: tableau des hanging nodes associés au arêtes c si ancienne relation c KMILIE(i,j)=n : le noeud global n est situe c au milieu de l'arete [i,k] c si nouvelle relation c KMILIE(i,j)=-1 : il faut construire un noeud c au milieu de l'arete [i,k] c - MELVA2: champ par elem valait 1 pour les élem à raffiner c et 0 pour les autres c - NACREE: nombre de noeuds à créer dans ipt2 c - KARAF: nombre d'élément à raffiner dans ipt2 c - JPLANS : tableau des faces JPLAN(i,j)=k c - JPLAN3 : tableau du nombre de relation de compat par face c - JPLCOM : tableau du nb de faces par noeuds c - JNOEFA : tableau des noeuds et éléments associé au faces c - JFARAF : tableau des hanging nodes associés au faces pour c les anciennes relations. c - XDEN : densité aux noeuds en notation globale c - KARET2: tableau du nombre d'éléments touchant les arêtes c sorties: c - KMILIE: tableau des hanging nodes associés au arêtes c (en sortie: à partir des relations donnée en c entrée et de celles crées c si anciene relation à garder c KMILIE(i,j)=n : le noeud global n est situe c au milieu de l'arete [i,k] c si anciene relation à supprimer c KMILIE(i,j)=0 : c si nouvelle relation c KMILIE(i,j)=n : le noeud n à été construit au c milieu de l'arete [i,k] c - JPLAN3 : tableau du nombre de relation de compat par face c mis a jour c - JFARAF : tableau des hanging nodes associés au faces pour c mis a jour. c - MELEME: Maillage élémentaire raffiné à partir de ipt2 c sans relations c - IPT7 : Maillage de relations créé incomplet c - XDEN : densité aux noeuds en notation globale c interpolé aux nouveaux noeuds. IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC CCGEOME -INC SMELEME -INC SMCHPOI -INC SMMODEL -INC SMCHAML C C====================================================================== C Declarations C====================================================================== SEGMENT ICPR(nbpts,2) SEGMENT XDEN(nbpts) SEGMENT KARETE(NBNDS,NCOL) SEGMENT KARET2(NBNDS,NCOL) SEGMENT KMILIE(NBNDS,NCOL) SEGMENT KARPOS(NBNDS) SEGMENT JPLANS(JPLA1,JPLA2) SEGMENT JPLAN3(JPLA1,JPLA2) SEGMENT JPLCOM(JPLA1) SEGMENT JNOEFA(JNBFA,5) SEGMENT JFARAF(JNBFA,LLLL) SEGMENT NUMNOE(INUMNO) SEGMENT IWORK2(JNBSOM) SEGMENT IWORK1(IJKLMN) SEGMENT IWORK4(JNBSOM) SEGMENT IWORK3(IJKLMN) REAL*8 LON68,LON510,LON47,LON29 C C====================================================================== C Initialisations C====================================================================== SEGACT JPLANS,JPLCOM,JNOEFA*MOD,JPLAN3*MOD, IPT3 SEGACT IPT2,ICPR,KARPOS,KARETE,KMILIE*MOD,MELVA2,JFARAF*MOD IJKLMN=4 ipt7=0 SEGINI IWORK1 IWORK1(1)=4 IWORK1(2)=8 IWORK1(3)=6 IWORK1(4)=10 C JNBFA=JNOEFA(/1) C LLLL=2 C IF (IPT2.ITYPEL.EQ.15) LLLL=10 C IF ((IPT2.ITYPEL.EQ.17).OR.(IPT2.ITYPEL.EQ.24)) LLLL=6 C LLLL=10 C SEGINI JFARAF NCOMPL=LNELM(1,(IPT2.ITYPEL-1)*2+2) IF (NCOMPL.NE.0) GOTO 999 NBNN=NBNNE(LNELM(2,(IPT2.ITYPEL-1)*2+1)) INELM=LNELM(1,(IPT2.ITYPEL-1)*2+1) C C Creation du squelette du maillage resultat NBELEM=IPT2.NUM(/2)-KARAF+INELM*KARAF NBSOUS=0 NBREF=0 SEGINI IPT4 IPT4.ITYPEL=LNELM(2,(IPT2.ITYPEL-1)*2+1) LPO2=LPOS2(IPT2.ITYPEL) segact mcoord*mod NBPT0=nbpts NBPT1=nbpts NBPTS=NBPT1+NACREE+KARAF*NBINTE(IPT2.ITYPEL) SEGADJ MCOORD,XDEN INUMNO=NBRAF(IPT2.ITYPEL) SEGINI NUMNOE LCOMP=1 C C====================================================================== C Phase de raffinement 3D C====================================================================== C======================================= C A) Boucle sur les elements a raffiner ! C======================================= DO 6 IARAF=1,MELVA2.VELCHE(/2) IF (MELVA2.VELCHE(1,IARAF).NE.1) THEN DO 1 IJKL=1,IPT2.NUM(/1) IPT4.NUM(IJKL,NBELEM)=IPT2.NUM(IJKL,IARAF) 1 CONTINUE NBELEM=NBELEM-1 GOTO 6 ENDIF JCOMPT=0 JPOS5=LPOS5(IPT2.ITYPEL) C C========================================================== C B) Boucle sur les noeuds a creer pour raffiner l'element ! C========================================================== DO 4 I=1,NBRAF(IPT2.ITYPEL) JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT LTYPNO=JTYPNO(JPOS5-1+I) c write(*,*) 'LTYPNO', LTYPNO C C------------------------------- C ** B.1 / On est sur une arete ! C------------------------------- IF (LTYPNO.EQ.0) THEN NPTA=IPT2.NUM(LISNOE(JLISN),IARAF) NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF) NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) DO 2 K=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K 2 CONTINUE IF (KMILIE(NMIN,NEXIST).GT.0) THEN c si le noeud existe deja on rempli numnoe avec Kmilie NUMNOE(I)=KMILIE(NMIN,NEXIST) JCOMPT=JCOMPT+JLONG GOTO 4 ELSE c sinon on ajoute un noeud a la fin de numnoe et on le met dans kmilie NBPT1=NBPT1+1 NUMNOE(I)=NBPT1 KMILIE(NMIN,NEXIST)=NBPT1 ENDIF c IF ((NBPT1.eq.10330).or.(NBPT1.eq.10341)) THEN c WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) c WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB c WRITE(*,*) 'NMIN', NMIN, 'NMAX', NMAX ,'NEXI', NEXIST c WRITE(*,*) 'KMILIE(NMIN,NEXIST)', KMILIE(NMIN,NEXIST) c WRITE(*,*) 'KARETE(NMIN,NEXIST)', KARETE(NMIN,NEXIST) c endif C C------------------------------ C ** B.2 / On est sur une face ! C------------------------------ C Il faut identifier cette face (numero, sommets...) ELSEIF ((LTYPNO.GT.0).AND.(LTYPNO.LT.7)) THEN C => a) On initialise toutes les donnees relatives a cette face JTYPEL=IPT2.ITYPEL JLTEL2=LTEL(1,JTYPEL)-1+LTYPNO JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO JLDEL1=LDEL(1,JLTEL2) JTYFAC=IWORK1(JLDEL1) JLDEL2=LDEL(2,JLTEL2) JNBSOM=NBSOM(JTYFAC) JSPOS=NSPOS(JTYFAC) C => b) On identifie les sommets de la face (n° global) SEGINI IWORK2 DO 10 IAA=1,JNBSOM NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF) IWORK2(IAA)=NGLOBA 10 CONTINUE C => c) On classe ces sommets par ordre croissant (NPTA < NPTB < NPTC) NPTA=nbpts+1 NPTB=NPTA+1 NPTC=NPTB+1 DO 11 ICC=1,JNBSOM IF (IWORK2(ICC).LT.NPTA) THEN NPTC=NPTB NPTB=NPTA NPTA=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTB) THEN NPTC=NPTB NPTB=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTC) THEN NPTC=IWORK2(ICC) ENDIF 11 CONTINUE C => d) On passe ces sommets en n° locale NPTA1=NPTA NPTB1=NPTB NPTC1=NPTC NPTA2=ICPR(NPTA,1) NPTB2=ICPR(NPTB,1) NPTC2=ICPR(NPTC,1) IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN NPTA=NPTA2 NPTB=MIN(NPTB2,NPTC2) NPTC=MAX(NPTB2,NPTC2) ENDIF IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN NPTA=NPTB2 NPTB=MIN(NPTA2,NPTC2) NPTC=MAX(NPTA2,NPTC2) ENDIF IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN NPTA=NPTC2 NPTB=MIN(NPTA2,NPTB2) NPTC=MAX(NPTA2,NPTB2) ENDIF C => e) On cherche le numero de la face NEXIS2=0 DO 12 IEE=1,JPLCOM(NPTA) MTMP=JPLANS(NPTA,IEE) JJ1=JNOEFA(MTMP,1) JJ2=JNOEFA(MTMP,2) JJ3=JNOEFA(MTMP,3) IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN NEXIS2=IEE ENDIF 12 CONTINUE JNUMFA=JPLANS(NPTA,NEXIS2) C C--------------------------------- C ** B.3 / Raffinement de la face ! C--------------------------------- KTEST1=JPLAN3(NPTA,NEXIS2) KTEST2=NBINTE(JTYFAC) C => a) Si la face est de type QUA4 (un seul noeud a creer, au milieu) IF (JTYFAC.EQ.8) THEN IF (KTEST1.LT.KTEST2) THEN c a.1) si le noeud milieu n'existe pas on l'ajoute a la fin de numnoe c et on rempli JFARAF et JPLAN3 NBPT1=NBPT1+1 NUMNOE(I)=NBPT1 JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 JFARAF(JNUMFA,1)=NBPT1 ELSEIF (KTEST1.EQ.KTEST2) THEN c a.2) si le noeud milieu existe deja on rempli numnoe grace à JFARAF c et on met celui-ci a zéro NUMNOE(I)=JFARAF(JNUMFA,1) JFARAF(JNUMFA,1)=0 JCOMPT=JCOMPT+JLONG GOTO 4 ELSE WRITE(*,*) 'ERREUR' ENDIF c IF ((NBPT1.eq.10248).or.(NBPT1.eq.10276)) THEN c IF ((NPTA.eq.235).or.(NPTA.eq.236).or.(NPTA.eq.237)) THEN c WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) c WRITE(*,*) 'NPTA', NPTA1, 'NPTB', NPTB1, 'NPTC', NPTC1 c WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB, 'NPTC', NPTC c WRITE(*,*) 'JNUMFA', JNUMFA, 'NEXI', NEXIS2 c WRITE(*,*) 'JPLAN3(NPTA,NEXIS2)', JPLAN3(NPTA,NEXIS2) c WRITE(*,*) 'JFARAF(JNUMFA,1)', JFARAF(JNUMFA,1) c WRITE(*,*) 'JNOEFA(JNUMFA,1)', JNOEFA(JNUMFA,1) c WRITE(*,*) 'JNOEFA(JNUMFA,2)', JNOEFA(JNUMFA,2) c WRITE(*,*) 'JNOEFA(JNUMFA,3)', JNOEFA(JNUMFA,3) c WRITE(*,*) 'JNOEFA(JNUMFA,4)', JNOEFA(JNUMFA,4) c WRITE(*,*) 'JNOEFA(JNUMFA,5)', JNOEFA(JNUMFA,5) c endif ENDIF C => b) Si la face n'est pas de type QUA4 (donc de type TRI6 ou QUA8) IF (JTYFAC.NE.8) THEN C b.1) Si il manque des noeuds dans la face C On ajoute un noeud à la fin de numnoe IF (KTEST1.LT.KTEST2) THEN NBPT1=NBPT1+1 NUMNOE(I)=NBPT1 JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 NEXIS3=0 XCO2=0.25 DO 13 KBB=1,JLONG XCO1=XCOEFF(JPOS1-1+KBB) IF (XCO1.EQ.XCO2) NEXIS3=KBB 13 CONTINUE JFARAF(JNUMFA,2*(KTEST1)+1)=NBPT1 IF (NEXIS3.NE.0) THEN C b.1.1) si le point est sur une arête du triangle interieur C (pour un TRI6) C si le point est sur une arête de la coie interieure C (pour un QUA8) C on renseigne la partie paire de JFARAF qui doit déja C exister ???? NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) JFARAF(JNUMFA,2*(KTEST1)+2)=NGLOB ELSE C b.1.2) si le point est le milieu d'un QUA8 C on met la partie paire de JFARAF à zéro JFARAF(JNUMFA,2*(KTEST1)+2)=0 ENDIF ELSEIF (KTEST1.EQ.KTEST2) THEN C b.2) Si il ne manque pas de noeuds dans la face NEXIS3=0 NEXIS4=0 XCO2=0.25 DO 14 KBB=1,JLONG XCO1=XCOEFF(JPOS1-1+KBB) IF (XCO1.EQ.XCO2) NEXIS3=KBB 14 CONTINUE IF (NEXIS3.NE.0) THEN C b.2.1) si le point est sur une arête du triangle interieur C (pour un TRI6) C si le point est sur une arête de la croie interieure C (pour un QUA8) NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) DO 15 KAA=2,JFARAF(/2),2 IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA 15 CONTINUE c on rempli numnoe grace à JFARAF et on met la partie c associée au noeud en question de celui-ci a 0 NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1) JFARAF(JNUMFA,NEXIS4-1)=0 ELSE C b.2.2) si le point est le milieu d'un QUA8 DO 16 KAA=2,JFARAF(/2),2 IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA 16 CONTINUE c on rempli numnoe grace à JFARAF et on met la partie c associée au noeud en question de celui-ci a 0 NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1) JFARAF(JNUMFA,NEXIS4-1)=0 ENDIF JCOMPT=JCOMPT+JLONG GOTO 4 ENDIF ENDIF C C------------------------------------------------------ C ** B.4 / On est a l'interieur du volume de l'element ! C------------------------------------------------------ ELSEIF (LTYPNO.EQ.7) THEN NBPT1=NBPT1+1 NUMNOE(I)=NBPT1 ENDIF C allocation de plus de mémoire dans xcoord et xden si nécéssaire IF (NBPT1.EQ.NBPTS) THEN NBPTS=NBPTS+200 SEGADJ MCOORD,XDEN ENDIF C C============================== C C) Creation du nouveau point ! C============================== C On continue ici que lorsque l'on doit creer un nouveau point XPT=0. YPT=0. ZPT=0. XDEN1=0.D0 XDENMIN = 1000. DO 3 J=1,JLONG NGLOB=IPT2.NUM(LISNOE(JLISN-1+J),IARAF) XINI=XCOOR((NGLOB-1)*(IDIM+1)+1) YINI=XCOOR((NGLOB-1)*(IDIM+1)+2) ZINI=XCOOR((NGLOB-1)*(IDIM+1)+3) XPT=XPT+XINI*XCOEFF(JPOS1-1+J) YPT=YPT+YINI*XCOEFF(JPOS1-1+J) ZPT=ZPT+ZINI*XCOEFF(JPOS1-1+J) XDEN1=XDEN1+XDEN(NGLOB)*XCOEFF(JPOS1-1+J) IF (XDEN(NGLOB).LT.XDENMIN) XDENMIN = XDEN(NGLOB) c IF (NBPT1.eq.10219) THEN c write (*,*) 'nglob', nglob ,'nbpt1', nbpt1 c write (*,*) 'xpt ', xpt c write (*,*) 'ypt ', ypt c write (*,*) 'zpt ', zpt c write (*,*) 'xden1 ', xden1 c write (*,*) 'xcoeff1 ', XCOEFF(JPOS1-1+J) c endif 3 CONTINUE c seuil de densité a xdenmin IF (XDEN1.LT.XDENMIN) XDEN1 = XDENMIN XCOOR((NBPT1-1)*(IDIM+1)+1)=XPT XCOOR((NBPT1-1)*(IDIM+1)+2)=YPT XCOOR((NBPT1-1)*(IDIM+1)+3)=ZPT XCOOR((NBPT1-1)*(IDIM+1)+4)=XDEN1 XDEN(NBPT1)=XDEN1 JCOMPT=JCOMPT+JLONG C====================================================================== 4 CONTINUE C====================================================================== JPOS4=LPOS4(IPT2.ITYPEL) C C Cas des tetraedres LONDIA=0 JTYPEL=IPT2.ITYPEL IF (JTYPEL.EQ.23) THEN LPT6=NUMNOE(6) LPT5=NUMNOE(5) LPT8=NUMNOE(8) LPT10=NUMNOE(10) XPT6=XCOOR((LPT6-1)*(IDIM+1)+1) XPT5=XCOOR((LPT5-1)*(IDIM+1)+1) XPT8=XCOOR((LPT8-1)*(IDIM+1)+1) XPT10=XCOOR((LPT10-1)*(IDIM+1)+1) YPT6=XCOOR((LPT6-1)*(IDIM+1)+2) YPT5=XCOOR((LPT5-1)*(IDIM+1)+2) YPT8=XCOOR((LPT8-1)*(IDIM+1)+2) YPT10=XCOOR((LPT10-1)*(IDIM+1)+2) ZPT6=XCOOR((LPT6-1)*(IDIM+1)+3) ZPT5=XCOOR((LPT5-1)*(IDIM+1)+3) ZPT8=XCOOR((LPT8-1)*(IDIM+1)+3) ZPT10=XCOOR((LPT10-1)*(IDIM+1)+3) LON68=SQRT((XPT8-XPT6)**2+(YPT8-YPT6)**2+(ZPT8-ZPT6)**2) LON510=SQRT((XPT10-XPT5)**2+(YPT10-YPT5)**2+(ZPT10-ZPT5)**2) IF (LON510.LT.LON68) LONDIA=4*4 ENDIF IF (JTYPEL.EQ.24) THEN LPT2=NUMNOE(2) LPT4=NUMNOE(4) LPT7=NUMNOE(7) LPT9=NUMNOE(9) XPT2=XCOOR((LPT2-1)*(IDIM+1)+1) XPT4=XCOOR((LPT4-1)*(IDIM+1)+1) XPT7=XCOOR((LPT7-1)*(IDIM+1)+1) XPT9=XCOOR((LPT9-1)*(IDIM+1)+1) YPT2=XCOOR((LPT2-1)*(IDIM+1)+2) YPT4=XCOOR((LPT4-1)*(IDIM+1)+2) YPT7=XCOOR((LPT7-1)*(IDIM+1)+2) YPT9=XCOOR((LPT9-1)*(IDIM+1)+2) ZPT2=XCOOR((LPT2-1)*(IDIM+1)+3) ZPT4=XCOOR((LPT4-1)*(IDIM+1)+3) ZPT7=XCOOR((LPT7-1)*(IDIM+1)+3) ZPT9=XCOOR((LPT9-1)*(IDIM+1)+3) LON47=SQRT((XPT7-XPT4)**2+(YPT7-YPT4)**2+(ZPT7-ZPT4)**2) LON29=SQRT((XPT9-XPT2)**2+(YPT9-YPT2)**2+(ZPT9-ZPT2)**2) IF (LON29.LT.LON47) LONDIA=4*10 ENDIF C C=================================== C D) Creation des nouveaux elements ! C=================================== C On remplit la portion de IPT4 relative aux elements crees a partir C de la division de l'element IARAF (indice de boucle 1). C Cette portion de IPT4 contient les colonnes dont la valeur s'etend C de INELM*(LCOMP-1)+1 a INELM*LCOMP. DO 5 J=1,INELM DO 5 I=1,IPT4.NUM(/1) NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I) IF (((JTYPEL.EQ.23).OR.(JTYPEL.EQ.24)).AND.(J.GT.4)) THEN NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I+LONDIA) ENDIF IF (NTEMP.GT.NBNN) THEN IPT4.NUM(I,INELM*(LCOMP-1)+J)=NUMNOE(NTEMP-NBNN) ELSE IPT4.NUM(I,INELM*(LCOMP-1)+J)=IPT2.NUM(NTEMP,IARAF) ENDIF 5 CONTINUE LCOMP=LCOMP+1 6 CONTINUE C NBPT =NBPT1 C SEGADJ MCOORD,XDEN C C======================================================================= C Preparation du maillage de relations C======================================================================= SEGSUP KARET2 NBNDS=KARETE(/1) NCOL=KARETE(/2) SEGINI KARET2 C================================== C A) Relations dues aux faces (3D) ! C================================== C Tous les noeuds qui restent dans JFARAF sont a creer en tant que C relations de conformite ou des bords du domaine C----------------------------------------------------------------------- C ** A.1 / Prise en compte des relation C d'arretes pour les arretes de chaque faces C----------------------------------------------------------------------- C----------------------------------------------------------------------- C** A.1.1/ faces raffinée C----------------------------------------------------------------------- c write (*,*) 'relations dues au faces' JTYPEL=IPT2.ITYPEL JLTEL1=LTEL(1,JTYPEL) JLTEL2=LTEL(2,JTYPEL) IF (JTYPEL.NE.14) GOTO 141 c Boucle sur les elements de l'ancien maillage DO 121 KELEM1=1,IPT2.NUM(/2) c Boucle sur les faces de l'ancien maillage DO 122 IBB1=1,JLTEL1 JLDEL1=LDEL(1,JLTEL2-1+IBB1) JTYFAC=IWORK1(JLDEL1) JLDEL2=LDEL(2,JLTEL2-1+IBB1) JNBSOM=NBSOM(JTYFAC) JSPOS=NSPOS(JTYFAC) SEGINI IWORK2 C C 1/ trouver le numéro de la face DO 123 IAA=1,JNBSOM NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM1) IWORK2(IAA)=NGLOBA 123 CONTINUE C C 1.1/ Classement des 3 sommets par ordre croissant de n° globale NPTA=NBPT +1 NPTB=NPTA+1 NPTC=NPTB+1 DO 124 ICC=1,JNBSOM IF (IWORK2(ICC).LT.NPTA) THEN NPTC=NPTB NPTB=NPTA NPTA=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTB) THEN NPTC=NPTB NPTB=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTC) THEN NPTC=IWORK2(ICC) ENDIF 124 CONTINUE C C1.2/ Classement des 3 sommets par ordre croissant de n° locale NPTA2=ICPR(NPTA,1) NPTB2=ICPR(NPTB,1) NPTC2=ICPR(NPTC,1) IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN NPTA=NPTA2 NPTB=MIN(NPTB2,NPTC2) NPTC=MAX(NPTB2,NPTC2) ENDIF IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN NPTA=NPTB2 NPTB=MIN(NPTA2,NPTC2) NPTC=MAX(NPTA2,NPTC2) ENDIF IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN NPTA=NPTC2 NPTB=MIN(NPTA2,NPTB2) NPTC=MAX(NPTA2,NPTB2) ENDIF C C 1.3/ Recherche du numero de la face NEXIS2=0 DO 125 IEE=1,JPLCOM(NPTA) MTMP=JPLANS(NPTA,IEE) JJ1=JNOEFA(MTMP,1) JJ2=JNOEFA(MTMP,2) JJ3=JNOEFA(MTMP,3) IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN NEXIS2=IEE ENDIF 125 CONTINUE JNUMF1=JPLANS(NPTA,NEXIS2) C C 2/ Exclusion des faces non voulues C on exlu les faces où un noeud existait déja IF (JFARAF(JNUMF1,1).EQ.0) GOTO 122 C on exlu les faces de bords IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 122 C 3/ Recherche des arretes de la faces pour y mettre des relations IF (JTYPEL.NE.14) GOTO 137 NBAR1=0 DO 138 ISS=1, (JNBSOM-1) DO 138 JSS=(ISS+1), JNBSOM NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1)) NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1)) DO 139 IRR=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,IRR).EQ.NMAX) THEN NBAR1=NBAR1+1 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1 c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR) c IF (JFARAF(JNUMF1,1).EQ.15748) THEN c write (*,*) ISS,JSS c write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!' c write (*,*) IWORK2(ISS), IWORK2(JSS) c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR) c write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) c ENDIF ENDIF 139 CONTINUE 138 CONTINUE c WRITE (*,*) 'Nbarrete', NBAR1 c IF (NBAR1.NE.4) THEN c write (*,*) IWORK2(1), IWORK2(2) c write (*,*) IWORK2(3), IWORK2(4) c ENDIF 137 CONTINUE 122 CONTINUE 121 CONTINUE 141 CONTINUE C----------------------------------------------------------------------- C** A.1.2/ relations de faces conservées dans le maillage précédent C----------------------------------------------------------------------- SEGINI IPT5 DO 670 MIE=1,MAX(1,IPT3.LISOUS(/1)) IF (IPT3.LISOUS(/1).EQ.0) THEN IPT5=IPT3 ELSE IPT5=IPT3.LISOUS(MIE) ENDIF SEGACT IPT5 c write (*,*) 'ipt3 raff 3D' , ipt3, IPT3.ITYPEL, IPT3.ICOLOR(1) c write (*,*) 'ipt5 raff 3D' , ipt5, IPT5.ITYPEL, IPT5.ICOLOR(1) IF (IPT5.ITYPEL.NE.48) GOTO 670 IF ((IPT5.ICOLOR(1).NE.3)) GOTO 670 LIGNUM=IPT2.NUM(/1) JNBSOM=4 SEGINI IWORK2 C C 1/ trouver le numéro de la face C 1.1/ Classement des 3 sommets par ordre croissant (num globale) DO 671,MAA=1,IPT5.NUM(/2) DO 673 MBB=1,JNBSOM c write (*,*) 'MAA' , MAA, 'MBB', MBB IWORK2(MBB)=IPT5.NUM(MBB+1,MAA) 673 CONTINUE NPTA=nbpt +1 NPTB=NPTA+1 NPTC=NPTB+1 DO 672 ICC=1,JNBSOM IF (IWORK2(ICC).LT.NPTA) THEN NPTC=NPTB NPTB=NPTA NPTA=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTB) THEN NPTC=NPTB NPTB=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTC) THEN NPTC=IWORK2(ICC) ENDIF 672 CONTINUE C C 1.2/ Classement des 3 sommets par ordre croissant (num locale) NPTA1=NPTA NPTB1=NPTB NPTC1=NPTC NPTA2=ICPR(NPTA,1) NPTB2=ICPR(NPTB,1) NPTC2=ICPR(NPTC,1) IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN NPTA=NPTA2 NPTB=MIN(NPTB2,NPTC2) NPTC=MAX(NPTB2,NPTC2) ENDIF IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN NPTA=NPTB2 NPTB=MIN(NPTA2,NPTC2) NPTC=MAX(NPTA2,NPTC2) ENDIF IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN NPTA=NPTC2 NPTB=MIN(NPTA2,NPTB2) NPTC=MAX(NPTA2,NPTB2) ENDIF C C 1.3/ Recherche du numero de la face NEXIS2=0 DO 674 IEE=1,JPLCOM(NPTA) MTMP=JPLANS(NPTA,IEE) JJ1=JNOEFA(MTMP,1) JJ2=JNOEFA(MTMP,2) JJ3=JNOEFA(MTMP,3) IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN NEXIS2=IEE ENDIF 674 CONTINUE JNUMF1=JPLANS(NPTA,NEXIS2) c IF (MAA.EQ.5) THEN c write (*,*) 'MAA', MAA, 'nexis2', nexis2 c write (*,*) IWORK2(1), IWORK2(2), IWORK2(3), IWORK2(4) c write (*,*) 'JFARAF(JNUMF1,1)', JFARAF(JNUMF1,1) c write (*,*) 'JNOEFA(JNUMF1,5)',JNOEFA(JNUMF1,5) c ENDIF C C 2/ Exclusion des faces non voulues IF (NEXIS2.EQ.0) THEN c WRITE (*,*)'FACE NON TROUVE SURE:', MAA GOTO 671 ENDIF C on exlu les faces où un noeud existait déja IF (JFARAF(JNUMF1,1).EQ.0) GOTO 671 C on exlu les faces de bords IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 671 NBAR1=0 C C 3/ Recherche des arretes de la faces pour y mettre des relations DO 135 ISS=1, (JNBSOM-1) DO 135 JSS=(ISS+1), JNBSOM NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1)) NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1)) DO 136 IRR=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,IRR).EQ.NMAX) THEN NBAR1=NBAR1+1 KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1 c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR) c IF (MAA.EQ.5) THEN c write (*,*) ISS,JSS c write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!' c write (*,*) IWORK2(ISS), IWORK2(JSS) c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR) c write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) c ENDIF ENDIF 136 CONTINUE 135 CONTINUE 671 CONTINUE 670 CONTINUE SEGDES IPT5 C--------------------------------- C ** A.2 / Initialisation de IPT5 ! C--------------------------------- C 1/ Comptage du nombre de noeuds soumis a des relations NBRELA=0 DO 114 IHF=1,JFARAF(/1) C on exlu les faces de bords IF (JNOEFA(IHF,5).EQ.0) GOTO 114 C on exlu les faces où un noeud existait déja IF (JFARAF(IHF,1).EQ.0) GOTO 114 DO 113 JF=1,JFARAF(/2),2 IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1 113 CONTINUE 114 CONTINUE C C 2/ Creation de IPT5 IF (NBRELA.EQ.0) THEN IPT5=0 GOTO 210 ENDIF C IF (IPT2.ITYPEL.EQ.14) NBNN=4 C IF (IPT2.ITYPEL.EQ.15) NBNN=8 C IF (IPT2.ITYPEL.EQ.17) NBNN=5 C IF (IPT2.ITYPEL.EQ.24) NBNN=5 NBNN=10 NBELEM=NBRELA NBSOUS=0 NBREF=0 SEGINI IPT5 IPT5.ITYPEL=48 C C 3/ Renseignement des noeuds support des relations NBRELA=0 DO 116 IPF=1,JFARAF(/1) IF (JNOEFA(IPF,5).EQ.0) GOTO 116 IF (JFARAF(IPF,1).EQ.0) GOTO 116 DO 115 JF=1,JFARAF(/2),2 IF (JFARAF(IPF,JF).EQ.0) GOTO 115 NBRELA=NBRELA+1 IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF) 115 CONTINUE 116 CONTINUE c write (*,*) 'NBRELA FACE ', NBRELA C C----------------------------------------------------- C ** A.3 / Recherche des noeuds formant les relations ! C----------------------------------------------------- C 1/ Boucle sur l'ensemble des noeuds a creer DO 200 IARAF=1,MELVA2.VELCHE(/2) IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200 JCOMPT=0 JPOS5=LPOS5(IPT2.ITYPEL) DO 190 I=1,NBRAF(IPT2.ITYPEL) JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT LTYPNO=JTYPNO(JPOS5-1+I) IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189 C C 2/ Preparation pour trouver le noeud et la face en question JTYPEL=IPT2.ITYPEL JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO JLDEL1=LDEL(1,JLTEL2) JTYFAC=IWORK1(JLDEL1) JLDEL2=LDEL(2,JLTEL2) JNBSOM=NBSOM(JTYFAC) JSPOS=NSPOS(JTYFAC) SEGINI IWORK2 C C 3/ Classement des 3 sommets par ordre croissant de n° globale DO 100 IAA=1,JNBSOM NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF) IWORK2(IAA)=NGLOBA 100 CONTINUE NPTA=nbpt +1 NPTB=NPTA+1 NPTC=NPTB+1 DO 110 ICC=1,JNBSOM IF (IWORK2(ICC).LT.NPTA) THEN NPTC=NPTB NPTB=NPTA NPTA=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTB) THEN NPTC=NPTB NPTB=IWORK2(ICC) ELSEIF (IWORK2(ICC).LT.NPTC) THEN NPTC=IWORK2(ICC) ENDIF 110 CONTINUE C C 4/ Classement des 3 sommets par ordre croissant de n° locale NPTA2=ICPR(NPTA,1) NPTB2=ICPR(NPTB,1) NPTC2=ICPR(NPTC,1) IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN NPTA=NPTA2 NPTB=MIN(NPTB2,NPTC2) NPTC=MAX(NPTB2,NPTC2) ENDIF IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN NPTA=NPTB2 NPTB=MIN(NPTA2,NPTC2) NPTC=MAX(NPTA2,NPTC2) ENDIF IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN NPTA=NPTC2 NPTB=MIN(NPTA2,NPTB2) NPTC=MAX(NPTA2,NPTB2) ENDIF C C 5/ Recherche du numero de la face NEXIS2=0 DO 120 IEE=1,JPLCOM(NPTA) MTMP=JPLANS(NPTA,IEE) JJ1=JNOEFA(MTMP,1) JJ2=JNOEFA(MTMP,2) JJ3=JNOEFA(MTMP,3) IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN NEXIS2=IEE ENDIF 120 CONTINUE JNUMFA=JPLANS(NPTA,NEXIS2) C C 6/ Recherche du numero global du point IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189 IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1) IF (JTYFAC.NE.8) THEN NEXIS3=0 NEXIS4=0 XCO2=0.25 DO 140 KBB=1,JLONG XCO1=XCOEFF(JPOS1-1+KBB) IF (XCO1.EQ.XCO2) NEXIS3=KBB 140 CONTINUE IF (NEXIS3.NE.0) THEN NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) DO 150 KAA=2,JFARAF(/2),2 IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA 150 CONTINUE INOEGL=JFARAF(JNUMFA,NEXIS4-1) ELSE DO 160 KAA=2,JFARAF(/2),2 IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA 160 CONTINUE INOEGL=JFARAF(JNUMFA,NEXIS4-1) ENDIF ENDIF C C------------------------------ C ** A.4 / Remplissage de IPT5 ! C------------------------------ C 1/ Recherche de la position du point dans IPT5 NEXIS5=0 DO 170 IGG=1,IPT5.NUM(/2) IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG 170 CONTINUE IF (NEXIS5.EQ.0) GOTO 189 IF (INOEGL.EQ.139) THEN c write (*,*) INOEGL, IWORK2(1), IWORK2(2), IWORK2(3) ENDIF C C 2/ Renseignement des points formant les relations DO 180 IHH=1,JLONG IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF) 180 CONTINUE IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3 IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4 IF (JLONG.EQ.8) THEN IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5 IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6 ENDIF 189 CONTINUE JCOMPT=JCOMPT+JLONG 190 CONTINUE 200 CONTINUE 210 CONTINUE C C=================================== C B) Relations dues aux aretes (2D) ! C=================================== C On cree un maillage IPT6 contenant tous les noeuds soumis a des C relations C c ILPL=LPL(IPT4.ITYPEL) c ILPT=LPT(IPT4.ITYPEL) c DO 317 J=1,IPT4.NUM(/2) c DO 317 K=1,ILPL*2-1,2 c NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J) c NPTB=IPT4.NUM(KSEGM(ILPT+K),J) c c if ((NPTB.eq.129).and.(NPTA.eq.9087)) then c write(*,*) '!ici!',NBPT0 c endif c IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN c GOTO 317 c ENDIF c NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) c NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) c if ((NMIN.eq.235).and.(NMAX.eq.236)) then c write(*,*) 'nmin', nmin, 'nmax', nmax c write(*,*) 'NPTA',NPTA, 'NPTB', NPTB c endif c NEXIST=0 c DO 316 I=1,MAX(1,KARPOS(NMIN)) c IF (KARETE(NMIN,I).EQ.NMAX) THEN c KARET2(NMIN,I)=KARET2(NMIN,I)+1 c c IF (KMILIE(NMIN,I).gt.0) then c prise en compte des arretes multiples sur le nouveau maillage c dans l'exemple ci dessous où les noeuds C et D sont des hanging nodes c L'arette [C B] n'a pas encore été prise en compte dans karet2 c'est c ce qui est fait ici. c c | | c | | c |A C D |B c x---------X----X----X c | | | | c | | | | c1/ Test sur la premiere arete : NMIN-KMILIE(NMIN,I) c NMILI = ICPR(KMILIE(NMIN,I),1) c NEXIS2=0 c NMIN2=MIN(NMIN,NMILI) c NMAX2=MAX(NMIN,NMILI) c DO 318 I2=1,KMILIE(/2) c IF (KMILIE(NMIN2,I2).GT.0) THEN c KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1 c ENDIF c 318 CONTINUE C C 2/ Test sur la deuxieme arete : NMAX-NPTC c NEXIS2=0 c NMIN2=MIN(NMAX,NMILI) c NMAX2=MAX(NMAX,NMILI) c DO 319 I2=1,KMILIE(/2) c IF (KMILIE(NMIN2,I2).GT.0) THEN c KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1 c ENDIF c 319 CONTINUE c ENDIF c ENDIF c 316 CONTINUE c 317 CONTINUE C 1/ Comptage du nombre de noeuds soumis a des relations NBELEM=0 DO 27 J=1,KMILIE(/2) DO 27 I=1,KMILIE(/1) c IF (KMILIE(I,J).eq.10330) write(*,*) '!kmilie', I, J, '=10330!' c IF (KMILIE(I,J).eq.10330) write(*,*) 'KARET2(I,J)', KARET2(I,J) c IF (KMILIE(I,J).eq.10341) write(*,*) '!kmilie', I, J, '=10341!' c IF (KMILIE(I,J).eq.10341) write(*,*) 'KARET2(I,J)', KARET2(I,J) IF (KARET2(I,J).EQ.0) GOTO 27 IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1 27 CONTINUE c write (*,*) 'NBELEM', NBELEM C C 2/ Creation de IPT6 IPT6=0 IF (NBELEM.EQ.0) GOTO 999 NBNN=5 NBREF=0 NBSOUS=0 SEGINI IPT6 IPT6.ITYPEL=48 C C 3/ Renseignement des noeuds support des relations DO 28 J=1,KMILIE(/2) DO 28 I=1,KMILIE(/1) IF (KARET2(I,J).EQ.0) GOTO 28 IF (KMILIE(I,J).GT.0) THEN NBREF=NBREF+1 IPT6.NUM(1,NBREF)=KMILIE(I,J) ENDIF 28 CONTINUE C WRITE (*,*) 'nombre de noeuds supports de rela seg', NBREF C C 4/ Recherche des noeuds formant les relations DO 24 IARAF=1,MELVA2.VELCHE(/2) IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24 JCOMPT=0 JPOS5=LPOS5(IPT2.ITYPEL) DO 23 I=1,NBRAF(IPT2.ITYPEL) JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT LTYPNO=JTYPNO(JPOS5-1+I) IF (LTYPNO.NE.0) GOTO 22 NPTA=IPT2.NUM(LISNOE(JLISN),IARAF) NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF) NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) DO 29 K=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K 29 CONTINUE IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22 IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22 NEXIS5=0 DO 20 MM=1,IPT6.NUM(/2) INOEGL=KMILIE(NMIN,NEXIST) INRELA=IPT6.NUM(1,MM) IF (INOEGL.EQ.INRELA) NEXIS5=MM 20 CONTINUE C C 5/ Renseignement des noeuds formant les relations DO 21 IHH=1,JLONG IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF) 21 CONTINUE IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1 IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2 22 CONTINUE JCOMPT=JCOMPT+JLONG 23 CONTINUE 24 CONTINUE 444 CONTINUE C C============================================ C C) Creation du maillage de relations final ! C============================================ IF (IPT5.EQ.0) THEN IPT7=IPT6 GOTO 999 ENDIF NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2) C NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1)) NBNN=10 NBREF=0 NBSOUS=0 SEGINI IPT7 IPT7.ITYPEL=48 DO 42 NEO=1,IPT5.NUM(/2) DO 42 MOR=1,10 IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO) IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO) 42 CONTINUE NN5 = IPT5.NUM(/2) DO 43 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2) DO 43 MOR=1,IPT6.NUM(/1) IF (MOR.LT.IPT6.NUM(/1)) THEN IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5) ELSE IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5) ENDIF IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5) 43 CONTINUE SEGDES IPT5, IPT6 C======================================================================= C Fin du programme C======================================================================= 999 CONTINUE c write (*,*) 'IPT2', IPT2 c write (*,*) 'IPT3', IPT3 c write (*,*) 'IPT4', IPT4 c write (*,*) 'IPT5', IPT5 c write (*,*) 'IPT6', IPT6 c write (*,*) 'IPT7', IPT7 RETURN END