C VOLUME SOURCE SP204843 25/02/19 21:15:04 12161 C MODIF : O.STAB / 25.03.97 / APPEL A VOLOS QUI C AUTORISE UN RACCORD DE 2 GRILLES NON IDENTIQUE C FABRICATION DE CUBES ET PRISMES PAR TRANSLATION ET ROTATION ET C ENTRE SURFACES OPPOSEES C MODIFICATION AOUT 1984 : MAILLAGE AUTOMATIQUE A L'INTERIEUR D'UNE C SURFACE ENVELOPPE C DECEMBRE 1984 VERIFICATION QUE LES TOPOLOGIES DU HAUT ET DU BAS C SONT SIMILAIRES (MEMES ICPR) C JANVIER 1985 NOMBRE DE COUCHES IMPOSE SI INBR NEGATIF C NOVEMBRE 1985 OPTION GENERATRICE (APPEL VOLUMG) c 12/97 KICH modif en liaison avec evolution de PROPER SUBROUTINE VOLUME IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC SMELEME -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC CCGEOME -INC CCREEL *- -INC SMLREEL -INC TDEMAIT logical ltelq SEGMENT TABPAR(NCOUCH) SEGMENT ICPR(NBNNEL,NBELEC) CHARACTER*4 MCLE(7) DATA MCLE /'TRAN','ROTA','DINI','DFIN','GENE','PROG','VERB'/ IDIMP1 = IDIM + 1 c MLREEL=0 IMPOI =0 IMPOF =0 ipt3 =0 ipt4 =0 IVERB =0 C Pour l'optimiseur (SIGSEV parfois en DEBUG) XDIS =0.D0 YDIS =0.D0 ZDIS =0.D0 DEN1I = 0.D0 DEN2I = 0.D0 c IF (ILCOUR.LT.14.OR.ILCOUR.GT.17) then IF (ILCOUR.LT.23.OR.ILCOUR.GT.26) CALL ERREUR(16) ENDIF IF (IERR.NE.0) RETURN INBR=0 ICLE=3 CALL LIRENT(INBR,0,IREINB) 80 CALL LIRMOT(MCLE,7,JCLE,0) IF (JCLE.EQ.0) GOTO 87 GOTO (81,82,83,84,79,78,75),JCLE C---- mot-clé 'TRAN' : 81 CONTINUE IF (ICLE.NE.3) GOTO 86 ICLE=JCLE GOTO 80 C---- mot-clé 'ROTA' : 82 CONTINUE IF (ICLE.NE.3) GOTO 86 C LECTURE N1 AU CAS OU IL SOIT DERRIERE LE MOT CLE ROTA C L'UTILISATEUR PEUT AUSSI DONNER L'ANGLE AVEC UN ENTIER IF (IREINB.EQ.0) THEN IREIN2=0 CALL LIRENT(INBR,0,IREIN2) IF (IREIN2.EQ.1) THEN CALL LIRREE(XXX,0,IRETOU) IF (IRETOU.EQ.0) THEN XXX=INBR ELSE IREINB = IREIN2 ENDIF ELSE CALL LIRREE(XXX,1,IRETOU) IF (IERR.NE.0) RETURN ENDIF ELSE CALL LIRREE(XXX,1,IRETOU) IF (IERR.NE.0) RETURN ENDIF ANGLI=XXX C write(6,*) 'ANGLI,INBR,DEN1I,DEN2I=',ANGLI,INBR,DEN1I,DEN2I CALL LIROBJ('POINT ',IP1,1,IRETOU) CALL LIROBJ('POINT ',IP2,1,IRETOU) IF (IERR.NE.0) RETURN ICLE=JCLE GOTO 80 C---- mot-clé 'DINI' : 83 IF(IMPOI.EQ.1) GOTO 86 IMPOI=1 CALL LIRREE(XXX,1,IRETOU) DEN1I=XXX IF (DEN1I.LE.0.) THEN CALL ERREUR(17) RETURN ENDIF IF (IERR.NE.0) RETURN GOTO 80 C---- mot-clé 'DFIN' : 84 IF (IMPOF.EQ.1) GOTO 86 IMPOF=1 CALL LIRREE(XXX,1,IRETOU) DEN2I=XXX IF (DEN2I.LE.0.) THEN CALL ERREUR(17) RETURN ENDIF IF (IERR.NE.0) RETURN GOTO 80 C---- mot-clé 'GENE' : 79 CONTINUE CALL VOLUMG RETURN C---- mot-clé 'PROG' : 78 CONTINUE CALL LIROBJ('LISTREEL',MLREEL,1,IRETOU) IF (IERR.NE.0) RETURN GOTO 80 C---- mot-clé 'VERB' : 75 CONTINUE IVERB=1 GOTO 80 86 CONTINUE CALL REFUS C ... Fin de lecture des mots clés ... 87 CONTINUE CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) C ... Si pas d'option TRAN ni ROTA on veut lire un deuxième maillage ... IF (ICLE.EQ.3) CALL LIROBJ('MAILLAGE',IPT2,0,IRETOU) C ... S'il n'y en a pas, on va remplir l'enveloppe .(ou épaisseur).. IF(IRETOU.EQ.0) GOTO 4400 C ================================== C ------- DEBUT DE MODIF - O.STAB 05.12.96 -------- C ================================== C C WRITE(IOIMP,*) ' ICLE = ',ICLE,DEN1I,DEN2I IF((ICLE.NE.1).AND.(ICLE.NE.2).AND.(ICLE.NE.5))THEN CALL LIROBJ('POINT ',IPO1,0,IRETOU) IF( IRETOU.EQ.0)GOTO 871 CALL LIROBJ('POINT ',IPO2,0,IRETOU) IF( IRETOU.EQ.0)GOTO 871 C C calcul des densités moyennes si pas données C DEN1D=1. IF(IMPOI.EQ.0.AND.INBR.EQ.0) THEN MELEME = IPT1 SEGACT MELEME NBNN=NUM(/1) NBELEM=NUM(/2) NPR=0 DEN1=0.D0 DO 710 I=1,NBNN DO 7101 J=1,NBELEM IR=NUM(I,J) IF (IR.EQ.0) GOTO 7101 IREF=(IR-1)*IDIMP1 NPR=NPR+1 DEN1=DEN1+XCOOR(IREF+4) 7101 CONTINUE 710 CONTINUE DEN1D=DEN1/NPR SEGDES MELEME ENDIF DEN2D=1. IF(IMPOF.EQ.0.AND.INBR.EQ.0) THEN MELEME = IPT2 SEGACT MELEME NBNN=NUM(/1) NBELEM=NUM(/2) NPR=0 DEN2=0.D0 DO 711 I=1,NBNN DO 7111 J=1,NBELEM IR=NUM(I,J) IF (IR.EQ.0) GOTO 7111 IREF=(IR-1)*IDIMP1 NPR=NPR+1 DEN2=DEN2+XCOOR(IREF+4) 7111 CONTINUE 711 CONTINUE DEN2D=DEN2/NPR SEGDES MELEME ENDIF CALL VOLOS(IPT1,IPT2,IPO1,IPO2,DEN1D,DEN2D,INBR) RETURN ENDIF C ================================== C ------- FIN DE MODIF - O.STAB 05.12.96 -------- C ================================== C 871 CONTINUE C ... Début du traitement ... ISVOL1=0 ISVOL2=0 SEGACT IPT1 C SI IPT1 VOLUME IL FAUT EN EXTRAIRE LA FACE 1 C ... dans la pratique le maillage initial soit n'a pas de C sous-maillages, soit il en a 2 (triangles et quadrilatères), C sinon la programmation ci-dessous poserait des problèmes : C en cas de IPT1.LISOUS(/1) > 2 un saut immédiat vers 3101 C provoquerait SEGDES IPT3,IPT4 qui ne sont pas encore initialisés ... 3100 IF (IPT1.LISOUS(/1).EQ.0) GOTO 1000 IF (IPT1.LISOUS(/1).NE.2) GOTO 3101 IDEUX=2 IPT3=IPT1.LISOUS(1) IPT4=IPT1.LISOUS(2) SEGACT IPT3,IPT4 IP=IPT3.ITYPEL*IPT4.ITYPEL c ... TRI3*QUA4 TRI6*QUA8 ... IF (IP.NE.32.AND.IP.NE.60) GOTO 3101 IS=IPT3.ITYPEL+IPT4.ITYPEL c ... TRI3+QUA4 TRI6+QUA8 ... IF (IS.NE.12.AND.IS.NE.16) GOTO 3101 c ... ici on a deux maillages : un composé de triangles, l'autre de quadrilatères ... INCR=1 IF (IS.EQ.16) INCR=2 c ... NBNNEL = nombre de noeuds / élément du maillage total ... NBNNEL=4*INCR C EN FAIT ON CREE UN SEGMENT QUI CONTIENT LES CUBES ET LES TRIANGLES C 0 DANS LA DERNIERE POSITION DU TRIANGLE NBSOUS=0 NBREF=0 NBNN=NBNNEL NBELE3=IPT3.NUM(/2) IF (IPT3.ITYPEL.LE.6) NBTRI=NBELE3 IF (IPT3.ITYPEL.GE.8) NBQUA=NBELE3 NBELE4=IPT4.NUM(/2) IF (IPT4.ITYPEL.LE.6) NBTRI=NBELE4 IF (IPT4.ITYPEL.GE.8) NBQUA=NBELE4 NBELEM=NBELE3+NBELE4 SEGINI MELEME C*C ... Initialisation du nouveau maillage à 0 (est-ce nécessaire ?) ... C* DO 1100 I=1,NBNN C* DO 1100 J=1,NBELEM C* NUM(I,J)=0 C* 1100 CONTINUE C ... On transvase IPT3 dans le nouveau maillage ... DO 1101 J=1,NBELE3 ICOLOR(J)=IPT3.ICOLOR(J) DO 11011 I=1,IPT3.NUM(/1) NUM(I,J)=IPT3.NUM(I,J) 11011 CONTINUE 1101 CONTINUE C ... On transvase IPT4 dans le nouveau maillage ... DO 1102 J=1,NBELE4 K=J+NBELE3 ICOLOR(K)=IPT4.ICOLOR(J) DO 11021 I=1,IPT4.NUM(/1) NUM(I,K)=IPT4.NUM(I,J) 11021 CONTINUE 1102 CONTINUE GOTO 1001 C RECHERCHE DE LA PREMIERE FACE DE IPT1 3101 if (ipt3.ne.0) segdes ipt3 if (ipt4.ne.0) segdes ipt4 3102 IF (IPT1.LISREF(/1).LT.2) CALL ERREUR(16) IF (IERR.NE.0) RETURN ISVOL1=IPT1 IAUX=IPT1.LISREF(2) SEGDES IPT1 IPT1=IAUX SEGACT IPT1 GOTO 3100 c ... On vient ici si le maillage est simple (homogène) ... 1000 CONTINUE IDEUX=1 NBNNEL=IPT1.NUM(/1) IF (IPT1.ITYPEL.NE.8.AND.IPT1.ITYPEL.NE.10.AND.IPT1.ITYPEL.NE.4 #.AND.IPT1.ITYPEL.NE.6) GOTO 3102 INCR=1 IF (KDEGRE(IPT1.ITYPEL).EQ.3) INCR=2 MELEME=IPT1 c ... ici MELEME est le maillage contanant tous les éléments de la surface initiale ... 1001 SEGACT MCOORD*mod c ... SI c'est 'ROTA' ... IF (ICLE.EQ.2) GOTO 1 c ... si ce n'est ni ROTA ni TRAN ... IF (ICLE.EQ.3) GOTO 2 c---- OPTION 'TRAN' C LECTURE DES ARGUMENTS CALL LIROBJ('POINT ',IVEC,1,IRETOU) IF (IERR.NE.0) RETURN C AU CAS OU LE DECOUPAGE (N1) EST DERRIERE LE MOT CLE TRAN IF (IREINB.EQ.0) CALL LIRENT(INBR,0,IREINB) C VECTEUR TRANSLATION IREF=(IVEC-1)*IDIMP1 XTRAN=XCOOR(IREF+1) YTRAN=XCOOR(IREF+2) ZTRAN=XCOOR(IREF+3) DEN2=XCOOR(IREF+4) DLONG=SQRT(XTRAN*XTRAN+YTRAN*YTRAN+ZTRAN*ZTRAN) XDIS=XTRAN YDIS=YTRAN ZDIS=ZTRAN GOTO 2 C---- OPTION 'ROTA' ... 1 CONTINUE ANGLE=ANGLI*XPI/180.D0 c ... et des deux points définissant l'axe de rotation ... IREF1=(IP1-1)*IDIMP1 IREF2=(IP2-1)*IDIMP1 XPT1=XCOOR(IREF1+1) YPT1=XCOOR(IREF1+2) ZPT1=XCOOR(IREF1+3) XVEC=XCOOR(IREF2+1)-XCOOR(IREF1+1) YVEC=XCOOR(IREF2+2)-XCOOR(IREF1+2) ZVEC=XCOOR(IREF2+3)-XCOOR(IREF1+3) DEN2=(XCOOR(IREF2+4)+XCOOR(IREF1+4))*0.5 RAY=SQRT(XVEC*XVEC+YVEC*YVEC+ZVEC*ZVEC) XVEC=XVEC/RAY YVEC=YVEC/RAY ZVEC=ZVEC/RAY IF (ANGLE.GE.0.) GOTO 2 ANGLE=-ANGLE XVEC=-XVEC YVEC=-YVEC ZVEC=-ZVEC c ... calcul des moyennes des densités et des coordonnées ... 2 CONTINUE NBNN=NUM(/1) NBELEM=NUM(/2) NPR=0 DEN1=0. XG=0. YG=0. ZG=0. DO 5 I=1,NBNN DO 51 J=1,NBELEM IR=NUM(I,J) IF (IR.EQ.0) GOTO 51 IREF=(IR-1)*IDIMP1 NPR=NPR+1 DEN1=DEN1+XCOOR(IREF+4) XG=XG+XCOOR(IREF+1) YG=YG+XCOOR(IREF+2) ZG=ZG+XCOOR(IREF+3) 51 CONTINUE 5 CONTINUE DEN1=DEN1/NPR XG=XG/NPR YG=YG/NPR ZG=ZG/NPR c ... cas 'TRAN' => GOTO 6 ... IF (ICLE.EQ.1) GOTO 6 c ... cas 'ROTA' => GOTO 3 ... IF (ICLE.EQ.2) GOTO 3 c ... cas du volume entre deux surfaces ... C COMPATIBILITE DU 2EME OBJET ET RECHERCHE DU CENTRE DE GRAVITE SEGACT IPT2 3150 IF (IPT2.LISOUS(/1).EQ.0) GOTO 1020 IF (IDEUX.NE.2) CALL ERREUR(21) IF (IERR.NE.0) RETURN IF (IPT2.LISOUS(/1).NE.2) GOTO 3151 IPT5=IPT2.LISOUS(1) IPT6=IPT2.LISOUS(2) SEGACT IPT5,IPT6 IF (IPT5.ITYPEL.NE.IPT3.ITYPEL) GOTO 3151 IF (IPT6.ITYPEL.NE.IPT4.ITYPEL) GOTO 3151 IF (IPT5.NUM(/2).NE.IPT3.NUM(/2)) CALL ERREUR(21) IF (IPT6.NUM(/2).NE.IPT4.NUM(/2)) CALL ERREUR(21) IF (IPT5.NUM(/1).NE.IPT3.NUM(/1)) CALL ERREUR(21) IF (IPT6.NUM(/1).NE.IPT4.NUM(/1)) CALL ERREUR(21) IF (IERR.NE.0) RETURN GOTO 1021 3151 SEGDES IPT5,IPT6 3152 IF (IPT2.LISREF(/1).LT.2) CALL ERREUR(21) IF (IERR.NE.0) RETURN IAUX=IPT2.LISREF(1) ISVOL2=IPT2 SEGDES IPT2 IPT2=IAUX SEGACT IPT2 GOTO 3150 c ... les deux maillages doivent être simples ... 1020 IF (IPT1.ITYPEL.NE.IPT2.ITYPEL) GOTO 3152 IF (IDEUX.NE.1) GOTO 3152 IF (IPT2.NUM(/1).NE.IPT1.NUM(/1)) CALL ERREUR(21) IF (IPT2.NUM(/2).NE.IPT1.NUM(/2)) CALL ERREUR(21) 1021 CONTINUE c ... calcul des densités et coordonnées moyennes de la deuxième surface ... NPR=0 DEN2=0. DLONG=0. IPT7=IPT2 IPT9=IPT1 M1031=2 IF (IDEUX.EQ.1) M1031=1 DO 1031 M=1,M1031 IF (M1031.NE.1) IPT7=IPT2.LISOUS(M) IF (M1031.NE.1) IPT9=IPT1.LISOUS(M) DO 4 I=1,IPT7.NUM(/1) DO 41 J=1,IPT7.NUM(/2) IREF=(IPT7.NUM(I,J)-1)*IDIMP1 IRE1=(IPT9.NUM(I,J)-1)*IDIMP1 DEN2=DEN2+XCOOR(IREF+4) XVAL=(XCOOR(IREF+1)-XCOOR(IRE1+1))**2 XVAL=(XCOOR(IREF+2)-XCOOR(IRE1+2))**2+XVAL XVAL=(XCOOR(IREF+3)-XCOOR(IRE1+3))**2+XVAL DLONG=DLONG+SQRT(XVAL) 41 CONTINUE 4 CONTINUE NPR=NPR+IPT7.NUM(/1)*IPT7.NUM(/2) 1031 CONTINUE DEN2=DEN2/NPR DLONG=DLONG/FLOAT(NPR) GOTO 6 c ... cas 'ROTA' ... 3 CONTINUE XV1=XG-XPT1 YV1=YG-YPT1 ZV1=ZG-ZPT1 PV1=XV1*XVEC+YV1*YVEC+ZV1*ZVEC XV1=XV1-PV1*XVEC YV1=YV1-PV1*YVEC ZV1=ZV1-PV1*ZVEC RAY=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) XV1=XV1/RAY YV1=YV1/RAY ZV1=ZV1/RAY XV2=YVEC*ZV1-ZVEC*YV1 YV2=ZVEC*XV1-XVEC*ZV1 ZV2=XVEC*YV1-YVEC*XV1 C RAYON MOYEN C ANGLE EN RADIANS D'OU LONGUEUR MOYENNE DLONG=RAY*ABS(ANGLE) c ... partie commune, recherche du nombre de couches à créér et des densités ... 6 CONTINUE IF (IMPOI.EQ.1) DEN1=DEN1I IF (IMPOF.EQ.1) DEN2=DEN2I C JE NE VOIS PAS DANS QUELS CAS CA INTERVIENT CALL LIRENT(INBR,0,IRETOU) C write(6,*) 'volume:DLONG,DEN1I,DEN2I =',DLONG,DEN1I,DEN2I IPT3=MELEME DENI=DEN1 DECA=DEN2-DEN1 DENM = ABS(MAX(DEN1,DEN2)) if (abs(dlong).lt.10.D0*XZPREC*DENM) dlong=1.d0 DEN1=DEN1/DLONG DEN2=DEN2/DLONG IF (MLREEL.NE.0)THEN SEGACT,MLREEL INBR=PROG(/1)-1 SEGDES,MLREEL ENDIF C write(6,*) 'volume:DLONG,DEN1, DEN2, INBR =',IMPOI,DEN1,DEN2,INBR CALL DECOUP(INBR,DEN1,DEN2,APROG,NCOUCH,DENI,DECA,DLONG) IF (IERR.NE.0) RETURN C write(6,*) 'NCOUCH =',NCOUCH NX=NCOUCH-1 IF (IIMPI.EQ.1) WRITE(IOIMP,9000) NCOUCH,APROG 9000 FORMAT(/' COUCHES ',I6,' RAISON ',G12.5) C ... Initialisation du nouveau maillage ... C ON FAIT TOUJOURS COMME SI IL N'Y AVAIT QU'UN TYPE D'ELEMENT NBSOUS=0 C MODIF POUR CONSTRUIRE TOUJOURS LE POURTOUR NBREF=3 IF (IPT1.LISREF(/1).NE.0) NBREF=3 NBNN=2*NBNNEL+(INCR-1)*(NBNNEL/2) NBNNV=NBNN NBASE=NBELEM NBELEM=NBELEM*NCOUCH SEGINI IPT7 IF (NBNNV.EQ.6 ) IPT7.ITYPEL=16 IF (NBNNV.EQ.15) IPT7.ITYPEL=17 IF (NBNNV.EQ.8 ) IPT7.ITYPEL=14 IF (NBNNV.EQ.20) IPT7.ITYPEL=15 IPT7.LISREF(1)=IPT1 C*c ... Mise à 0 des connectivités ... C* DO 1040 I=1,NBNN C* DO 1040 J=1,NBELEM C* IPT7.NUM(I,J)=0 C* 1040 CONTINUE SEGINI TABPAR c ... si ce n'est ni TRAN ni ROTA on saute ... IF (ICLE.EQ.3) GOTO 16 IF (ICLE.EQ.2) GOTO 10 c ... cas TRAN, on fait appel à l'opérateur PLUS ... IOPTG=1 CALL ECROBJ('POINT ',IVEC) CALL ECROBJ('MAILLAGE',IPT1) CALL PROPER(IOPTG) GOTO 11 c ... cas ROTA, on fait appel à l'opérateur TOURner ... 10 XXX=ANGLI CALL ECRREE(XXX) CALL ECROBJ('POINT ',IP2) CALL ECROBJ('POINT ',IP1) CALL ECROBJ('MAILLAGE',IPT1) CALL TOURNE 11 CONTINUE c ... puis on lit le second MAILLAGE ... CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) IF (IERR.NE.0) RETURN C IPT3 ET IPT4 ONT ETE DESCENDU DANS L'OPERATION AINSI QUE MCOORD/REFPO 16 SEGACT IPT1,IPT2,MCOORD IPT4=IPT2 c ... si le 1er maillage est simple on suppose que le deuxième le sera aussi ... IF (IDEUX.EQ.1) GOTO 15 IPT5=IPT2.LISOUS(1) IPT6=IPT2.LISOUS(2) SEGACT IPT5,IPT6 C ON FAIT COMME POUR LE BAS NBSOUS=0 NBREF=0 NBNN=4*INCR c ... qui ont les mêmes nombres d'éléments que les sous-maillages du premier ... NBELEM=NBELE3+NBELE4 SEGINI MELEME c ... et que l'on transvase succesivement dans une nouvelle entité (MELEME) ... C* DO 1110 J=1,NBELEM C* DO 1110 I=1,NBNN C* NUM(I,J)=0 C* 1110 CONTINUE c ... d'abord IPT5 ... DO 1111 J=1,NBELE3 ICOLOR(J)=IPT5.ICOLOR(J) DO 11111 I=1,IPT5.NUM(/1) NUM(I,J)=IPT5.NUM(I,J) 11111 CONTINUE 1111 CONTINUE c ... puis IPT6 ... DO 1112 J=1,NBELE4 K=J+NBELE3 ICOLOR(K)=IPT6.ICOLOR(J) DO 11121 I=1,IPT6.NUM(/1) NUM(I,K)=IPT6.NUM(I,J) 11121 CONTINUE 1112 CONTINUE SEGDES IPT5,IPT6,IPT2 c ... et qui remplace le maillage lu ... IPT4=MELEME 15 IPT7.LISREF(2)=IPT2 C CONSTRUCTION DE LA TABLE DES POINTS EFFECTIFS c ... IPT3 = maillage (parfois bâtard) contenant toutes les facettes c de la surface initiale (?) ... NBELEC=IPT3.NUM(/2) c ... ICPR(ligne = nombre maxi de noeuds / facette, colonne = nb facettes) ... SEGINI ICPR C*c ... mise à 0 ... C* DO 12 I=1,NBNNEL C* DO 12 J=1,NBELEC C* 12 ICPR(I,J)=0 c ... on parcourt les 2 maillages ... DO 13 J=1,NBELEC DO 131 I=1,NBNNEL c ... IR = N° du noeud (ou 0) du 1er ... IR=IPT3.NUM(I,J) c ... IR2 = N° du noeud (ou 0) equivalent du 2nd ... IR2=IPT4.NUM(I,J) c ... si le 1er est absent ... IF (IR.EQ.0) GOTO 1120 c ... sinon, si son equivalent est nul => kk !!!!!!! ... IF (IR2.EQ.0) GOTO 8833 I1=IR I1R2=IR2 c ... si ce n'est pas le 1er élément ... IF (J.EQ.1) GOTO 131 c ... on va vérifier que l'equivalence est la même pour c tous les éléments précédents ... JM1=J-1 DO 14 JJ=1,JM1 DO 141 II=1,NBNNEL IR=IPT3.NUM(II,JJ) IR2=IPT4.NUM(II,JJ) IF (IR.EQ.0) GOTO 141 IF (IR.NE.I1) GOTO 8834 IF (IR2.NE.I1R2) GOTO 8833 c ... on met dans ICPR(n° noeud,n° élt) la valeur de II+(JJ-1)*8 c tq. le noeud II de l'élt JJ de IPT3 est le même que c le noeud I de l'élt J (toutefois JJ < J donc aucun noeud ne c pointe sur lui même) ... ICPR(I,J)=II+(JJ-1)*8 GOTO 131 8834 IF (IR2.EQ.I1R2) GOTO 8833 141 CONTINUE 14 CONTINUE GOTO 131 c ... si le 1er est absent (suite), on met ICPR correspondant à -1 ... 1120 ICPR(I,J)=-1 c ... si son equivalent est non nul => kk !!!!!!! ... IF (IR2.NE.0) GOTO 8833 131 CONTINUE 13 CONTINUE GOTO 8835 8833 CONTINUE C LES TOPOLOGIES SONT DIFFERENTES SEGSUP ICPR CALL ERREUR(21) RETURN 8835 CONTINUE C ON FABRIQUE POUR LE MOMENT DES CUBES A 8 OU 20 NOEUDS ET DES PRISMES C A 6 OU 15 NOEUDS C D'ABORD LES POINTS DU BAS DIN=DEN1 DO 20 I=1,NBELEC c ... On commence par donner la couleur ... c ... Si c'est TRAN ou ROTA ce sera celle du maillage d'origine ... IF (ICLE.NE.3) THEN IPT7.ICOLOR(I)=IPT3.ICOLOR(I) c ... sinon, une <> au sens de ITABM ... ELSE ICOLI=IPT3.ICOLOR(I) C ... CORRECTION PROBLEME SAUSSAIS 29 NOVEMBRE 1985 ICOLJ=IPT4.ICOLOR(I) IPT7.ICOLOR(I)=ITABM(ICOLI,ICOLJ) ENDIF c ... Puis on commence par transvaser les connectivités de la surface initiale ... DO 201 J=1,NBNNEL IR=IPT3.NUM(J,I) IF (IR.EQ.0) GOTO 201 IPT7.NUM(J,I)=IR 201 CONTINUE 20 CONTINUE IBASE=nbpts C ON FABRIQUE ENSUITE LES COUCHES C ON AFFECTE SEULEMENT LES NUMEROS DE NOEUDS c ... IDIF = nombre de noeuds placés <> les couches ... IDIF=(INCR-1)*(NBNNEL/2) NX=NCOUCH-1 DO 21 ICOUCH=1,NCOUCH DIN=DIN*APROG TABPAR(ICOUCH)=DIN IF (ICOUCH.EQ.NCOUCH) GOTO 21 JBASE=(ICOUCH-1)*NBELEC IF (INCR.EQ.1) GOTO 2000 C ON FABRIQUE D'ABORD LA COUCHE INTERMEDIAIRE DO 2001 J=1,NBELEC IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) DO 20011 IA=1,(NBNNEL/2) I=2*IA-1 IF (ICPR(I,J).EQ.-1) GOTO 20011 IF (ICPR(I,J).NE.0) GOTO 2002 IBASE=IBASE+1 IPT7.NUM(IA+NBNNEL,J+JBASE)=IBASE GOTO 20011 2002 IAUX=ICPR(I,J) JJ=(IAUX-1)/8+1 II=IAUX-8*JJ+8 IIA=(II+1)/2 IPT7.NUM(IA+NBNNEL,J+JBASE)=IPT7.NUM(IIA+NBNNEL,JJ+JBASE) 20011 CONTINUE 2001 CONTINUE 2000 CONTINUE DO 22 J=1,NBELEC IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) DO 221 I=1,NBNNEL IF (ICPR(I,J).EQ.-1) GOTO 221 IF (ICPR(I,J).NE.0) GOTO 23 IBASE=IBASE+1 IPT7.NUM(I,J+JBASE+NBELEC)=IBASE IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IBASE GOTO 221 23 IAUX=ICPR(I,J) JJ=(IAUX-1)/8+1 II=IAUX-8*JJ+8 IPT7.NUM(I,J+JBASE+NBELEC)=IPT7.NUM(II,JJ+JBASE+NBELEC) IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IPT7.NUM(II+NBNNEL+IDIF, & JJ+JBASE) 221 CONTINUE 22 CONTINUE 21 CONTINUE IF (MLREEL.NE.0)THEN SEGACT,MLREEL DPROG=PROG(NCOUCH+1)-PROG(1) DO 12345 ICOUCH=1,NCOUCH TABPAR(ICOUCH)=(PROG(ICOUCH+1)-PROG(ICOUCH))/DPROG 12345 CONTINUE SEGDES,MLREEL ENDIF 25 CONTINUE C ON FAIT LES POINTS DU HAUT ET EVENTUELLEMENT LA COUCHE INTERMEDIAIRE C PRECEDENTE JBASE=NBELEC*NX IF (INCR.EQ.1) GOTO 2003 DO 2004 J=1,NBELEC IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) DO 20041 IA=1,(NBNNEL/2) I=2*IA-1 IF (ICPR(I,J).EQ.-1) GOTO 20041 IF (ICPR(I,J).NE.0) GOTO 2005 IBASE=IBASE+1 IPT7.NUM(IA+NBNNEL,J+JBASE)=IBASE GOTO 20041 2005 IAUX=ICPR(I,J) JJ=(IAUX-1)/8+1 II=IAUX-8*JJ+8 IIA=(II+1)/2 IPT7.NUM(IA+NBNNEL,J+JBASE)=IPT7.NUM(IIA+NBNNEL,JJ+JBASE) 20041 CONTINUE 2004 CONTINUE 2003 CONTINUE DO 30 J=1,NBELEC IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) DO 301 I=1,NBNNEL IF (ICPR(I,J).EQ.-1) GOTO 301 IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IPT4.NUM(I,J) 301 CONTINUE 30 CONTINUE DPAR=0. C CREATION DES POINTS IADR=nbpts NBPTS=IADR+NCOUCH*INCR*NBELEC*NBNNEL SEGADJ MCOORD DO 61 ICOUCH=1,NCOUCH DIN=TABPAR(ICOUCH) DO 610 IC=1,INCR IC1=INCR+1-IC DPAR=DPAR+DIN/INCR UMDPAR=1.-DPAR DEN610=DENI+DECA*DPAR IF (ICOUCH.EQ.NCOUCH.AND.IC.EQ.INCR) GOTO 610 IF (ICLE.NE.2) GOTO 63 ANG=DPAR*DLONG/RAY SI=SIN(ANG) CO=COS(ANG) 63 CONTINUE DO 620 J=1,NBELEC DO 62 I=1,NBNNEL,IC1 IF (ICPR(I,J).NE.0) GOTO 62 IREF=4*IPT3.NUM(I,J)-4 C write(6,*) 'XCOOR(/1)=',XCOOR(/1) C write(6,*) 'IADR,IREF,DPAR,XDIS=',IADR,IREF,DPAR,XDIS GOTO (67,64,66),ICLE 67 XCOOR(IADR*IDIMP1+1)=XCOOR(IREF+1)+DPAR*XDIS XCOOR(IADR*IDIMP1+2)=XCOOR(IREF+2)+DPAR*YDIS XCOOR(IADR*IDIMP1+3)=XCOOR(IREF+3)+DPAR*ZDIS GOTO 65 66 IREF2=4*IPT4.NUM(I,J)-4 XCOOR(IADR*IDIMP1+1)=UMDPAR*XCOOR(IREF+1)+DPAR*XCOOR(IREF2+1) XCOOR(IADR*IDIMP1+2)=UMDPAR*XCOOR(IREF+2)+DPAR*XCOOR(IREF2+2) XCOOR(IADR*IDIMP1+3)=UMDPAR*XCOOR(IREF+3)+DPAR*XCOOR(IREF2+3) GOTO 65 64 X1=XCOOR(IREF+1)-XPT1 Y1=XCOOR(IREF+2)-YPT1 Z1=XCOOR(IREF+3)-ZPT1 XV=X1*XV1+Y1*YV1+Z1*ZV1 YV=X1*XV2+Y1*YV2+Z1*ZV2 ZV=X1*XVEC+Y1*YVEC+Z1*ZVEC XD=XV*CO-YV*SI YD=XV*SI+YV*CO ZD=ZV XCOOR(IADR*IDIMP1+1)=XD*XV1+YD*XV2+ZD*XVEC+XPT1 XCOOR(IADR*IDIMP1+2)=XD*YV1+YD*YV2+ZD*YVEC+YPT1 XCOOR(IADR*IDIMP1+3)=XD*ZV1+YD*ZV2+ZD*ZVEC+ZPT1 GOTO 65 65 CONTINUE XCOOR((IADR+1)*IDIMP1)=DEN610 IADR=IADR+1 62 CONTINUE 620 CONTINUE 610 CONTINUE 61 CONTINUE NBPTS=IADR SEGADJ MCOORD 60 CONTINUE C C'EST FINI C IL RESTE DANS LE CAS OU ON A DES CUBES ET DES PRISMES A LES SEPARER C ET A SUPPRIMER LES SEGMENTS SUPPLEMENTAIRES DE TRAVAIL C D'ABORD FAIRE LE POURTOUR A PARTIR DU CONTOUR IF (IPT7.LISREF(/1).EQ.2) GOTO 3000 CALL ECROBJ('MAILLAGE',IPT1) CALL ECRCHA('NOID') CALL PRCONT CALL LIROBJ('MAILLAGE',IPT5,1,IRETOU) IF (IERR.NE.0) GOTO 3000 C IPT5 LE CONTOUR IPT6 SERA LE POURTOUR SEGACT IPT5 NBASE=IPT5.NUM(/2) NBNN=INCR*4 NBELEM=NBASE*NCOUCH NBSOUS=0 NBREF=0 SEGINI IPT6 IPT6.ITYPEL=6+2*INCR SEGACT IPT3 DO 3001 IEL=1,NBASE DO 30011 IP=1,INCR+1 INP=IPT5.NUM(IP,IEL) DO 3003 IELS=1,NBELEC DO 30031 IPS=1,NBNNEL IPSP=IPT3.NUM(IPS,IELS) IF (IPSP.EQ.0) GOTO 30031 IF (IPSP.EQ.INP) GOTO 3002 30031 CONTINUE 3003 CONTINUE GOTO 3000 3002 CONTINUE DO 3004 IC=1,NCOUCH IBASE=(IC-1)*NBASE JBASE=(IC-1)*NBELEC C PTS DU BAS IPT6.NUM(IP,IEL+IBASE)=IPT7.NUM(IPS,IELS+JBASE) C PTS DU HAUT IPT6.NUM(NBNN+2-INCR-IP,IEL+IBASE)= # IPT7.NUM(IPS+NBNNEL+IDIF,IELS+JBASE) C EVENTUELLEMENT PTS MILIEUX IF (INCR.EQ.1.OR.IP.EQ.2) GOTO 3004 IPT6.NUM(10-2*IP,IEL+IBASE)=IPT7.NUM((IPS+1)/2+NBNNEL,IELS+JBASE) 3004 CONTINUE 30011 CONTINUE 3001 CONTINUE DO 3005 I=1,NCOUCH DO 30051 J=1,NBASE IPT6.ICOLOR(J+(I-1)*NBASE)=IPT5.ICOLOR(J) 30051 CONTINUE 3005 CONTINUE SEGDES IPT5,IPT6 IPT7.LISREF(3)=IPT6 3000 CONTINUE * cas ou on a saute la creation de ipt7.lisref(3) avec le goto 3000 if (ipt7.lisref(ipt7.lisref(/1)).eq.0) then nbnn=ipt7.num(/1) nbelem=ipt7.num(/2) nbsous=ipt7.lisous(/1) nbref=ipt7.lisref(/1)-1 segadj ipt7 endif IF (IDEUX.EQ.1) GOTO 1500 SEGSUP IPT3,IPT4 MELEME=IPT7 NBSOUS=2 NBREF=LISREF(/1) NBNN=0 NBELEM=0 SEGINI IPT7 IPT7.LISREF(1)=LISREF(1) IPT7.LISREF(2)=LISREF(2) IF (NBREF.EQ.3) IPT7.LISREF(3)=LISREF(3) NBSOUS=0 NBREF=0 NBNN=6 IF (INCR.EQ.2) NBNN=15 NBELEM=NBTRI*NCOUCH SEGINI IPT3 IPT3.ITYPEL=16 IF (INCR.EQ.2) IPT3.ITYPEL=17 IPT7.LISOUS(1)=IPT3 NBNN=8 IF (INCR.EQ.2) NBNN=20 NBELEM=NBQUA*NCOUCH SEGINI IPT4 IPT4.ITYPEL=14 IF (INCR.EQ.2) IPT4.ITYPEL=15 IPT7.LISOUS(2)=IPT4 IT=0 IQ=0 DO 1501 J=1,NUM(/2) IF (NUM(NBNNV,J).EQ.0) GOTO 1502 C C'EST UN CUBE IQ=IQ+1 IPT4.ICOLOR(IQ)=ICOLOR(J) DO 1503 K=1,IPT4.NUM(/1) IPT4.NUM(K,IQ)=NUM(K,J) 1503 CONTINUE GOTO 1501 1502 IT=IT+1 IPT3.ICOLOR(IT)=ICOLOR(J) C C'EST UN PRISME IF (INCR.EQ.2) GOTO 2020 IPT3.NUM(1,IT)=NUM(1,J) IPT3.NUM(2,IT)=NUM(2,J) IPT3.NUM(3,IT)=NUM(3,J) IPT3.NUM(4,IT)=NUM(NBNNEL+1,J) IPT3.NUM(5,IT)=NUM(NBNNEL+2,J) IPT3.NUM(6,IT)=NUM(NBNNEL+3,J) GOTO 1501 2020 CONTINUE DO 2021 L=1,6 IPT3.NUM(L,IT)=NUM(L,J) 2021 CONTINUE IPT3.NUM(7,IT)=NUM(NBNNEL+1,J) IPT3.NUM(8,IT)=NUM(NBNNEL+2,J) IPT3.NUM(9,IT)=NUM(NBNNEL+3,J) DO 2022 L=1,6 IPT3.NUM(L+9,IT)=NUM(NBNNEL+IDIF+L,J) 2022 CONTINUE 1501 CONTINUE SEGDES IPT3,IPT4 SEGSUP MELEME 1500 SEGDES IPT1,IPT2 SEGSUP ICPR,TABPAR IF (ISVOL1.EQ.0) GOTO 3200 IPT8=ISVOL1 SEGACT IPT8 ltelq=.false. CALL FUSE(IPT8,IPT7,IRET,ltelq) SEGDES IPT7,IPT8 IPT7=IRET 3200 CONTINUE IF (ISVOL2.EQ.0) GOTO 3201 IPT8=ISVOL2 SEGACT IPT8 ltelq=.false. CALL FUSE(IPT7,IPT8,IRET,ltelq) SEGDES IPT7,IPT8 IPT7=IRET 3201 CONTINUE SEGDES IPT7 CALL ECROBJ('MAILLAGE',IPT7) RETURN 4400 CONTINUE mchpoi=0 epai=0.d0 Call LIRREE(EPAI,0,iretou) if(iretou.eq.0) call lirobj('CHPOINT ' , MCHPOI,0,iretch) if(iretou+iretch.eq.0) then C MAILLAGE AUTOMATIQUE DE VOLUME IF (IVERB.EQ.1) write(IOIMP,*) ' appel a demete' CALL DEMETE(IPT1) IF (IERR.NE.0) RETURN IPT7=IPT1 GOTO 3201 else call lirobj('POINT ',ip1,1,iretpt) if(ierr.ne.0) return call volshb(ipt1,epai,mchpoi,ip1,ipt7) if(ierr.ne.0) return go to 3201 endif END