C PROT SOURCE CB215821 24/04/12 21:16:58 11897 SUBROUTINE PROT(IPMODE,IPCHT,IPCHE,ITPR) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) *--------------------------------------------------------------------* * * * Sous-programme associé à l'opérateur CALP * * ____________________________________________ * * * * Projection d'un chamelem de temperature sur une géometrie * constituée de coques * * * * * * Auteur, date de création: * * ------------------------- * * * * Bruno VIGAN, le 26 février 1997. * * * *--------------------------------------------------------------------* * -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMMODEL -INC SMCHAML -INC SMELEME -INC SMCHPOI -INC SMLCHPO -INC SMLMOTS -INC TMTRAV SEGMENT VECT REAL*8 VEC1(IDIM) REAL*8 VEC2(IDIM) REAL*8 VECN(IDIM) ENDSEGMENT SEGINI VECT * SEGMENT ICPR(NBNOE,NCHAM) SEGMENT NKON(IKOUR) SEGMENT NUIN(IKOUR) * SEGMENT ICARAC REAL*8 XEPAI(NCHAM) REAL*8 XEXCE(NCHAM) ENDSEGMENT SEGMENT NCARAC(NCHAM) * MMODEL=IPMODE SEGACT,MMODEL NMOD=KMODEL(/1) * MCHELM=IPCHE SEGACT,MCHELM NCHAM=ICHAML(/1) * segact mcoord*mod NBNOE = nbpts SEGINI ICPR SEGINI ICARAC SEGINI NCARAC * DO 10, I = 1, NCHAM ICARAC.XEPAI(I)= 0. ICARAC.XEXCE(I)= 0. DO 10, J=1, NBNOE ICPR(J,I)=0 CONTINUE 10 CONTINUE NBCAR = 0 * * Création du maillage principal * NBSOUS = 0 NBREF = 0 NBELEM = 0 NBNN = 0 SEGINI IPT2 IKOUR=0 c listmots des phases ilphmo = -1 jgn = 8 jgm = nmod segini mlmots ilphmo = mlmots jgm = 1 * * Boucle sur l'ensemble des sous zones du modeles * DO 30, NUMO = 1, NMOD * IMODEL = KMODEL(NUMO) SEGACT,IMODEL * * Test si le modele est une coque * NUFOR = NUMMFR(NEFMOD) IF (NUFOR.EQ.3 .OR. NUFOR.EQ.5 .OR. NUFOR.EQ.9)THEN * * Recherche du chamemlem de caracteristique assossiée * NUCHA = 0 DO 15, NUCH = 1, NCHAM * IF ( CONCHE(NUCH).EQ.CONMOD.AND. C IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH * 15 CONTINUE * IF (NUCHA.NE.0) THEN MCHAML=ICHAML(NUCHA) SEGACT,MCHAML * XEXCE1 = 0. XEPAI1 = 0. NCOMP = IELVAL(/1) DO 20, I = 1, NCOMP IF (NOMCHE(I).EQ.'EPAI') THEN MELVAL = IELVAL(I) SEGACT, MELVAL XEPAI1 = VELCHE(1,1) ELSEIF (NOMCHE(I).EQ.'EXCE') THEN MELVAL = IELVAL(I) SEGACT, MELVAL XEXCE1 = VELCHE(1,1) ENDIF 20 CONTINUE * * recherche du numero de caracteristique associe * a l'epaisseur et l'excentricitee * NUCAR = 0 DO 22, I = 1, NBCAR IF (ICARAC.XEPAI(I).EQ.XEPAI1.AND. C ICARAC.XEXCE(I).EQ.XEXCE1) NUCAR = I 22 CONTINUE * IF (NUCAR.EQ.0) THEN NUCAR = NBCAR+1 ICARAC.XEPAI(NUCAR)=XEPAI1 ICARAC.XEXCE(NUCAR)=XEXCE1 NBCAR = NUCAR ENDIF NCARAC(NUCHA)=NUCAR * MELEME = IMAMOD SEGACT MELEME * * recherche du nombre de noeuds * DO 25 I=1, NUM(/1) DO 25 J=1, NUM(/2) ITH= NUM(I,J) IF (ICPR(ITH,NUCAR).EQ.0) THEN IKOUR=IKOUR+1 ICPR(ITH,NUCAR)=IKOUR ENDIF 25 CONTINUE ENDIF ENDIF * if (numo.eq.1) then mots(1) = conmod(17:24) else do ipl = 1,jgm if (mots(ipl).eq.conmod(17:24)) goto 27 enddo jgm = jgm + 1 mots(jgm) = conmod(17:24) 27 continue endif C 30 CONTINUE * segadj mlmots * Augmentation du tableau de coordonnées * NBPTS = NBNOE+3*IKOUR SEGADJ MCOORD * NNNO = IKOUR SEGINI NKON SEGINI NUIN * DO 40, I = 1, NNNO NKON(I)=0 DO 40, K = 1, IDIM XCOOR((NBNOE+I-1)*(IDIM+1)+K) = 0. XCOOR((NBNOE+I-1+IKOUR)*(IDIM+1)+K) = 0. XCOOR((NBNOE+I-1+2*IKOUR)*(IDIM+1)+K) = 0. 40 CONTINUE * * Boucle sur l'ensemble des sous zones du modeles * DO 100, NUMO = 1, NMOD IMODEL = KMODEL(NUMO) * * Test si le modele est une coque * NUFOR = NUMMFR(NEFMOD) IF (NUFOR.EQ.3 .OR. NUFOR.EQ.5 .OR. NUFOR.EQ.9) THEN * * Recherche du chamemlem de caracteristique assossiée * NUCHA = 0 DO 50, NUCH = 1, NCHAM * IF ( CONCHE(NUCH).EQ.CONMOD.AND. C IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH * 50 CONTINUE * IF (NUCHA.NE.0) THEN * NUCAR = NCARAC(NUCHA) MELEME = IMAMOD * * création du nouveau maillage * NBSOUS = 0 NBREF = 0 NBELE1 = NUM(/2) NBELEM = 3* NBELE1 NBNN = NUM(/1) SEGINI IPT1 IPT1.ITYPEL = ITYPEL * DO 95 J=1, NBELE1 IPT1.ICOLOR(J) = ICOLOR(J) IPT1.ICOLOR(J+NBELE1) = ICOLOR(J) IPT1.ICOLOR(J+2*NBELE1) = ICOLOR(J) * * Recherche d'une normale a l'element courant * XNORM = 0. DO 55, K = 1, IDIM VECN(K) = 0. 55 CONTINUE IF (IDIM.EQ.2) THEN ICO1 = NUM(NBNN,J) ICO2 = NUM(1,J) DO 57, K = 1, IDIM VEC1(K) = XCOOR((ICO1-1)*(IDIM+1)+K)- C XCOOR((ICO2-1)*(IDIM+1)+K) K1 = K+1 IF (K1.GT.IDIM) K1 = 1 VECN(K) = VEC1(K1)*(-1)**K XNORM = XNORM +VECN(K)*VECN(K) 57 CONTINUE ENDIF IF (IDIM.EQ.3) THEN ICO1 = NUM(NBNN-1,J) ICO2 = NUM(NBNN,J) * DO 65 I=1, NBNN ICO3 = NUM(I,J) DO 60, K = 1, IDIM VEC1(K) = XCOOR((ICO1-1)*(IDIM+1)+K)- C XCOOR((ICO2-1)*(IDIM+1)+K) VEC2(K) = XCOOR((ICO2-1)*(IDIM+1)+K)- C XCOOR((ICO3-1)*(IDIM+1)+K) 60 CONTINUE * ICO1 = ICO2 ICO2 = ICO3 DO 65, K = 1, IDIM K1 = K+1 K2 = K+2 IF (K1.GT.IDIM) K1 = K1 - IDIM IF (K2.GT.IDIM) K2 = K2 - IDIM VECN(K) = VEC1(K1)*VEC2(K2) -VEC2(K1)*VEC1(K2) C + VECN(K) IF (I.EQ.NBNN) XNORM = XNORM + VECN(K)*VECN(K) 65 CONTINUE ENDIF XNORM = SQRT(XNORM) * DO 70, K = 1, IDIM VECN(K) = VECN(K)/XNORM 70 CONTINUE DO 95 I=1, NBNN * ICOU = NUM(I,J) IKOUR = ICPR(ICOU,NUCAR) NKON(IKOUR) = NKON(IKOUR)+1 NUIN(IKOUR) = ICOU IPT1.NUM(I,J)= NBNOE+IKOUR IPT1.NUM(I,J+NBELE1)= NBNOE+IKOUR+NNNO IPT1.NUM(I,J+2*NBELE1)= NBNOE+IKOUR+2*NNNO * * Calcul des coordonées des nouveaux points * DO 90, K = 1, IDIM XCOOR((IPT1.NUM(I,J)-1)*(IDIM+1)+K) = C XCOOR((IPT1.NUM(I,J)-1)*(IDIM+1)+K) + C VECN(K)*ICARAC.XEXCE(NUCAR) XCOOR((IPT1.NUM(I,J)+NNNO-1)*(IDIM+1)+K) = C XCOOR((IPT1.NUM(I,J)+NNNO-1)*(IDIM+1)+K) + C VECN(K)*(ICARAC.XEXCE(NUCAR)+ICARAC.XEPAI(NUCAR)/2) XCOOR((IPT1.NUM(I,J)+2*NNNO-1)*(IDIM+1)+K) = C XCOOR((IPT1.NUM(I,J)+2*NNNO-1)*(IDIM+1)+K) + C VECN(K)*(ICARAC.XEXCE(NUCAR)-ICARAC.XEPAI(NUCAR)/2) 90 CONTINUE 95 CONTINUE * * Ajustement du pointeur maillage principal * NBSOUS = IPT2.LISOUS(/1)+1 NBNN = 0 NBREF = 0 NBELEM = 0 SEGADJ IPT2 IPT2.LISOUS(NBSOUS) = IPT1 ENDIF ENDIF * 100 CONTINUE * DO 110 I=1, NNNO DO 110, K=1, IDIM XCOOR((NBNOE+I-1)*(IDIM+1)+K) = C XCOOR((NBNOE+I-1)*(IDIM+1)+K)/NKON(I) + C XCOOR((NUIN(I)-1)*(IDIM+1)+K) XCOOR((NBNOE+I+NNNO-1)*(IDIM+1)+K) = C XCOOR((NBNOE+I+NNNO-1)*(IDIM+1)+K)/NKON(I) + C XCOOR((NUIN(I)-1)*(IDIM+1)+K) XCOOR((NBNOE+I+2*NNNO-1)*(IDIM+1)+K) = C XCOOR((NBNOE+I+2*NNNO-1)*(IDIM+1)+K)/NKON(I) + C XCOOR((NUIN(I)-1)*(IDIM+1)+K) 110 CONTINUE * SEGSUP ICARAC SEGSUP NKON SEGSUP NUIN SEGSUP VECT NMAILL = IPT2.LISOUS(/1) IF (NMAILL.GE.1) THEN IF (NMAILL.EQ.1) THEN IPT3 = IPT2.LISOUS(1) SEGSUP IPT2 IPT2 = IPT3 ENDIF * * appel a PRO2 pour projeter les temperature sur le maillage * cree. isort= 1 * CALL PRO2(IPT2,IPCHT,isort,IPOUT,ilphmo) if (ierr.ne.0) return * * Recopie des valeurs du champoint dans un Chamelem image * de la geometrie initiale de la coque * mlchpo = ipout segact mlchpo * kich : pour la projection du champ de temperature on n attend qu une phase if (ichpoi(/1).ne.1) call erreur(5) MCHPOI = ICHPOI(1) SEGACT MCHPOI * * Creation du Chamelem * N1 = NMAILL N3 = 6 L1 = 12 SEGINI MCHEL1 MCHEL1.TITCHE='SCALAIRE' MCHEL1.IFOCHE=IFOUR NUCHAM = 0 * * Boucle sur l'ensemble des sous zones du modeles * DO 200, NUMO = 1, NMOD * IMODEL = KMODEL(NUMO) SEGACT IMODEL * * Test si le modele est une coque * NUFOR = NUMMFR(NEFMOD) IF (NUFOR.EQ.3 .OR. NUFOR.EQ.5 .OR. NUFOR.EQ.9) THEN * * Recherche du chamemlem de caracteristique assossiée * NUCHA = 0 DO 120, NUCH = 1, NCHAM * IF ( CONCHE(NUCH).EQ.CONMOD.AND. C IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH * 120 CONTINUE * IF (NUCHA.NE.0) THEN * NUCAR = NCARAC(NUCHA) MELEME = IMAMOD SEGACT MELEME * * création du nouveau segment MCHAML * N2 = 3 SEGINI MCHAML NUCHAM = NUCHAM+1 MCHEL1.IMACHE(NUCHAM)=MELEME MCHEL1.ICHAML(NUCHAM)=MCHAML MCHEL1.CONCHE(NUCHAM)=CONMOD MCHEL1.INFCHE(NUCHAM,1)=0 MCHEL1.INFCHE(NUCHAM,2)=0 MCHEL1.INFCHE(NUCHAM,3)=0 MCHEL1.INFCHE(NUCHAM,4)=0 MCHEL1.INFCHE(NUCHAM,5)=0 MCHEL1.INFCHE(NUCHAM,6)=1 * N1PTEL = NUM(/1) N1EL = NUM(/2) N2PTEL = 0 N2EL = 0 * DO 170, IPOS = 1, N2 * SEGINI MELVAL IF (IPOS.EQ.1) THEN NOMCHE(IPOS) = 'T' IMUL = 0 ELSEIF (IPOS.EQ.2) THEN NOMCHE(IPOS) = 'TSUP' IMUL = 1 ELSEIF (IPOS.EQ.3) THEN NOMCHE(IPOS) = 'TINF' IMUL = 2 ENDIF IELVAL(IPOS) = MELVAL TYPCHE(IPOS) = 'REAL*8' * DO 160 NUEL=1, N1EL * DO 160 NUPT=1, N1PTEL * ICO3 = NUM(NUPT,NUEL) IKOUR = ICPR(ICO3,NUCAR) * * * Boucle sur les sous-zones du champoint * DO 150, I = 1, IPCHP(/1) * MSOUPO = IPCHP(I) SEGACT MSOUPO MPOVAL = IPOVAL SEGACT MPOVAL IPT1 = IGEOC SEGACT IPT1 * * Boucle sur les composantes du champoint * DO 140, J = 1, NOCOMP(/2) * IF (NOCOMP(J).EQ.'T ') THEN * * Boucle sur les points * DO 130, K = 1, IPT1.NUM(/2) * * Comparaison des numeros de points * entre le champoint et la geometrie creee * IF (IPT1.NUM(1,K).EQ.IKOUR+NBNOE+IMUL*NNNO) C THEN VELCHE(NUPT,NUEL) = VPOCHA(K,J) GOTO 160 ENDIF * 130 CONTINUE ENDIF 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE ENDIF ENDIF * 200 CONTINUE * * Suppression du champoint * DO 220, I = 1, IPCHP(/1) * MSOUPO = IPCHP(I) MPOVAL = IPOVAL IPT1 = IGEOC ***** SEGSUP IPT1 SEGSUP MPOVAL SEGSUP MSOUPO * 220 CONTINUE SEGSUP MCHPOI * * Suppression du maillage intermediaire * SEGACT IPT2 * DO 240, IOB =1, IPT2.LISOUS(/1) * IPT1 = IPT2.LISOUS(IOB) SEGSUP IPT1 * 240 CONTINUE ***** SEGSUP IPT2 * * Reajustement du tableau de coordonées * NBPTS = NBNOE SEGADJ MCOORD * * RESTITUTION DU CHAMP DE SORTIE * ITPR= MCHEL1 * ELSE CALL ERREUR(704) ENDIF * SEGSUP ICPR SEGSUP NCARAC END