C RAFF SOURCE OF166741 26/04/22 21:15:02 12520 SUBROUTINE RAFF C----------------------------------------------------------------------C C Operateur RAFF C C----------------------------------------------------------------------C C C Syntaxe 1 : GEO2 = RAFF GEO1 CHPO1 ; C C Syntaxe 2 : LRE2 = RAFF LRE1 ENT1 ; C C----------------------------------------------------------------------C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC CCGEOME c-INC CCREEL -INC SMCOORD -INC SMELEME -INC SMCHPOI -INC SMMODEL -INC SMCHAML INTEGER MCOL,NCOL C-------------------------------- C B) VARIABLES RELATIVES A LA 2D ! C-------------------------------- SEGMENT KARETE(NBNDS,NCOL) SEGMENT KARET2(NBNDS,NCOL) SEGMENT KMILIE(NBNDS,NCOL) SEGMENT KARPOS(NBNDS) C C-------------------------------- C C) VARIABLES RELATIVES A LA 3D ! C-------------------------------- SEGMENT IRELA(NBNDS) SEGMENT JPLANS(JPLA1,JPLA2) SEGMENT JPLAN2(JPLA1,JPLA2) SEGMENT JPLAN3(JPLA1,JPLA2) SEGMENT JPLCOM(JPLA1) SEGMENT JNOEFA(JNBFA,5) SEGMENT IWORK2(JNBSOM) SEGMENT JFARAF(JNBFA,LLLL) SEGMENT IWORK1(IJKLMN) C-------------------------- C A) VARIABLES 'GENERALES' ! C-------------------------- SEGMENT ICPR(nbpts,2) SEGMENT LTYPE(NBLIS,2) SEGMENT LTYPE2(NBLIS2,4) SEGMENT LTYPE3(NBLIS3,2) SEGMENT LTYPE4(NBLIS3,2) SEGMENT IDENSI REAL*8 DENSI(NBNDS) ENDSEGMENT SEGMENT XDEN(nbpts) XPETI= 1E-30 C C--------------------------------------------------- C D) VARIABLES RELATIVES AUX MAILLAGES DE RELATIONS ! C--------------------------------------------------- POINTEUR IPT10.MELEME POINTEUR IPT11.MELEME POINTEUR IPT12.MELEME POINTEUR IPT13.MELEME POINTEUR IPT14.MELEME POINTEUR IPT15.MELEME POINTEUR IPT16.MELEME SEGMENT MAIREL(NBMAIL) C======================================================================= C Lecture des donnees utilisateurs C======================================================================= C C SYNTAXE 2 : cas du LISTREEL C --------------------------- CALL LIROBJ('LISTREEL',ILRE1,0,IRETOU) IF (IERR.NE.0) RETURN C IF (IRETOU.NE.0) THEN ICAS = 0 CALL LIRENT(KR1,0,IRETOU) IF (IERR.NE.0) RETURN IF (IRETOU.EQ.0) THEN CALL LIRREE(XR1,1,IRETOU) IF (IERR.NE.0) RETURN ICAS = 2 ELSE ICAS = 1 ENDIF CALL RAFLRE(ILRE1,KR1,XR1,ICAS,ILRE2) IF (IERR.NE.0) RETURN IF (ILRE2.NE.0) THEN CALL ECROBJ('LISTREEL',ILRE2) ELSE CALL ERREUR(5) ENDIF RETURN ENDIF C SYNTAXE 1 : cas du MAILLAGE C --------------------------- C On lit le maillage a raffiner et le champoint de densite a affecter c ipt1 : maillage initial a raffiner c mchpo1 : chpoint de densite voulu (qui devrait etre defini sur ipt1) CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) CALL ACTOBJ('MAILLAGE',IPT1,1) IF (IERR.NE.0) RETURN CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU) CALL ACTOBJ('CHPOINT',MCHPO1,1) IF (IERR.NE.0) RETURN C======================================================================= C Initialisations - Declarations C======================================================================= NCOL=0 MCOL=0 NBOUCL=0 NPTINI=nbpts IJKLMN=4 SEGINI IWORK1 IWORK1(1)=4 IWORK1(2)=8 IWORK1(3)=6 IWORK1(4)=10 SEGINI,XDEN C======================================================================= C Tests sur les donnees C======================================================================= C On teste pour savoir si le maillage IPT1 (maillage 1) a un type C d'elements implemente. c TRI3 TRI6 QUA4 QUA8 CUB8 CUB26 PRY6 PR15 TET4 TET10 SURE C les pyramides 5 et 13 ont été débranchés c creation du tableau NBSOU1 = IPT1.LISOUS(/1) NBLIS = MAX(1,NBSOU1) NSURE=0 SEGINI LTYPE IPT2=IPT1 DO 1 INB=1,NBLIS IF (NBSOU1.NE.0) IPT2=IPT1.LISOUS(INB) NT=IPT2.ITYPEL IF (NT.EQ.48) NSURE=NSURE+1 LTYPE(INB,1)=INB LTYPE(INB,2)=NT IF (NT.EQ.16) THEN write (*,*) 'Attention le raffinement PRI6 a debrancher !' ENDIF IF (NT.EQ.23) THEN write (*,*) 'Attention le raffinement TET4 a debrancher !' ENDIF IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND. + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND. + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48)) THEN WRITE(*,*) 'Erreur, type d element non pris en compte' RETURN ENDIF c write(*,*) 'type d elem du maillage de depart', LTYPE(INB,1), c + LTYPE(INB,2) C IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND. C + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND. C + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48).AND.(NT.NE.25).AND. C + (NT.NE.26)) RETURN IF ((NBSOU1.EQ.0).AND.(NT.EQ.48)) RETURN 1 CONTINUE C======================================================================= C Separation du maillage initial C======================================================================= C Si le maillage initial comporte des lisous de type 48, on le separe C en deux maillages dont un ne sera compose que des lisous de type 48. C Maillage d'elements de type 48 : IPT6 C Maillage sans elements de type 48 : IPT5 C IPT1 reste le maillage initial complet LCONF2=0 IF (NSURE.EQ.0) GOTO 7 LCONF2=1 7 CONTINUE C------------------------------- C A) Creation du maillage IPT6 ! C------------------------------- C A.1/** Cas ou il n'y a qu'un lisous de type 48 IF (NSURE.EQ.1) THEN DO 2 K=1,LTYPE(/1) IF (LTYPE(K,2).EQ.48) IPT6=IPT1.LISOUS(LTYPE(K,1)) 2 CONTINUE ENDIF C C A.2/** Cas ou il y a plusieurs lisous de type 48 IF (NSURE.GT.1) THEN NBREF=0 NBSOUS=NSURE NBNN=0 NBELEM=0 SEGINI IPT6 KLISOU=0 DO 3 K=1,LTYPE(/1) IF (LTYPE(K,2).EQ.48) THEN KLISOU=KLISOU+1 IPT6.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1)) ENDIF 3 CONTINUE ENDIF IF (NSURE.EQ.0) THEN NBREF=0 NBSOUS=NSURE NBNN=0 NBELEM=0 SEGINI IPT6 ENDIF C C------------------------------ C B) Creation de IPT5 ! C------------------------------ C NBSOUS=LTYPE(/1)-NSURE NBREF=0 NBNN=0 NBELEM=0 C B.0/** Cas ou il n'y a qu'un lisous en tout IF ((NBSOUS.EQ.1).AND.(NSURE.EQ.0)) THEN IPT5=IPT1 ENDIF C B.1/** Cas ou il n'y a qu'un lisous de type different de 48 IF ((NBSOUS.EQ.1).AND.(NSURE.NE.0)) THEN DO 4 K=1,LTYPE(/1) IF (LTYPE(K,2).NE.48) IPT5=IPT1.LISOUS(LTYPE(K,1)) 4 CONTINUE ENDIF C C B.2/** Cas ou il y a plusieurs lisous de types differents de 48 IF (NBSOUS.GT.1) THEN SEGINI IPT5 KLISOU=0 DO 5 K=1,LTYPE(/1) IF (LTYPE(K,2).NE.48) THEN KLISOU=KLISOU+1 IPT5.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1)) ENDIF 5 CONTINUE ENDIF c compteur de relaton globale LCONF3=LCONF2 C======================================================================= C Debut de la boucle sur les itérations de Raffinement C======================================================================= 9 CONTINUE c write(*,*) 'RAFF: ITERATION boucle 9',NBOUCL c write(*,*) 'ipt5 : ', ipt5 c write(*,*) 'ipt6 : ', ipt6 c write(*,*) 'lconf3 : ', lconf3 C C======================================================================= C Renumerotation locale des noeuds du maillage C======================================================================= C ICPR(I,1)=J : Le Ieme noeud global est le Jeme local C ICPR(I,2)=J : Le Ieme noeud global appartient a J elements c* SEGACT IPT5 NBOUCL=NBOUCL+1 SEGINI ICPR NBNDS=0 NCOL=0 NBS5 = IPT5.LISOUS(/1) IPT2=IPT5 DO 11 INB=1,MAX(1,NBS5) IF (NBS5.NE.0) THEN IPT2=IPT5.LISOUS(INB) c* SEGACT IPT2 ENDIF IF (IPT2.ITYPEL.EQ.48) GOTO 11 DO 10 J=1,IPT2.NUM(/2) DO 110 I=1,IPT2.NUM(/1) ITMP=IPT2.NUM(I,J) IF (ICPR(ITMP,1).NE.0) THEN ICPR(ITMP,2)=ICPR(ITMP,2)+1 NCOL=MAX(NCOL,ICPR(ITMP,2)) ELSE NBNDS=NBNDS+1 ICPR(ITMP,1)=NBNDS ICPR(ITMP,2)=1 NCOL=MAX(NCOL,ICPR(ITMP,2)) ENDIF 110 CONTINUE 10 CONTINUE 11 CONTINUE IF (IDIM.EQ.2) NCOL=NCOL*2 IF (IDIM.EQ.3) NCOL=NCOL*4 C======================================================================= C Prise en compte des noeuds de relations C======================================================================= C IRELA(I,1) = 1 si le noeud I (en numerotation locale) est un noeud C de relation ("hanging node") C IRELA(I,1) = 0 sinon SEGINI IRELA IF (LCONF3.EQ.0) GOTO 1200 SEGACT IPT6 DO 1218 LAA=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT2=IPT6 ELSE IPT2=IPT6.LISOUS(LAA) ENDIF SEGACT IPT2 IF (IPT2.ITYPEL.NE.48) GOTO 1218 DO 1220 J=1,IPT2.NUM(/2) IRELA(ICPR(IPT2.NUM(1,J),1))=1 1220 CONTINUE 1218 CONTINUE 1200 CONTINUE C C======================================================================= C Prise en compte du champoint de densite C======================================================================= C On affecte une densite a chacun des noeuds du maillage. On calcule C cette densite a partir des densites existantes dans le champoint, a C l'aide de la fonction RAFDEN. C Remarque : on ne prend pas en compte les densites qui pourraient C exister dans le tableau XCOOR. C------------------------------------------ C A) Recuperation des valeurs du champoint ! C------------------------------------------ C XDEN(I,1) = densite affectee au noeud I (en numerotation globale) IF (NBOUCL.EQ.1) THEN SEGACT MCHPO1 IF (MCHPO1.IPCHP(/1).NE.1) RETURN MSOUPO=MCHPO1.IPCHP(1) SEGACT MSOUPO IPT3=IGEOC MPOVA1=IPOVAL SEGACT IPT3 SEGACT MPOVA1 IF (MPOVA1.VPOCHA(/2).NE.1) RETURN NBPTCH=IPT3.NUM(/2) c SEGINI ICHPO DO 13 I=1,NBPTCH J=IPT3.NUM(1,I) c XR(I)=XCOOR((J-1)*(IDIM+1)+1) c YR(I)=XCOOR((J-1)*(IDIM+1)+2) c IF (IDIM.EQ.3) ZR(I)=XCOOR((J-1)*(IDIM+1)+3) c DEN(I)=MPOVA1.VPOCHA(I,1) c if(J.eq.115.or.J.eq.321.or. c . J.eq.6261.or.J.eq.5815) then C write(*,*) 'RAFF: #',J,'DEN(',I,'/',NBPTCH,')=',DEN(I) c endif XDEN(J)=MPOVA1.VPOCHA(I,1) 13 CONTINUE C Rem ;seulement pour les noeud "peres", C pour les "fils" raff2d remplit XDEN ENDIF C C-------------------------------------------------------------- C B) Affectation d'une densite a chacun des noeuds du maillage ! C-------------------------------------------------------------- C IDENSI(I,1) = densite affectee au noeud I (en numerotation locale) SEGINI IDENSI DO 15 INB=1,MAX(1,IPT5.LISOUS(/1)) IF (IPT5.LISOUS(/1).EQ.0) THEN IPT2=IPT5 ELSE IPT2=IPT5.LISOUS(INB) SEGACT IPT2 ENDIF c write(*,*) 'ipt2 : ', ipt2, 'type', IPT2.ITYPEL c write(*,*) IPT2.NUM(/2), IPT2.NUM(/1) c IF (IPT2.ITYPEL.EQ.48) GOTO 15 DO 14 J=1,IPT2.NUM(/2) DO 114 I=1,IPT2.NUM(/1) NGLOB=IPT2.NUM(I,J) NLOC=ICPR(NGLOB,1) c XPT=XCOOR((NGLOB-1)*(IDIM+1)+1) c YPT=XCOOR((NGLOB-1)*(IDIM+1)+2) c ZPT=0 c IF (IDIM.EQ.3) ZPT=XCOOR((NGLOB-1)*(IDIM+1)+3) c CALL RAFDEN(ICHPO,XPT,YPT,ZPT,IDIM,DENPT) c DENSI(NLOC)=DENPT DENSI(NLOC)=XDEN(NGLOB) c WRITE (*,*) 'densite au pt', NGLOB, XDEN(NGLOB), XPETI IF (XDEN(NGLOB).LT.XPETI) THEN WRITE (*,*) 'Champ de densité en entrée nul au point', + NGLOB, 'impossible de raffiner' RETURN ENDIF 114 CONTINUE 14 CONTINUE 15 CONTINUE C C======================================================================= C Remplissage du tableau KARETE C======================================================================= C On remplit le tableau KARETE avec les donnees provenant de l'analyse C du maillage. C KARETE(i,j)=k : le ieme noeud local forme une arete C avec le keme noeud local (k>i) C KARPOS(i)=k : compteur de position (k=nb d'aretes touchant i) C KARET2(i,j)=k : l'arete i-j appartient a k elements JPLA1=NBNDS JPLA2=NCOL SEGINI KARETE,KARPOS,KARET2 NBARET=0 NCOL=1 DO 18 INB=1,MAX(1,IPT5.LISOUS(/1)) IF (IPT5.LISOUS(/1).EQ.0) THEN IPT2=IPT5 ELSE IPT2=IPT5.LISOUS(INB) SEGACT IPT2 ENDIF c IF (IPT2.ITYPEL.EQ.48) GOTO 18 ILPL=LPL(IPT2.ITYPEL) ILPT=LPT(IPT2.ITYPEL) DO 17 J=1,IPT2.NUM(/2) DO 117 K=1,ILPL*2-1,2 NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J) NPTB=IPT2.NUM(KSEGM(ILPT+K),J) NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) NEXIST=0 DO 16 I=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I 16 CONTINUE c création d'une nouvelle arrète IF (NEXIST.EQ.0) THEN IF ((KARPOS(NMIN)+1).GT.NCOL) THEN NCOL=KARPOS(NMIN)+1 SEGADJ KARETE, KARET2 ENDIF KARPOS(NMIN)=KARPOS(NMIN)+1 KARETE(NMIN,KARPOS(NMIN))=NMAX KARET2(NMIN,KARPOS(NMIN))=1 NBARET=NBARET+1 c NCOL=MAX(NCOL,KARPOS(NMIN)) ENDIF IF (NEXIST.NE.0) THEN KARET2(NMIN,NEXIST)=KARET2(NMIN,NEXIST)+1 ENDIF 117 CONTINUE 17 CONTINUE 18 CONTINUE c SEGADJ KARETE,KARET2 SEGDES KARETE,KARPOS,KARET2 C C======================================================================= C Remplissage des tableaux JPLANS et JNOEFA C======================================================================= C On remplit les tableaux de description des faces dans le cas 3D. C JPLANS(i,j)=k : la jeme face dont le plus petit sommet en C numerotation locale est le sommet i est la face k C JNOEFA(i,1-3)=j-k-l : Les 3 plus petits sommets en numerotation C locale de la face i sont les noeuds j,k et l C JNOEFA(i,4-5)=j-k : La face identifiee par les 3 sommets precedents C appartient aux lisous de type j et k C JPLCOM(i)=k : compteur de position (k=nb de faces touchant i) IF (IDIM.NE.3) GOTO 40 C---------------------------------------------------- C A) Comptage du nombre d'elements total du maillage ! C---------------------------------------------------- NELTOT=0 DO 50 INB=1,MAX(1,IPT5.LISOUS(/1)) IF (IPT5.LISOUS(/1).EQ.0) THEN IPT2=IPT5 ELSE IPT2=IPT5.LISOUS(INB) SEGACT IPT2 ENDIF NELTOT=NELTOT+IPT2.NUM(/2) 50 CONTINUE C C------------------------------------------------- C B) Initialisation des tableaux JPLANS et JNOEFA ! C------------------------------------------------- JNBFA=NELTOT*6 JMAXPL=0 JNUMFA=0 SEGINI JPLANS,JPLCOM,JNOEFA C C----------------------------------------------- C C) Remplissage des tableaux JPLANS et JNOEFA ! C----------------------------------------------- DO 51 INB=1,MAX(1,IPT5.LISOUS(/1)) IF (IPT5.LISOUS(/1).EQ.0) THEN IPT2=IPT5 ELSE IPT2=IPT5.LISOUS(INB) c* SEGACT IPT2 ENDIF JTYP2=IPT2.ITYPEL C IF (JTYP2.EQ.48) GOTO 51 DO 52 KELEM=1,IPT2.NUM(/2) JLTEL1=LTEL(1,JTYP2) JLTEL2=LTEL(2,JTYP2) DO 53 IBB=1,JLTEL1 JLDEL1=LDEL(1,JLTEL2-1+IBB) JTYFAC=IWORK1(JLDEL1) JLDEL2=LDEL(2,JLTEL2-1+IBB) JNBSOM=NBSOM(JTYFAC) JSPOS=NSPOS(JTYFAC) SEGINI IWORK2 DO 54 IAA=1,JNBSOM NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM) IWORK2(IAA)=NGLOBA 54 CONTINUE C C---------------------------------------------------- C D) Tri des sommets de la face par ordre croissant ! C---------------------------------------------------- C On va ranger les 3 plus petits sommets de la face (en numerotation C globale) dans les variables NPTA, NPTB, NPTC. On les classe par ordre C croissant. C C 1/ Classement des 3 sommets par ordre croissant (num globale) NPTA=nbpts+1 NPTB=NPTA+1 NPTC=NPTB+1 DO 55 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 55 CONTINUE C C 2/ Classement des 3 sommets par ordre croissant (num 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---------------------------------------- C E) Renseignement des nouvelles faces ! C---------------------------------------- SEGACT LTYPE c IHERED=0 c DO 1022 MOI=1,IPT2.NUM(/1) c IHERED=IHERED+IRELA(ICPR(IPT2.NUM(MOI,KELEM),1)) c 1022 CONTINUE C Si NPTA Ne correspond a aucune face : On cree une face IF (JPLCOM(NPTA).EQ.0) THEN JNUMFA=JNUMFA+1 JPLCOM(NPTA)=JPLCOM(NPTA)+1 JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA JNOEFA(JNUMFA,1)=NPTA JNOEFA(JNUMFA,2)=NPTB JNOEFA(JNUMFA,3)=NPTC JNOEFA(JNUMFA,4)=LTYPE(INB,2) c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2) C On signale les sous face de faces supportant un sure C JNOEFA(JNUMFA,5)= 48 IF (LCONF3.NE.0) THEN SEGACT IPT6 DO 501 LAA1=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT4=IPT6 ELSE IPT4=IPT6.LISOUS(LAA1) ENDIF SEGACT IPT4 IF (IPT4.ITYPEL.NE.48) GOTO 501 IF (IPT4.ICOLOR(1).NE.3) GOTO 501 DO 502,MAA1=1,IPT4.NUM(/2) NEXI1=0 NEXI2=0 DO 503 ICC1=1,JNBSOM IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN NEXI1=NEXI1+1 ENDIF 503 CONTINUE IF (NEXI1.GT.0) THEN DO 504 ICC2=1,JNBSOM DO 1504 INN2=2,IPT4.NUM(/1) IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN NEXI2=NEXI2+1 ENDIF 1504 CONTINUE 504 CONTINUE IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48 ENDIF 502 CONTINUE 501 CONTINUE ENDIF C Si NPTA correspond a aumoins une face ELSE NEXIST=0 DO 56 IDD=1,JPLCOM(NPTA) ITEMP1=JPLANS(NPTA,IDD) NAA=JNOEFA(ITEMP1,1) NBB=JNOEFA(ITEMP1,2) NCC=JNOEFA(ITEMP1,3) C Si la face que actuelle existe déja : on renseigne JNOEFA(NEXIST,5) IF (NAA.EQ.NPTA.AND.NBB.EQ.NPTB.AND.NCC.EQ.NPTC) THEN NEXIST=ITEMP1 JNOEFA(NEXIST,5)=LTYPE(INB,2) ENDIF 56 CONTINUE C Si la face que actuelle n'existe pas : On cree une face IF (NEXIST.EQ.0) THEN JNUMFA=JNUMFA+1 JPLCOM(NPTA)=JPLCOM(NPTA)+1 JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA JNOEFA(JNUMFA,1)=NPTA JNOEFA(JNUMFA,2)=NPTB JNOEFA(JNUMFA,3)=NPTC JNOEFA(JNUMFA,4)=LTYPE(INB,2) c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2) C On signale les sous face de faces supportant un sure C JNOEFA(JNUMFA,5)= 48 IF (LCONF3.NE.0) THEN SEGACT IPT6 DO 505 LAA1=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT4=IPT6 ELSE IPT4=IPT6.LISOUS(LAA1) ENDIF SEGACT IPT4 IF (IPT4.ITYPEL.NE.48) GOTO 505 IF (IPT4.ICOLOR(1).NE.3) GOTO 505 DO 506,MAA1=1,IPT4.NUM(/2) NEXI1=0 NEXI2=0 DO 507 ICC1=1,JNBSOM IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN NEXI1=NEXI1+1 ENDIF 507 CONTINUE IF (NEXI1.GT.0) THEN DO 508 ICC2=1,JNBSOM DO 1508 INN2=2,IPT4.NUM(/1) IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN NEXI2=NEXI2+1 ENDIF 1508 CONTINUE 508 CONTINUE IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48 ENDIF 506 CONTINUE 505 CONTINUE ENDIF ENDIF ENDIF JMAXPL=MAX(JMAXPL,JPLCOM(NPTA)) 53 CONTINUE 52 CONTINUE 51 CONTINUE JPLA2=JMAXPL JNBFA=JNUMFA SEGADJ JPLANS,JPLCOM,JNOEFA SEGDES LTYPE C C======================================================================= C Pre-remplissage du tableau KMILIE C======================================================================= C On remplit le tableau lorsqu'il y a des maillages de relations de C type 48 (IPT6, LCONF3 non nul) C KARETE(i,j)=k ET KMILIE(i,j)=n : le noeud global n est situe C au milieu de l'arete [i,k] 40 CONTINUE SEGINI KMILIE SEGACT KARETE*MOD,KARPOS*MOD,KARET2*MOD IF (LCONF3.EQ.0) THEN GOTO 28 ENDIF 22 CONTINUE NBRELA=0 SEGACT IPT6 DO 218 LAA=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT2=IPT6 ELSE IPT2=IPT6.LISOUS(LAA) ENDIF SEGACT IPT2 IF (IPT2.ITYPEL.NE.48) GOTO 218 IF ((IPT2.ICOLOR(1).NE.1).AND.(IPT2.ICOLOR(1).NE.2)) GOTO 218 DO 220 J=1,IPT2.NUM(/2) NPTA=IPT2.NUM(2,J) NPTB=IPT2.NUM(3,J) NPTC=IPT2.NUM(1,J) NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) NEXIST=0 DO 219 I=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I 219 CONTINUE IF (NEXIST.NE.0) then KMILIE(NMIN,NEXIST)=NPTC c WRITE (*,*) 'REPORT ANCIENNE RELA AU NOEUD:', NPTC c write (*,*) 'karret2(nmin nexist)', KARET2(NMIN,NEXIST) KARET2(NMIN,NEXIST)=2 ELSE IF (IDIM.EQ.2) goto 27 IF ((KARPOS(NMIN)+1).GT.NCOL) THEN NCOL=KARPOS(NMIN)+1 SEGADJ KARETE, KARET2, KMILIE ENDIF JPOSI = KARPOS(NMIN) KARPOS(NMIN)=JPOSI+1 KARETE(NMIN,JPOSI+1)=NMAX KARET2(NMIN,JPOSI+1)=2 KMILIE(NMIN,JPOSI+1)=NPTC c WRITE (*,*) 'création de l arrete:', NMIN, NMAX ENDIF C C----------------------------------------------- C Prise en compte des relations "hereditaires" ! C----------------------------------------------- C Cette partie sert a traiter le cas ou un element serait raffine C plusieurs fois alors qu'a sa frontiere les autres elements ne le C seraient pas. Il y aurait donc plusieurs points sur une arete ayant C des relations de conformite et plus seulement le noeud milieu. C c dans l'exenple ci dessous ou les noeuds C et D sont des hanging nodes c L'arete [C B] n'a pas encore ete creee c'est ce qui est fait ici. c c | | c | | c |A C D |B c x---------X----X----X c | | | | c | | | | 27 CONTINUE IF (IDIM.EQ.2) THEN C 1/ Test sur la premiere arete : NMIN-NPTC IF (NEXIST.NE.0) THEN NEXIS2=0 NMIN2=MIN(NMIN,ICPR(NPTC,1)) NMAX2=MAX(NMIN,ICPR(NPTC,1)) DO 24 I=1,MAX(1,KARPOS(NMIN2)) IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I 24 CONTINUE IF (NEXIS2.EQ.0) THEN NBRELA=NBRELA+1 JNBCOL=KARETE(/2) JPOSI=KARPOS(NMIN2) IF (JNBCOL.GT.JPOSI) THEN KARPOS(NMIN2)=JPOSI+1 KARETE(NMIN2,JPOSI+1)=NMAX2 KARET2(NMIN2,JPOSI+1)=2 ELSEIF (JNBCOL.EQ.JPOSI) THEN NCOL=JNBCOL+1 SEGADJ KARETE,KMILIE,KARET2 KARETE(NMIN2,JPOSI+1)=NMAX2 KARPOS(NMIN2)=JPOSI+1 KARET2(NMIN2,JPOSI+1)=2 ENDIF c WRITE (*,*) 'creation de l arrete:', NMIN2, NMAX2 ELSE KARET2(NMIN2,NEXIS2)=2 ENDIF C C 2/ Test sur la deuxieme arete : NMAX-NPTC NEXIS2=0 NMIN2=MIN(NMAX,ICPR(NPTC,1)) NMAX2=MAX(NMAX,ICPR(NPTC,1)) DO 25 I=1,MAX(1,KARPOS(NMIN2)) IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I 25 CONTINUE IF (NEXIS2.EQ.0) THEN NBRELA=NBRELA+1 JNBCOL=KARETE(/2) JPOSI=KARPOS(NMIN2) IF (JNBCOL.GT.JPOSI) THEN KARPOS(NMIN2)=JPOSI+1 KARETE(NMIN2,JPOSI+1)=NMAX2 KARET2(NMIN2,JPOSI+1)=2 ELSEIF (JNBCOL.EQ.JPOSI) THEN NCOL=JNBCOL+1 SEGADJ KARETE,KMILIE,KARET2 KARETE(NMIN2,JPOSI+1)=NMAX2 KARPOS(NMIN2)=JPOSI+1 KARET2(NMIN2,JPOSI+1)=2 ENDIF c WRITE(*,*) 'creation de l arrete:', NMIN2, NMAX2 ELSE KARET2(NMIN2,NEXIS2)=2 ENDIF ENDIF ENDIF 220 CONTINUE 218 CONTINUE IF (NBRELA.NE.0) GOTO 22 28 CONTINUE C======================================================================= C Pre-remplissage du tableau JFARAF C======================================================================= C On remplit le tableau JFARAF et JPLAN3 lorsqu'il y a des maillages de C relations de type 48 (IPT6, LCONF3 non nul). c C Si JPLANS(i,j)=k et JFARAF(k , (l-1)*2+1)= n : le noeud n (en C numérotation globale) est le hanging node au centre de la jieme face C dont le plus petit sommet en numerotation locale est le sommet i. C C Si de plus la face en question est un TRI6 et JFARAF(k , l*2)=m : C le noeud m (en numérotation globale) est le 4ieme noeud de l'element C 48 coul 4 C C Si de plus la face en question est un QUA8-1 et JFARAF(k , l*2)=m : C le noeud m (en numérotation globale) est le 5eme noeud de l'element c 48 coul 4 C C si JPLAN3(i,j)=k il y a k relation de compatibilité de faces associé a C la jieme face dont le plus petit sommet en numerotation locale est le C sommet i. IF (IDIM.NE.3) GOTO 711 JNBFA=JNOEFA(/1) LLLL=10 SEGINI JFARAF SEGINI JPLAN2,JPLAN3 IF (LCONF3.EQ.0) GOTO 711 MPOS25=0 JNBSOM=0 1287 CONTINUE NBRELA=0 DO 670 MIE=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT2=IPT6 ELSE IPT2=IPT6.LISOUS(MIE) ENDIF SEGACT IPT2 IF (IPT2.ITYPEL.NE.48) GOTO 670 IF ((IPT2.ICOLOR(1).EQ.1).OR.(IPT2.ICOLOR(1).EQ.2)) GOTO 670 LIGNUM=IPT2.NUM(/1) IF (IPT2.ICOLOR(1).EQ.4) MPOS25=4 IF (IPT2.ICOLOR(1).EQ.5) MPOS25=5 IF (IPT2.ICOLOR(1).EQ.3) JNBSOM=4 IF (IPT2.ICOLOR(1).EQ.4) JNBSOM=3 IF (IPT2.ICOLOR(1).EQ.5) JNBSOM=4 IF (IPT2.ICOLOR(1).EQ.6) JNBSOM=4 SEGINI IWORK2 C C------------------------- C A) Recherche de la face ! C------------------------- C 1/ Classement des 3 sommets par ordre croissant (num globale) DO 671,MAA=1,IPT2.NUM(/2) DO 673 MBB=1,JNBSOM IWORK2(MBB)=IPT2.NUM(LIGNUM-JNBSOM+MBB,MAA) 673 CONTINUE NPTA=nbpts+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 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 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 c 3.1/ si la face en question n'existe pas, dans le cas général on c s'arrète IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).NE.3)) THEN c WRITE (*,*) 'IPT2.ICOLOR(1)', IPT2.ICOLOR(1) WRITE (*,*) 'Raffinement impossible, le champ de densite + presente des variations trop brusques' WRITE (*,*) 'sure concerne', MAA WRITE (*,*) 'premier nd du sure concerne', IPT2.NUM(1,MAA) WRITE (*,*) 'numero du 1er nd de la face concernee', NPTA KARAF2=0 c GOTO 679 ENDIF c 3.2/ si la face en question n'existe pas dans le cas d'une relation c de face QUA4 il suffie de la créer IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).EQ.3)) THEN JPLCOM(NPTA)= JPLCOM(NPTA)+1 JNUMFA=JNUMFA+1 JNUMFS=JNUMFA JNUMF1=JPLCOM(NPTA) JMAXPL=MAX(JMAXPL,JPLCOM(NPTA)) c 3.2.1/ ajustement des segments JPLA2=JMAXPL JNBFA=JNUMFA SEGADJ JPLANS SEGADJ JNOEFA SEGADJ JPLAN3 SEGADJ JFARAF C 3.2.1 Renseignement des iformation de la nouvelle face dans JPLANS et c JNOEFA, initialisation de JPLAN3 JPLANS(NPTA,JNUMF1)=JNUMFA JPLAN3(NPTA,JNUMF1)=0 JNOEFA(JNUMFA,1)= NPTA JNOEFA(JNUMFA,2)= NPTB JNOEFA(JNUMFA,3)= NPTC JNOEFA(JNUMFA,4)= 48 JNOEFA(JNUMFA,5)= 48 ENDIF IF (NEXIS2.NE.0) THEN JNUMFS=JPLANS(NPTA,NEXIS2) ENDIF C C 4/ Recherche de la position de la face dans JPLANS NEXIS2=0 DO 675 MAC=1,JPLCOM(NPTA) IF (JPLANS(NPTA,MAC).EQ.JNUMFS) NEXIS2=MAC 675 CONTINUE C C------------------------------------- C B) Renseignement de JFARAF et JPLAN3! C------------------------------------- JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 LTEMPO=(JPLAN3(NPTA,NEXIS2)-1)*2+1 JFARAF(JNUMFS,LTEMPO)=IPT2.NUM(1,MAA) IF (MPOS25.NE.0) JFARAF(JNUMFS,LTEMPO+1)=IPT2.NUM(MPOS25,MAA) IF (JNOEFA(JNUMFS,5).EQ.0) JNOEFA(JNUMFS,5)=48 671 CONTINUE 670 CONTINUE C IF (NBRELA.NE.0) GOTO 1287 711 CONTINUE C======================================================================= C Calcul du critere et raffinement si necessaire C======================================================================= C Pour calculer le critere, on a besoin d'un modele et d'un chamelem C constant egal a 1. C C------------------------------------------ C A) Creation d'un modele simple : MMODE1 ! C------------------------------------------ MN3 =12 NFOR=1 NMAT=2 NPARMO=0 NOBMOD=0 N1=1 SEGINI MMODE1 SEGINI IMODE1 MMODE1.KMODEL(1)=IMODE1 IMODE1.FORMOD(1)='MECANIQUE' IMODE1.MATMOD(1)='ELASTIQUE' IMODE1.MATMOD(2)='ISOTROPE' IMODE1.CONMOD =' ' KARAF2=0 C======================================================================= C Boucle sur les eventuels sous-domaines de IPT5 C======================================================================= c----------------------------------------------------------------------- c initialisation de LTYPE3 qui va contenir les tableau LTYP2 c correspondant à chaque sous domaine de IPT5. Ces tableau LTYP2 c contiendront eux même les maillages raffinés issue des c différents sous domaine de IPT5. c LTYPE3(I,1)=LTYP2 pour le sous domaine I de IPT5 c LTYPE3(I,2)=IRR le permier indice d'un maillage 48 dans LTYP2 c----------------------------------------------------------------------- NBSOU5=IPT5.LISOUS(/1) NBLIS3=MAX(1,NBSOU5) SEGINI LTYPE3 MCONF=0 NSUR = 0 DO 33 INB=1,NBLIS3 IF (NBSOU5.EQ.0) THEN IPT2=IPT5 ELSE IPT2=IPT5.LISOUS(INB) ENDIF IMODE1.IMAMOD=IPT2 IMODE1.NEFMOD=IPT2.ITYPEL CALL PRQUOI(IMODE1) c PRQUOI n'active que le dernier MINTE de INFMOD qui n'est pas forcement c celui utilise par INTGCA. Appel a ACTOBJ pour activer tous les MINTE. CALL ACTOBJ('MMODEL ',MMODE1,1) C------------------------------------------------------- C B) Creation d'un chamelem constant egal a 1 : MCHEL1 ! C------------------------------------------------------- CALL ZEROP(MMODE1,'NOEUD ',MCHEL1) c* SEGACT MCHEL1 MCHAM1=MCHEL1.ICHAML(1) SEGACT MCHAM1 MELVA1=MCHAM1.IELVAL(1) SEGACT MELVA1*MOD MELVA1.VELCHE(1,1)=1.D0 C C------------------------------------------------------------- C C) Calcul de la surface ou du volume des elements : MCHEL2 ! C------------------------------------------------------------- ichm2 = 0 CALL INTGCA(MMODE1,MCHEL1,ichm2,1,IRET,XRET,MCHEL2) IF (IRET.NE.1) RETURN IF (XRET.LT.XPETI) THEN WRITE (*,*) 'Volume(resp surf) de la zone a raffiner nul' RETURN ENDIF c* SEGACT IPT2 C--------------------------------------------------------------------- C D) Creation d'un chamelem indiquant les elements a raffiner MCHAM2 ! C--------------------------------------------------------------------- SEGACT KMILIE*MOD SEGACT MCHEL2 MCHAM2=MCHEL2.ICHAML(1) SEGACT MCHAM2 MELVA2=MCHAM2.IELVAL(1) SEGACT MELVA2*MOD NACREE=0 KARAF=0 XDIM=1.D0/IDIM NBNN2=IPT2.NUM(/1) JTYP2=IPT2.ITYPEL c boucles sur les éléments de ipt2 DO 32 J=1,IPT2.NUM(/2) XMOY=0.D0 DO 29 I=1,NBNN2 NGLOB=IPT2.NUM(I,J) NLOC=ICPR(NGLOB,1) DENPT=DENSI(NLOC) XMOY=XMOY+DENPT/NBNN2 29 CONTINUE XINT=MELVA2.VELCHE(1,J) C Un element est a raffiner si la moyenne des densites affectees a C l'ensemble de ses noeuds est inférieure a la racine carree de sa C surface (en 2D) ou a la racine cubique de son volume (en 3D) C WRITE (*,*) 'XMOY', XMOY, '(XINT**XDIM)',(XINT**XDIM) IF ((XMOY-(XINT**XDIM)).LT.0) THEN c c y a t'il des +epsilon...? c XDIFF = XMOY-(XINT**XDIM) c if(XDIFF.ge.0.and.XDIFF.lt.0.001*XMOY) then c write(*,*) 'RAFF : XMOY,XDIFF=',XMOY,XDIFF, c & 'pour l EF',J,(IPT2.NUM(iou,J),iou=1,IPT2.NUM(/1)) c endif c IF (XDIFF.LT.0) THEN MELVA2.VELCHE(1,J)=1 c KARAF nombre d'éléments à rafiner dans le sous domaine KARAF=KARAF+1 ILPT=LPT(JTYP2) ILPL=LPL(JTYP2) c boucle sur les segments de l'élément à raffiner DO 31 K=1,ILPL*2-1,2 NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J) NPTB=IPT2.NUM(KSEGM(ILPT+K),J) NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) NEXIST=0 DO 30 I=1,MAX(1,KARPOS(NMIN)) IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I 30 CONTINUE c On met dans kmilie les arêtes où il faudra mettre un nouveau noeud c On leur associe la valeur -1 IF (KMILIE(NMIN,NEXIST).EQ.0) THEN NACREE=NACREE+1 KMILIE(NMIN,NEXIST)=-1 ENDIF 31 CONTINUE C C------------------------------------------- C E) Comptage du nombre de faces a raffiner ! C------------------------------------------- C 1/ Initialisations IF (IDIM.NE.3) GOTO 32 JLTEL1=LTEL(1,JTYP2) JLTEL2=LTEL(2,JTYP2) DO 73 IBB=1,JLTEL1 JLDEL1=LDEL(1,JLTEL2-1+IBB) JTYFAC=IWORK1(JLDEL1) JLDEL2=LDEL(2,JLTEL2-1+IBB) C JNBPTS=NBNNE(JTYFAC) JNBSOM=NBSOM(JTYFAC) JSPOS=NSPOS(JTYFAC) SEGINI IWORK2 DO 74 IAA=1,JNBSOM NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),J) IWORK2(IAA)=NGLOBA 74 CONTINUE C C 2/ Classement des 3 sommets par ordre croissant (num globale) NPTA=nbpts+1 NPTB=NPTA+1 NPTC=NPTB+1 DO 75 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 75 CONTINUE C C 3/ Classement des 3 sommets par ordre croissant (num 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 4/ Calcul du nombre de points a creer pour raffiner la face NEXIS2=0 DO 76 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 76 CONTINUE IF (NEXIS2.NE.0) THEN KFARAF=JPLAN2(NPTA,NEXIS2) IF (KFARAF.EQ.0) THEN NACREE=NACREE+NBINTE(JTYFAC) JPLAN2(NPTA,NEXIS2)=JPLAN2(NPTA,NEXIS2)+1 ENDIF ENDIF 73 CONTINUE ELSE MELVA2.VELCHE(1,J)=0 ENDIF 32 CONTINUE c KARAF2 nombre d'éléments à rafiner sur tout les sous domaines KARAF2=KARAF2+KARAF c si aucun élément à raffiner on stoque directement IPT2 dans LTYPE3 c on stoque également les anciennes relations et on passe au c sous-domaine suivant. c ATTENTION SI NBOUCL=1 il est possible que des relations surabondantes c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c forcage : nb max itération de raff c IF (NBOUCL.EQ.3) THEN c KARAF2=0 c KARAF=0 c ENDIF c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (KARAF.EQ.0) THEN IF (NBOUCL.EQ.1) THEN IF (LCONF3.GT.0) THEN SEGACT IPT6 NBLIS2=1+MAX(1,IPT6.LISOUS(/1)) SEGINI LTYPE2 LTYPE2(1,1)=1 LTYPE2(1,2)=IPT2.ITYPEL LTYPE2(1,3)=IPT2 LTYPE2(1,4)=IPT2.ICOLOR(1) LTYPE3(INB,2)=0 DO 676 MIE=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT4=IPT6 ELSE IPT4=IPT6.LISOUS(MIE) ENDIF NSUR=1 SEGACT IPT4 LTYPE2(1+MIE,1)=1+MIE LTYPE2(1+MIE,2)=IPT4.ITYPEL LTYPE2(1+MIE,3)=IPT4 LTYPE2(1+MIE,4)=IPT4.ICOLOR(1) LTYPE3(INB,2)=2 676 CONTINUE ELSE NBLIS2=1 SEGINI LTYPE2 LTYPE2(1,1)=1 LTYPE2(1,2)=IPT2.ITYPEL LTYPE2(1,3)=IPT2 LTYPE2(1,4)=IPT2.ICOLOR(1) LTYPE3(INB,2)=0 ENDIF LTYPE3(INB,1)=LTYPE2 ELSE SEGACT LTYPE4 LTYPE3(INB,1)=LTYPE4(INB,1) LTYPE3(INB,2)=LTYPE4(INB,2) ENDIF GOTO 33 ENDIF C C---------------------------- C F) Raffinement du maillage ! C---------------------------- NML=2 C C 1/ Raffinement d'un maillage 2D C On recupere dans MELEME le maillage raffine sans ses relations de C conformite et dans IPT8 les relations de ce maillage C (IPT8 encore incomplet en sortie de RAFF2D) c entrees de raff2D: c - IPT2: maillage élémentaire à raffiner c (ici sous domaine de ipt5) c - ICPR: tableau de passage noeuds local/global c (ici à partir de ipt5 interne à la boucle 9) c - KARPOS: tableau du nb d'arretes par noeuds c (ici à partir de ipt5, mis a jour avec ipt6) c - KARETE: tableau des d'arretes KARETE(i,j)=k c (ici à partir de ipt5, mis a jour avec ipt6) c - KMILIE: tableau des hanging nodes associés au arrètes c en entrée: à partir des relations de ipt6 et c des nouveau noeuds à construires c si relation issue de ipt6 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 (ici à partir de ipt2 et de densi) c - NACREE: nombre de noeuds à créer dans ipt2 c - KARAF: nombre d'élément à raffiner dans ipt2 c - IDCP : ??? variable non utilisé ??? c - XDEN : densité aux noeuds en notation globale c (ici à partir de mchpo1) c - KARET2: tableau du nombre d'éléments touchant les arretes c (ici à partir de ipt5, mis a jour avec ipt6) csorties de raff2D : c - KMILIE: tableau des hanging nodes associés au arrè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 à suprimer 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 - MELEME: Maillage élémentaire raffiné à partir de ipt2 c sans relations c - IPT8 : Maillage de relation incomplet c si J nouvel élément 48: c IPT8.NUM(1,j) = nouveau hanging node c IPT8.NUM(2:4,j) = autres noeuds de la relation c IPT8.NUM(5,j) = 1 si arête a 2 noeuds c 2 si arête a 3 noeuds c si J nouvel élément 48: c IPT8.NUM(1,j) = nouveau hanging node c IPT8.NUM(2:5,j) = 0 c - XDEN : densité aux noeuds en notation globale c interpolé aux nouveaux noeuds. IF (IDIM.EQ.2) CALL RAFF2D(IPT2,ICPR,KARPOS,KARETE,KMILIE, + MELVA2,NACREE,KARAF,MELEME,IDCP,IPT8,KARET2,XDEN) C C 2/ Raffinement d'un maillage 3D IF (IDIM.EQ.3) THEN NML=0 IF ((IPT2.ITYPEL.EQ.25).OR.(IPT2.ITYPEL.EQ.26)) NML=1 ENDIF C a) Maillage 3D contenant des pyramides IPT9=0 IF (NML.EQ.1) CALL RAFPYR(IPT2,ICPR,KARPOS,KARETE,KMILIE, + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA, + IPT8,JFARAF,KARET2,IPT9,XDEN) C b) Maillage 3D sans pyramides IF (NML.EQ.0) CALL RAFF3D(IPT2,IPT6,ICPR,KARPOS,KARETE,KMILIE, + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA, + IPT8,JFARAF,KARET2,XDEN) C C======================================================================= C Creation du maillage de relations de conformite C======================================================================= C-------------------- C A) Initialisations ! C-------------------- IPT10=0 IPT11=0 IPT12=0 IPT13=0 IPT14=0 IPT15=0 IPT16=0 MCONF=0 c IF (KARAF2.EQ.0) GOTO 33 IF (IPT8.NE.0) THEN C C--------------------------------------- C B) Traitement des anciennes relations ! C--------------------------------------- C A ce stade, IPT8 contient les relations crees lors du dernier C raffinement (noeud de relation + noeuds formant la relation) et les C anciennes relations a reporter (noeud de relation seulement) C 1/ Recherche de l'existence d'anciennes relations a reporter NB10=0 MLONG=IPT8.NUM(/1) DO 775 MM=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM).EQ.0) NB10=NB10+1 775 CONTINUE IF (NB10.NE.0) THEN C C 2/ Creation de IPT10 c MCOL=0 c NBNN=9 c NBELEM=NB10 c NBREF=0 c NBSOUS=0 c SEGINI IPT10 c IPT10.ITYPEL=48 C C 3/ Recherche des anciennes relations a reporter C Cette partie permet de compléter IPT8 avec les noeuds formant les C relations, pour les anciennes relations à reporter DO 628 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.0) GOTO 628 INRELA=IPT8.NUM(1,MM11) c On active pous le cas ou IPT6=IPT2. SEGACT IPT6 DO 625 KOU=1,MAX(1,IPT6.LISOUS(/1)) IF (IPT6.LISOUS(/1).EQ.0) THEN IPT2=IPT6 ELSE IPT2=IPT6.LISOUS(KOU) ENDIF SEGACT IPT2 IF (IPT2.ITYPEL.NE.48) GOTO 625 DO 626 JKG=1,IPT2.NUM(/2) INOEGL=IPT2.NUM(1,JKG) IF (INOEGL.NE.INRELA) GOTO 626 MCOL=MCOL+1 DO 627 MTMP=1,IPT2.NUM(/1) IPT8.NUM(MTMP,MM11)=IPT2.NUM(MTMP,JKG) IPT8.NUM(MLONG,MM11)= IPT2.ICOLOR(JKG) 627 CONTINUE GOTO 628 626 CONTINUE 625 CONTINUE 628 CONTINUE ENDIF c à ce stade le maillage de relation ipt8 est complet C C--------------------------------------- C C) Traitement des nouvelles relations ! C--------------------------------------- C Ici, IPT8 contient toutes les relations du nouveau maillage C On les trie par type de relations C 1/ Recherche de l'existence de nouvelles relations NB11=0 NB12=0 NB13=0 NB14=0 NB15=0 NB16=0 MLONG=IPT8.NUM(/1) DO 776 MM=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM).EQ.1) NB11=NB11+1 IF (IPT8.NUM(MLONG,MM).EQ.2) NB12=NB12+1 IF (IPT8.NUM(MLONG,MM).EQ.3) NB13=NB13+1 IF (IPT8.NUM(MLONG,MM).EQ.4) NB14=NB14+1 IF (IPT8.NUM(MLONG,MM).EQ.5) NB15=NB15+1 IF (IPT8.NUM(MLONG,MM).EQ.6) NB16=NB16+1 776 CONTINUE C C 2/ Cas des relations 'aretes a 2 noeuds' IF (NB11.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=3 NBELEM=NB11 NBREF=0 NBSOUS=0 SEGINI IPT11 IPT11.ITYPEL=48 DO 778 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.1) GOTO 778 MCOL=MCOL+1 DO 777 MTMP=1,3 IPT11.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT11.ICOLOR(MCOL)=1 777 CONTINUE 778 CONTINUE ENDIF C C 3/ Cas des relations 'aretes a 3 noeuds' IF (NB12.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=4 NBELEM=NB12 NBREF=0 NBSOUS=0 SEGINI IPT12 IPT12.ITYPEL=48 DO 668 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.2) GOTO 668 MCOL=MCOL+1 DO 667 MTMP=1,4 IPT12.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT12.ICOLOR(MCOL)=2 667 CONTINUE 668 CONTINUE ENDIF C C 4/ Cas des relations 'face QUA4' IF (NB13.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=5 NBELEM=NB13 NBREF=0 NBSOUS=0 SEGINI IPT13 IPT13.ITYPEL=48 DO 658 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.3) GOTO 658 MCOL=MCOL+1 DO 657 MTMP=1,5 IPT13.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT13.ICOLOR(MCOL)=3 657 CONTINUE 658 CONTINUE ENDIF C C 5/ Cas des relations 'face TRI6' IF (NB14.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=7 NBELEM=NB14 NBREF=0 NBSOUS=0 SEGINI IPT14 IPT14.ITYPEL=48 DO 648 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.4) GOTO 648 MCOL=MCOL+1 DO 647 MTMP=1,7 IPT14.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT14.ICOLOR(MCOL)=4 647 CONTINUE 648 CONTINUE ENDIF C C 6/ Cas des relations 'face QUA8-1' IF (NB15.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=9 NBELEM=NB15 NBREF=0 NBSOUS=0 SEGINI IPT15 IPT15.ITYPEL=48 DO 638 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.5) GOTO 638 MCOL=MCOL+1 DO 637 MTMP=1,9 IPT15.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT15.ICOLOR(MCOL)=5 637 CONTINUE 638 CONTINUE ENDIF C C 7/ Cas des relations 'face QUA8-2' IF (NB16.NE.0) THEN MCONF=MCONF+1 MCOL=0 NBNN=9 NBELEM=NB16 NBREF=0 NBSOUS=0 SEGINI IPT16 IPT16.ITYPEL=48 DO 678 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.6) GOTO 678 MCOL=MCOL+1 DO 677 MTMP=1,9 IPT16.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11) IPT16.ICOLOR(MCOL)=6 677 CONTINUE 678 CONTINUE ENDIF C ENDIF C C-------------------------------------- C D) Creation du maillage de relations ! C-------------------------------------- C 1/ Initialisation de MAIREL C MAIREL contient tous les maillages de relations, certains étant vides C car il n'y a aucune relation du type considéré dans le nveau maillage NBMAIL=6 SEGINI MAIREL MAIREL(1)=IPT11 MAIREL(2)=IPT12 MAIREL(3)=IPT13 MAIREL(4)=IPT14 MAIREL(5)=IPT15 MAIREL(6)=IPT16 C C 2/ Initialisation de LTYPE2 C LTYPE2 va contenir tous les sous-objets du nouveau maillage, C relations et non relations NBLIS2=MCONF+MAX(1,MELEME.LISOUS(/1)) SEGINI LTYPE2 C C 3/ Renseignement des lisous non 48 dans LTYPE2 C S'il y a N sous-objets de non-relations dans le nouveau maillage, C les N-ièmes premières colonnes de LTYPE2 contiennent ces sous-objets DO 685 IZZ=1,MAX(1,MELEME.LISOUS(/1)) IF (MELEME.LISOUS(/1).EQ.0) THEN IPT2=MELEME ELSE IPT2=MELEME.LISOUS(IZZ) SEGACT IPT2 ENDIF LTYPE2(IZZ,1)=IZZ LTYPE2(IZZ,2)=IPT2.ITYPEL LTYPE2(IZZ,3)=IPT2 LTYPE2(IZZ,4)=0 685 CONTINUE C C 4/ Renseignement des lisous 48 dans LTYPE2 C Les colonnes suivantes contiennent les sous-objets de relations, C chaque sous-objet correspondant a un type de relations de conformité IZZ=NBLIS2-MCONF DO 686 IXX=1,6 IF (MAIREL(IXX).EQ.0) GOTO 686 IPT2=MAIREL(IXX) IZZ=IZZ+1 LTYPE2(IZZ,1)=IZZ LTYPE2(IZZ,2)=IPT2.ITYPEL LTYPE2(IZZ,3)=IPT2 LTYPE2(IZZ,4)=IPT2.ICOLOR(1) 686 CONTINUE C-------------------------------------- C E) Renseignement de LTYPE2 interne a c la boucle 33 dans LTYPE3 externe a c la boucle 33 C--------------------------------------C C LTYPE3(INB,1)=LTYPE2 IF (MCONF.EQ.0) THEN LTYPE3(INB,2)=0 ELSE LTYPE3(INB,2)=NBLIS2-MCONF+1 ENDIF c fin de boucle sur les sous domaines d'ipt5 33 CONTINUE c SEGSUP IPT2 C Mise a jour des maillages en entrée sortie de la boucle 9 C IPT5, IPT6, IPT7 C======================================================================= 679 CONTINUE C--------------------------------------C C A) initialisation C--------------------------------------C c 1/ initialisaton de IPT7 maillage stockant tout les sous maillage c à la fin d'une itération de la boucle 9 NBS7=0 DO 331 I3=1,LTYPE3(/1) LTYPE2=LTYPE3(I3, 1) SEGACT LTYPE2 NBS7=NBS7+LTYPE2(/1)+1 331 CONTINUE IF (NBS7.EQ.1) NBS7=0 NBNN=0 NBELEM=0 NBREF=0 NBSOUS=NBS7 SEGINI IPT7 c 2/ initialisaton de IPT6 maillage stockant tout les sous maillage c de relation à la fin d'une itération de la boucle 9 NBS6=0 DO 332 I3=1,LTYPE3(/1) LTYPE2=LTYPE3(I3, 1) SEGACT LTYPE2 IF (LTYPE3(I3,2).NE.0) NBS6=NBS6+(LTYPE2(/1)-LTYPE3(I3,2)+1) 332 CONTINUE IF (NBS6.EQ.1) NBS6=0 NBNN=0 NBELEM=0 NBREF=0 NBSOUS=NBS6 SEGINI IPT6 c 3/ initialisaton de IPT5 maillage stockant tout les sous maillage c non relation à la fin d'une itération de la boucle 9 NBS5=0 DO 333 I3=1,LTYPE3(/1) LTYPE2=LTYPE3(I3, 1) SEGACT LTYPE2 IF (LTYPE3(I3,2).EQ.0) THEN NBS5=NBS5+LTYPE2(/1) ELSE NBS5=NBS5+(LTYPE3(I3,2)-1) ENDIF 333 CONTINUE IF (NBS5.EQ.1) NBS5=0 NBNN=0 NBELEM=0 NBREF=0 NBSOUS=NBS5 SEGINI IPT5 C--------------------------------------C C B) assemblage des nouveaux maillages a c partir des sous maillages stockes c dans LTYP3 C--------------------------------------C IND5=0 IND6=0 IND7=0 c 1/assemblage des nouveaux maillages non relation a partir des sous c maillages stockes dans LTYP3 DO 35 I3=1,LTYPE3(/1) c extraction d'un LTYPE2 de LTYPE3 LTYPE2=LTYPE3(I3, 1) SEGACT LTYPE2 c NBNN=0 c NBELEM=0 c NBREF=0 c NBSOUS=0 DO 613 I2=1,MAX(1,(LTYPE3(I3,2)-1)) c extraction d'un MAILLAGE non relation de LTYPE2 IPT2=LTYPE2(I2,3) SEGACT IPT2 IND5=IND5+1 IF (NBS5.EQ.0) THEN IPT5=IPT2 c write (*,*) 'IPT5.' , IPT5 c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ELSE IPT5.LISOUS(IND5)=IPT2 c write (*,*) 'IPT5.LISOUS(',IND5,')' , IPT5.LISOUS(IND5) c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ENDIF IND7=IND7+1 IF (NBS7.EQ.0) THEN IPT7=IPT2 ELSE IPT7.LISOUS(IND7)=IPT2 c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7) c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ENDIF 613 CONTINUE 35 CONTINUE c 2/assemblage des nouveaux maillage de relation à partir des derniers c sous maillages stockés dans LTYPE3 DO 36 I3=1,LTYPE3(/1) c extraction d'un LTYPE2 de LTYPE3 c write (*,*) 'I3', I3 c write (*,*) 'LTYPE3(I3,2)', LTYPE3(I3,2) IF (LTYPE3(I3,2).EQ.0) GOTO 36 SEGACT LTYPE2 LTYPE2=LTYPE3(I3, 1) NBNN=0 NBELEM=0 NBREF=0 NBSOUS=0 DO 614 I2=LTYPE3(I3,2), LTYPE2(/1) c extraction d'un MAILLAGE de relation de LTYPE2 c write (*,*) 'I2', I2 SEGINI IPT2 IPT2=LTYPE2(I2,3) SEGACT IPT2 c write (*,*) 'IPT2: Type', IPT2.ITYPEL,'nbelem', IPT2.NUM(/2) c write (*,*) 'coul IPT2' , IPT2.ICOLOR(1) IND6=IND6+1 IF (NBS6.EQ.0) THEN IPT6=IPT2 c write (*,*) 'IPT6' , IPT6 c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ELSE IPT6.LISOUS(IND6)=IPT2 c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6) c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ENDIF 614 CONTINUE 36 CONTINUE c 3/ Concatenation des sous maillages de relation de même couleur c GG 11 juin 2018 : attention il peut encore y avoir des problemes a la c jonctions entre deux sous zones dans le cas ou on raffine plusieurs c types d'element differents (passer cette etape de concatenation dans c la boucle 33 ? IF (NBS6.EQ.0) then IND7=IND7+1 IPT7.LISOUS(IND7) = IPT6 c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7) c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) goto 370 ENDIF NBS60 = 0 SEGACT IPT6 c boucle sur les couleurs NBELEM = 0 NBNN= 0 NBREF= 0 NBSOUS=6 SEGINI IPT4 DO 37 ICOUL = 1, 6 c boucle sur les sous maillages de ipt6 NCOl = 0 DO 38 I6 = 1, NBS6 IPT3 = IPT6.LISOUS(I6) SEGACT IPT3 IF ((IPT3.icolor(1)).EQ.ICOUL) then NCOL = NCOL+1 IF (NCOL.EQ.1) THEN SEGINI IPT2 IPT2 = IPT6.LISOUS(I6) ENDIF IF (NCOL.GE.2) THEN c IPT3 = IPT6.LISOUS(I6) c SEGACT IPT3 c Plusieurs sous maillages de même couleur on concatene IPT2 et IPT3 NBELEM = IPT2.NUM(/2) NBNN= IPT2.NUM(/1) NBREF=0 NBSOUS=0 NBELE0 = NBELEM DO 381 IE3 = 1, IPT3.NUM(/2) DO 382 IE2 = 1, NBELE0 IF ((IPT3.NUM(1,IE3)).EQ.(IPT2.NUM(1,IE2))) GOTO 381 382 CONTINUE NBELEM = NBELEM +1 SEGADJ IPT2 DO IN2 = 1, NBNN IPT2.NUM(IN2 , NBELEM) = IPT3.NUM(IN2, IE3) IPT2.icolor(NBELEM) = IPT3.icolor(IE3) ENDDO 381 CONTINUE ENDIF ENDIF 38 CONTINUE IF (NCOL.GE.1) THEN IPT4.LISOUS(ICOUL) = IPT2 NBS60= NBS60+1 ENDIF 37 CONTINUE c mise a jour de ipt6 et ipt7 IF (NBS60.EQ.1) then IPT6 = IPT2 segact ipt6 c write (*,*) 'IPT6' , IPT6 c write (*,*) 'type', IPT6.ITYPEL, 'Nbel' , IPT6.num(/2) IND7=IND7+1 IPT7.LISOUS(IND7) = IPT2 c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7) c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2) ELSE SEGACT IPT4 NBELEM = 0 NBNN = 0 NBREF = 0 NBSOUS = NBS60 SEGINI IPT6 IS6=0 DO ICOUL = 1, 6 IF ((IPT4.LISOUS(ICOUL)).NE.0) THEN IS6 = IS6+1 IPT6.LISOUS(IS6)= IPT4.LISOUS(ICOUL) c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6) c write (*,*) 'type', IPT4.ITYPEL, 'Nbel' , IPT4.num(/2) IND7=IND7+1 IPT7.LISOUS(IND7) = IPT4.LISOUS(ICOUL) ENDIF ENDDO ENDIF 370 CONTINUE NBELEM = 0 NBNN = 0 NBREF = 0 NBSOUS = IND7 SEGADJ IPT7 SEGINI LTYPE4 LTYPE4 = LTYPE3 c mise a jour du nombre global de relations LCONF3 = MCONF C======================================================================= C Test sur la progression du raffinement C======================================================================= C C CB215821 Retrait du caractere *MOD le cas echeant SEGACT IPT5,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE c SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI c SEGSUP ICPR,KARETE,KARPOS,KARET2,IDENSI SEGSUP ICPR,KARETE,KARPOS,KARET2 SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1 SEGSUP MMODE1,MCHEL1,MCHEL2 IF ((IDIM.EQ.3).AND.(NBOUCL.NE.1)) SEGSUP JPLANS,JPLAN2,JPLAN3, +JPLCOM,JNOEFA,IWORK2,JFARAF, IPT9 C S'il y a eu a raffiner des elements lors de l'iteration courante, C on passe a l'iteration suivante C c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c forcage : nb max itération de raff c IF (NBOUCL.EQ.3) THEN c KARAF2=0 c KARAF=0 c ENDIF c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (KARAF2.NE.0) GOTO 9 C ### SORTIE TEMPORAIRE DU PROGRAMME ### C Si on est a la premiere iteration de RAFF et qu'il n'y a pas C d'elements a raffiner, on retourne en sortie le maillage de depart, C donné en entree IF ((KARAF2.EQ.0).AND.(NBOUCL.EQ.1)) THEN WRITE(*,*) 'il n est pas necessaire de raffiner le maillage' MELEME=IPT1 GOTO 998 ENDIF C======================================================================= C Mises a jour diverses C======================================================================= C ** On affecte une densite a chaque noeud que l'on a cree c DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1) c XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1)) c 37 CONTINUE C======================================================================= C Fin du programme C======================================================================= MELEME=IPT7 SEGSUP XDEN 998 CONTINUE CALL ACTOBJ('MAILLAGE',MELEME,1) CALL ECROBJ('MAILLAGE',MELEME) c RETURN END