C RAFF SOURCE CHAT 11/03/16 21:30:07 6902 SUBROUTINE RAFF IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) -INC CCOPTIO -INC SMCOORD -INC CCGEOME -INC SMELEME -INC SMCHPOI -INC SMMODEL -INC SMCHAML C 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((XCOOR(/1)/(IDIM+1)),2) SEGMENT LTYPE(NBLIS,2) SEGMENT LTYPE2(NBLIS2,4) SEGMENT ICHPO REAL*8 XR(NBPTCH),YR(NBPTCH),ZR(NBPTCH),DEN(NBPTCH) ENDSEGMENT SEGMENT IDENSI REAL*8 DENSI(NBNDS) ENDSEGMENT 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 On lit le maillage a raffiner et le champoint de densite a affecter CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) IF (IERR.NE.0) RETURN CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU) IF (IERR.NE.0) RETURN C C======================================================================= C Initialisations - Declarations C======================================================================= NBOUCL=0 NPTINI=(XCOOR(/1)/(IDIM+1)) SEGACT IPT1 IJKLMN=4 SEGINI IWORK1 IWORK1(1)=4 IWORK1(2)=8 IWORK1(3)=6 IWORK1(4)=10 C C======================================================================= C Tests sur les donnees C======================================================================= C On teste pour savoir si le maillage IPT1 (maillage 1) a un type C d'elements implemante. NBLIS=MAX(1,IPT1.LISOUS(/1)) NSURE=0 SEGINI LTYPE DO 1 INB=1,MAX(1,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF NT=IPT2.ITYPEL IF (NT.EQ.48) NSURE=NSURE+1 LTYPE(INB,1)=INB LTYPE(INB,2)=NT 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 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 ((IPT1.LISOUS(/1).EQ.0).AND.(NT.EQ.48)) RETURN IF (IPT1.LISOUS(/1).NE.0) SEGDES IPT2 1 CONTINUE C 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 : IPT1 LCONF2=0 IF (NSURE.EQ.0) GOTO 7 LCONF2=1 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 C C------------------------------ C B) Creation du nouveau IPT1 ! C------------------------------ NBSOUS=LTYPE(/1)-NSURE NBREF=0 C C B.1/** Cas ou il n'y a qu'un lisous de type different de 48 IF (NBSOUS.EQ.1) 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 NBNN=0 NBELEM=0 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 7 CONTINUE 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 9 CONTINUE SEGACT IPT1 NBOUCL=NBOUCL+1 SEGINI ICPR NBNDS=0 NCOL=0 DO 11 INB=1,MAX(1,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF IF (IPT2.ITYPEL.EQ.48) GOTO 11 DO 10 J=1,IPT2.NUM(/2) DO 10 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 10 CONTINUE IF (IPT1.LISOUS(/1).NE.0) SEGDES IPT2 11 CONTINUE IF (IDIM.EQ.2) NCOL=NCOL*2 IF (IDIM.EQ.3) NCOL=NCOL*4 C C C C C======================================================================= C Prise en compte des noeuds de relations C======================================================================= SEGINI IRELA IF (LCONF2.EQ.0) GOTO 1200 SEGACT IPT6 DO 1218 LAA=1,IPT6.LISOUS(/1) IPT2=IPT6.LISOUS(LAA) 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 SEGDES IPT2 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 existees dans le tableau XCOOR. C C------------------------------------------ C A) Recuperation des valeurs du champoint ! C------------------------------------------ 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) SEGINI ICHPO DO 13 I=1,NBPTCH J=IPT3.NUM(1,I) XR(I)=XCOOR((J-1)*(IDIM+1)+1) YR(I)=XCOOR((J-1)*(IDIM+1)+2) IF (IDIM.EQ.3) ZR(I)=XCOOR((J-1)*(IDIM+1)+3) DEN(I)=MPOVA1.VPOCHA(I,1) 13 CONTINUE SEGINI IDENSI C C-------------------------------------------------------------- C B) Affectation d'une densite a chacun des noeuds du maillage ! C-------------------------------------------------------------- DO 15 INB=1,MAX(1,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF IF (IPT2.ITYPEL.EQ.48) GOTO 15 DO 14 J=1,IPT2.NUM(/2) DO 14 I=1,IPT2.NUM(/1) NGLOB=IPT2.NUM(I,J) NLOC=ICPR(NGLOB,1) XPT=XCOOR((NGLOB-1)*(IDIM+1)+1) YPT=XCOOR((NGLOB-1)*(IDIM+1)+2) ZPT=0 IF (IDIM.EQ.3) ZPT=XCOOR((NGLOB-1)*(IDIM+1)+3) CALL RAFDEN(ICHPO,XPT,YPT,ZPT,IDIM,DENPT) DENSI(NLOC)=DENPT 14 CONTINUE IF (IPT1.LISOUS(/1).NE.0) SEGDES IPT2 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,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF IF (IPT2.ITYPEL.EQ.48) GOTO 18 ILPL=LPL(IPT2.ITYPEL) ILPT=LPT(IPT2.ITYPEL) DO 17 J=1,IPT2.NUM(/2) DO 17 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 IF (NEXIST.EQ.0) THEN KARPOS(NMIN)=KARPOS(NMIN)+1 KARETE(NMIN,KARPOS(NMIN))=NMAX KARET2(NMIN,KARPOS(NMIN))=1 NBARET=NBARET+1 NCOL=MAX(NCOL,KARPOS(NMIN)) ENDIF IF (NEXIST.NE.0) THEN KARET2(NMIN,NEXIST)=KARET2(NMIN,NEXIST)+1 ENDIF 17 CONTINUE IF (IPT1.LISOUS(/1).NE.0) SEGDES IPT2 18 CONTINUE 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 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,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.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,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF IF (IPT2.ITYPEL.EQ.48) GOTO 51 DO 52 KELEM=1,IPT2.NUM(/2) JTYPEL=IPT2.ITYPEL JLTEL1=LTEL(1,JTYPEL) JLTEL2=LTEL(2,JTYPEL) 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=(XCOOR(/1)/(IDIM+1))+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 IHERED=0 DO 1022 MOI=1,IPT2.NUM(/1) IHERED=IHERED+IRELA(ICPR(IPT2.NUM(MOI,KELEM),1)) 1022 CONTINUE 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) IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2) 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) IF (NAA.EQ.NPTA.AND.NBB.EQ.NPTB.AND.NCC.EQ.NPTC) THEN NEXIST=ITEMP1 JNOEFA(NEXIST,5)=LTYPE(INB,2) END IF 56 CONTINUE 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) IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2) 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, LCONF2 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 (LCONF2.EQ.0) THEN GOTO 28 ENDIF 22 CONTINUE NBRELA=0 SEGACT IPT6 DO 218 LAA=1,IPT6.LISOUS(/1) IPT2=IPT6.LISOUS(LAA) 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) KMILIE(NMIN,NEXIST)=NPTC IF (NEXIST.NE.0) KARET2(NMIN,NEXIST)=2 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 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 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 ELSE KARET2(NMIN2,NEXIS2)=2 ENDIF ENDIF C 220 CONTINUE SEGDES IPT2 218 CONTINUE IF (NBRELA.NE.0) GOTO 22 28 CONTINUE C C C======================================================================= C Pre-remplissage du tableau JFARAF C======================================================================= IF (IDIM.NE.3) GOTO 711 JNBFA=JNOEFA(/1) LLLL=10 SEGINI JFARAF SEGINI JPLAN2,JPLAN3 IF (LCONF2.EQ.0) GOTO 711 MPOS25=0 JNBSOM=0 1287 CONTINUE NBRELA=0 DO 670 MIE=1,IPT6.LISOUS(/1) IPT2=IPT6.LISOUS(MIE) 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=(XCOOR(/1)/(IDIM+1))+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) 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 IF (NEXIS2.EQ.0) THEN WRITE(*,*) 'Raffinement impossible, le champ de densite presente + des variations trop brusques' GOTO 998 ENDIF JNUMFA=JPLANS(NPTA,NEXIS2) 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.JNUMFA) NEXIS2=MAC 675 CONTINUE C C---------------------------- C B) Renseignement de JFARAF ! C---------------------------- JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 LTEMPO=(JPLAN3(NPTA,NEXIS2)-1)*2+1 JFARAF(JNUMFA,LTEMPO)=IPT2.NUM(1,MAA) IF (MPOS25.NE.0) JFARAF(JNUMFA,LTEMPO+1)=IPT2.NUM(MPOS25,MAA) IF (JNOEFA(JNUMFA,5).EQ.0) JNOEFA(JNUMFA,5)=48 671 CONTINUE 670 CONTINUE C IF (NBRELA.NE.0) GOTO 1287 711 CONTINUE C 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=0 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' SEGDES IMODE1 C KARAF2=0 DO 33 INB=1,MAX(1,IPT1.LISOUS(/1)) IF (IPT1.LISOUS(/1).EQ.0) THEN IPT2=IPT1 ELSE IPT2=IPT1.LISOUS(INB) SEGACT IPT2 ENDIF SEGACT IMODE1*MOD IMODE1.IMAMOD=IPT2 IMODE1.NEFMOD=IPT2.ITYPEL SEGDES IMODE1 C C------------------------------------------------------- C B) Creation d'un chamelem constant egal a 1 : MCHEL1 ! C------------------------------------------------------- CALL ZEROP(MMODE1,'NOEUD ',MCHEL1) SEGACT MCHEL1*MOD MCHEL1.IMACHE(1)=IPT2 IF (IDIM.EQ.2) MCHEL1.IFOCHE=-1 IF (IDIM.EQ.3) MCHEL1.IFOCHE=2 MCHAM1=MCHEL1.ICHAML(1) SEGACT MCHAM1 MELVA1=MCHAM1.IELVAL(1) SEGACT MELVA1*MOD MELVA1.VELCHE(1,1)=1. SEGDES MELVA1 SEGDES MCHAM1 C C------------------------------------------------------------- C C) Calcul de la surface ou du volume des elements : MCHEL2 ! C------------------------------------------------------------- CALL INTGCA(MMODE1,MCHEL1,0,1,IRET,XRET,MCHEL2) SEGACT IPT2 C 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 DO 32 J=1,IPT2.NUM(/2) XMOY=0 DO 29 I=1,IPT2.NUM(/1) NGLOB=IPT2.NUM(I,J) NLOC=ICPR(NGLOB,1) DENPT=DENSI(NLOC) XMOY=XMOY+DENPT/(IPT2.NUM(/1)) 29 CONTINUE XINT=MELVA2.VELCHE(1,J) XDIM=IDIM*1.0 IF ((XMOY-(XINT**(1/XDIM))).LT.0) THEN MELVA2.VELCHE(1,J)=1 KARAF=KARAF+1 ILPT=LPT(IPT2.ITYPEL) ILPL=LPL(IPT2.ITYPEL) 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 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 JTYPEL=IPT2.ITYPEL JLTEL1=LTEL(1,JTYPEL) JLTEL2=LTEL(2,JTYPEL) 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=(XCOOR(/1)/(IDIM+1))+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 C---------------------------- C F) Raffinement du maillage ! C---------------------------- NML=2 KARAF2=KARAF2+KARAF IF (KARAF.EQ.0) GOTO 33 C C 1/ Raffinement d'un maillage 2D IF (IDIM.EQ.2) CALL RAFF2D(IPT2,ICPR,KARPOS,KARETE,KMILIE, + MELVA2,NACREE,KARAF,MELEME,IDCP,IPT8,KARET2) 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) C b) Maillage 3D sans pyramides IF (NML.EQ.0) CALL RAFF3D(IPT2,ICPR,KARPOS,KARETE,KMILIE, + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA, + IPT8,JFARAF,KARET2) 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 IF (KARAF2.EQ.0) GOTO 33 IF (IPT8.NE.0) THEN C C--------------------------------------- C B) Traitement des anciennes relations ! C--------------------------------------- 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 C C 2/ Creation de IPT10 IF (NB10.NE.0) THEN MCOL=0 NBNN=9 NBELEM=NB10 NBREF=0 NBSOUS=0 SEGINI IPT10 IPT10.ITYPEL=48 C C 3/ Recherche des anciennes relations a reporter DO 628 MM11=1,IPT8.NUM(/2) IF (IPT8.NUM(MLONG,MM11).NE.0) GOTO 628 INRELA=IPT8.NUM(1,MM11) NEXIS6=0 DO 625 KOU=1,IPT6.LISOUS(/1) IPT2=IPT6.LISOUS(KOU) 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 C--------------------------------------- C C) Traitement des nouvelles relations ! C--------------------------------------- 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 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 NBLIS2=MCONF+MAX(1,MELEME.LISOUS(/1)) SEGINI LTYPE2 C C 3/ Renseignement des lisous non 48 dans LTYPE2 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 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 5/ Creation et renseignement de IPT6 IF (MCONF.NE.0) THEN LCONF2=1 NBNN=0 NBELEM=0 NBREF=0 NBSOUS=LTYPE2(/1) SEGINI IPT6 DO 613 IRR=1,LTYPE2(/1) IPT2=LTYPE2(IRR,3) SEGACT IPT2 IPT6.LISOUS(IRR)=IPT2 613 CONTINUE ELSE LCONF2=0 IPT6=MELEME ENDIF C C======================================================================= C Test sur la progression du raffinement C======================================================================= C IF (IPT1.LISOUS(/1).NE.0) SEGDES IPT2 C SEGDES IPT1,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE SEGSUP IPT9 SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1 SEGSUP MMODE1,MCHEL1,MCHEL2 IF (IDIM.EQ.3) SEGSUP JPLANS,JPLAN2,JPLAN3,JPLCOM,JNOEFA,IWORK2, +JFARAF C IF (KARAF2.NE.0) THEN IPT1=MELEME GOTO 9 ENDIF C 33 CONTINUE C C ### SORTIE TEMPORAIRE DU PROGRAMME ### 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 IF (MCONF.NE.0) MELEME=IPT6 C GOTO 998 C C======================================================================= C Mises a jour diverses C======================================================================= C ** On affecte une densite a chaque noeud que l'on a cree DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1) XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1)) 37 CONTINUE C C======================================================================= C Fin du programme C======================================================================= 998 CONTINUE CALL ECROBJ('MAILLAGE',MELEME) RETURN END