C KALPRE SOURCE OF166741 24/10/03 21:15:23 12022 SUBROUTINE KALPRE (MYMOD,INFOEL,SKFACE,XDEC,NELD) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C---------------------------------------------------------------------- C Facteurs de forme OPTION CACHE 3D C Calcul du decoupage des faces et initialisations C SP appelle par facgen C entree: C MYMOD : pointeur de l'objet modèle C INFOEL : informations concernant le type des éléments des maillages C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF C SHC3D : pointeur sur la surface de projection C XDEC : paramŠtre pour le decoupage des faces C NELD : nombre d'éléments virtuels C (un élément COQ -> 2 facteurs de forme <=> 2 éléments virtuels) C sortie: C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF C SHC3D : pointeur sur la surface de projection C---------------------------------------------------------------------- -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMMODEL -INC SMELEME -INC TFFOR3D C ----------------------------------- SEGMENT SDECOU INTEGER KDECOU(MFACE) ENDSEGMENT POINTEUR MYMOD.MMODEL C ----------------------------------- C Stockage d'informations concernant le type des éléments des maillages SEGMENT ,INFOEL LOGICAL KCOQ(N1),KQUAD(N1) ENDSEGMENT C ----------------------------------- DIMENSION A(3,4),UU(4),IN(3) DIMENSION A1(3,3),A2(3,3),U1(4),U2(4),G1(3),G2(3) DIMENSION GG(3,400),SS(400) LOGICAL ICOQ C C-------------------------------------------------------------------- C C DMIN=1.E-5 DMAX=1./DMIN C... nombre max de triangles de subdivision NPM=20 C... quadratiques NSPA=1 IF(INFOEL.EQ.0) THEN ICOQ=.FALSE. ELSE ICOQ=.TRUE. SEGACT INFOEL ENDIF C C>>> CREATION DE L'OBJET FACE C NFACE = 0 NELD = 0 C MTYP = MYMOD.KMODEL(/1) SEGACT INFOEL C C On va compter le nombre d'éléments virtuels (égal au nombre de C facteurs de forme par ligne) : NELD , et le nombre de faces : NFACE . DO 10 ITYP=1,MTYP IMODE1 = MYMOD.KMODEL(ITYP) IPT1 = IMODE1.IMAMOD NEL = IPT1.NUM(/2) NPGEO = IPT1.NUM(/1) C C un élément COQ -> 2 facteurs de formes IF (KCOQ(ITYP)) NEL=2*NEL NELD = NELD + NEL C C un élément carré -> 2 éléments triangles IF(NPGEO.EQ.4.OR.NPGEO.EQ.8) NEL = 2*NEL NFACE=NFACE+NEL 10 CONTINUE C IF (KIMP.GE.1) THEN WRITE( 6,*) ' NOMBRE TOTAL D ELEMENTS ',NELD,' DE FACES ',NFACE ENDIF C C MFACE = NFACE MS = 3 SEGINI SKFACE SEGINI SDECOU C C KEL : INDICE ELEMENT DECOUPAGE (-> numéro global de la facette) C IEL : INDICE ELEMENT D'ORIGINE (-> numero global de l'élément virtuel) C Ces indices sont différents car on découpe les quadrangles en triangles C KEL = 1 IEL = 0 C DO 100 ITYP = 1, MTYP IMODE1 = MYMOD.KMODEL(ITYP) IPT1 = IMODE1.IMAMOD NSGEO = IPT1.NUM(/1) IF(ICOQ.AND.KQUAD(ITYP)) THEN NS = NSGEO/2 ELSE NS=NSGEO ENDIF IF(ICOQ.AND.KQUAD(ITYP)) NSPA=2 NEL = IPT1.NUM(/2) C C... TRI3 ou TRI6 IF (NS.EQ.3) THEN *//////// * WRITE (6,*) 'Le maillage ',IPT1, * # 'est constitué d éléments triangles' *//////// DO 110 I = 1,NEL IF (KCOQ(ITYP)) THEN IEL1 = IEL + (2*I-1) ELSE IEL1 = IEL + I ENDIF KCORR(KEL) = IEL1 DO 111 IS = 1,NSGEO,NSPA LS=IS IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 C** NUK(IS,KEL) = IPT1.NUM(IS,I) NUK(LS,KEL) = IPT1.NUM(IS,I) IREF = (IDIM+1)*(NUK(LS,KEL)-1) DO 112 K = 1,KES C** A(K,LS) = XCOOR(IREF+K) A(K,LS) = XCOOR(IREF+K) 112 CONTINUE *//////// * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES) *//////// 111 CONTINUE IF (KIMP.GE.3) THEN WRITE (6,*) 'La facette ',KEL,'repose sur les points', # (IPT1.NUM(IS,I),IS=1,NS) ENDIF CALL KNORM(KES,A,NS,UU,S(KEL)) DO 113 L = 1,4 U(L,KEL) = UU(L) 113 CONTINUE C WRITE(6,*) ' KEL S U ',KEL,S(KEL),(U(K,NEL),K=1,4) KEL = KEL + 1 *//////// IF (KCOQ(ITYP)) THEN C On remplit KCOOR , NUK , U et S pour l'élément inverse KCORR(KEL) = IEL1 + 1 NUK(1,KEL) = NUK(3,KEL-1) NUK(2,KEL) = NUK(2,KEL-1) NUK(3,KEL) = NUK(1,KEL-1) DO 114 L = 1,4 U(L,KEL) = - UU(L) 114 CONTINUE S(KEL) = S(KEL-1) KEL = KEL + 1 ENDIF C 110 CONTINUE C IF (KCOQ(ITYP)) THEN IEL = IEL + NEL*2 ELSE IEL = IEL + NEL ENDIF C C... QUA4 ou QUA8 ELSE *//////// * WRITE (6,*) 'Le maillage ',IPT1, * # 'est constitué d éléments carrés' *//////// C DO 120 I = 1,NEL C IF (KCOQ(ITYP)) THEN IEL1 = IEL + (2*I-1) ELSE IEL1 = IEL + I ENDIF DO 121 IS = 1,NSGEO,NSPA LS=IS IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 IREF = (IDIM+1)*(IPT1.NUM(IS,I)-1) *//////// * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) * WRITE (6,*) 'IREF =',IREF *//////// DO 122 K = 1,KES C** A(K,IS) = XCOOR(IREF+K) A(K,LS) = XCOOR(IREF+K) 122 CONTINUE *//////// * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES) *//////// 121 CONTINUE C C On regarde comment on va découper le quadrangle DIA1 = 0. DO 123 K = 1,KES DIA1 = DIA1 + (A(K,1)-A(K,3))**2 123 CONTINUE DIA2 = 0. DO 124 K = 1,KES DIA2 = DIA2 + (A(K,2)-A(K,4))**2 124 CONTINUE *////// * WRITE(6,*) ' DIA1 DIA2 ',DIA1,DIA2 *////// C stockage des connectivités des traiangles IF(ICOQ.AND.KQUAD(ITYP)) THEN * cas des qua8 IF (DIA1.LE.DIA2) THEN NUK(1,KEL) = IPT1.NUM(1,I) NUK(2,KEL) = IPT1.NUM(3,I) NUK(3,KEL) = IPT1.NUM(5,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(5,I) NUK(2,KEL) = IPT1.NUM(3,I) NUK(3,KEL) = IPT1.NUM(1,I) KEL = KEL + 1 ENDIF NUK(1,KEL) = IPT1.NUM(5,I) NUK(2,KEL) = IPT1.NUM(7,I) NUK(3,KEL) = IPT1.NUM(1,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(1,I) NUK(2,KEL) = IPT1.NUM(7,I) NUK(3,KEL) = IPT1.NUM(5,I) KEL = KEL + 1 ENDIF ELSE NUK(1,KEL) = IPT1.NUM(3,I) NUK(2,KEL) = IPT1.NUM(5,I) NUK(3,KEL) = IPT1.NUM(7,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(7,I) NUK(2,KEL) = IPT1.NUM(5,I) NUK(3,KEL) = IPT1.NUM(3,I) KEL = KEL + 1 ENDIF NUK(1,KEL) = IPT1.NUM(7,I) NUK(2,KEL) = IPT1.NUM(1,I) NUK(3,KEL) = IPT1.NUM(3,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(3,I) NUK(2,KEL) = IPT1.NUM(1,I) NUK(3,KEL) = IPT1.NUM(7,I) KEL = KEL + 1 ENDIF ENDIF ELSE * cas des qua4 IF (DIA1.LE.DIA2) THEN NUK(1,KEL) = IPT1.NUM(1,I) NUK(2,KEL) = IPT1.NUM(2,I) NUK(3,KEL) = IPT1.NUM(3,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(3,I) NUK(2,KEL) = IPT1.NUM(2,I) NUK(3,KEL) = IPT1.NUM(1,I) KEL = KEL + 1 ENDIF NUK(1,KEL) = IPT1.NUM(3,I) NUK(2,KEL) = IPT1.NUM(4,I) NUK(3,KEL) = IPT1.NUM(1,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(1,I) NUK(2,KEL) = IPT1.NUM(4,I) NUK(3,KEL) = IPT1.NUM(3,I) KEL = KEL + 1 ENDIF ELSE NUK(1,KEL) = IPT1.NUM(2,I) NUK(2,KEL) = IPT1.NUM(3,I) NUK(3,KEL) = IPT1.NUM(4,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(4,I) NUK(2,KEL) = IPT1.NUM(3,I) NUK(3,KEL) = IPT1.NUM(2,I) KEL = KEL + 1 ENDIF NUK(1,KEL) = IPT1.NUM(4,I) NUK(2,KEL) = IPT1.NUM(1,I) NUK(3,KEL) = IPT1.NUM(2,I) KEL = KEL + 1 IF (KCOQ(ITYP)) THEN NUK(1,KEL) = IPT1.NUM(2,I) NUK(2,KEL) = IPT1.NUM(1,I) NUK(3,KEL) = IPT1.NUM(4,I) KEL = KEL + 1 ENDIF ENDIF ENDIF C DO 125 IL = 1,2 C IF (KCOQ(ITYP)) THEN INDICE = KEL - IL*2 ELSE INDICE = KEL - IL ENDIF KCORR(INDICE) = IEL1 C DO 126 IS = 1,3 IREF = (IDIM+1)*(NUK(IS,INDICE)-1) DO 127 K = 1,KES A(K,IS) = XCOOR(IREF+K) C WRITE(6,*) 'A ',K,IS,A(K,IS) 127 CONTINUE 126 CONTINUE CALL KNORM(KES,A,3,UU,S(INDICE)) DO 128 L = 1,4 U(L,INDICE) = UU(L) 128 CONTINUE CC WRITE(6,*) ' KEL S U ',INDICE,S(INDICE),(U(K,INDICE),K=1,4) 125 CONTINUE C IF (KCOQ(ITYP)) THEN C On remplit KCOOR , U et S pour les éléments inverses DO 135 IL = 1,2 INDICE = KEL - IL*2 +1 KCORR(INDICE) = IEL1 +1 DO 138 L = 1,4 U(L,INDICE) = - U(L,INDICE-1) 138 CONTINUE S(INDICE)=S(INDICE-1) 135 CONTINUE ENDIF C 120 CONTINUE C IF (KCOQ(ITYP)) THEN IEL = IEL + NEL*2 ELSE IEL = IEL + NEL ENDIF C C ENDIF C 100 CONTINUE C C On a pris toute l'information nécessaire concernant le maillage SEGDES INFOEL *//////// IF(KIMP.GE.2) THEN WRITE (6,*) 'Subdivision' ENDIF *//////// C----------------------------------------------------------- DO 200 K1=1,MFACE DO 202 K = 1,KES G1(K)=0. DO 201 ISS = 1,MS IREF = (IDIM+1)*(NUK(ISS,K1)-1) A1(K,ISS) = XCOOR(IREF+K) G1(K) = G1(K) + A1(K,ISS) 201 CONTINUE G1(K)=G1(K)/3. 202 CONTINUE DO 203 K=1,KES+1 U1(K) = U(K,K1) 203 CONTINUE C IF (XDEC.GE.0.01) THEN DK1 = DMAX DO 400 K2 = 1,MFACE C DO 412 K=1,KES+1 U2(K) = U(K,K2) 412 CONTINUE DO 212 K = 1,KES G2(K)=0. DO 211 ISS = 1,MS IREF = (IDIM+1)*(NUK(ISS,K2)-1) A2(K,ISS) = XCOOR(IREF+K) G2(K)=G2(K)+A2(K,ISS) 211 CONTINUE G2(K)=G2(K)/3. 212 CONTINUE C>>> VISIBILITE A PRIORI C CALL KPRIOR(KES,MS,MS,A1,A2,U1,U2,KVU) IF (KVU.NE.0) THEN D1=G1(1)-G2(1) D2=G1(2)-G2(2) D3=G1(3)-G2(3) DKK1 = SQRT(D1*D1+D2*D2+D3*D3) C WRITE(6,*) ' DKK1 CK1 ',DKK1,CK1 IF(DKK1.GE.DMIN) THEN CK1 = ABS(U1(1)*D1+ U1(2)*D2 + U1(3)*D3)/DKK1 IF(CK1.GE.0.01) THEN DK1 =MIN(DK1,DKK1) ENDIF ENDIF ENDIF 400 CONTINUE DR = DK1/XDEC ELSE DR = DMAX ENDIF C 901 CONTINUE CALL KCREPA(DR,A1,S(K1),KES,MS,NPAT,GG,SS) IF(NPAT.GT.NPM) THEN DR=DR*2. IF (KIMP.GE.2) THEN WRITE(6,900) K1,NPAT 900 FORMAT(1X,' FACE ',I4,' : SUBDIVISION EXCESSIVE ',I4) ENDIF GOTO 901 ENDIF NPATCH = NPAT KDECOU(K1)=NPAT MSP = KES SEGINI IPATCH DO 130 IP = 1,NPAT SP(IP) = SS(IP)/S(K1) DO 129 K=1,KES GP(K,IP) = GG(K,IP) 129 CONTINUE C WRITE(6,*) ' SP GP ',SP(IP),GP(1,IP),GP(2,IP),GP(3,IP) 130 CONTINUE KPATCH(K1) = IPATCH SEGDES IPATCH 200 CONTINUE C----------------------------------------------------------- IF(KIMP.GE.2) THEN WRITE(6,*) WRITE(6,*) 'NOMBRE DE SUBDIVISIONS PAR FACE ' WRITE(6,1000) (KDECOU(I),I=1,MFACE) 1000 FORMAT(1X,10(I4)) ENDIF C SEGDES SKFACE SEGSUP SDECOU C RETURN END