C CALJBR SOURCE MAGN 08/10/10 21:15:00 6185 SUBROUTINE CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, *NES,ND,NP,NPG,IAXI,AIRE,AJ,SGN) C*********************************************************************** C Ce Sp calcule le gradient des fonctions de formes pour un element C courant a partir des gradients de l'element de reference et des C coordonnees des points de l'element C Il calcule le produit PgSq poids d'integration x element d'aire et C le Jacobien C Pour les elements coques ou filaires les normales et tangentes C sont donnees a la place du Jacobien C les fonctions de forme FN sont donnees en entree pour calculer C les elements d'aire en axisymetrique C C Entree : C NES Dimension d'espace de l'element de reference C (different de ND pour les elements coques ou filaires) C ND Dimension d'espace C NP Nombre de noeuds de l'element C NPG Nombre de points d'integration C IAXI =0 en 3D et en 2D PLAN C IAXI NE 0 en axi symetrique axe de symetrie x C FN(NP,NPG) Valeurs des fonctions de forme aux points d'integration C GR(NES,NP,NPG) Valeurs des gradients dans l'element de reference C PG(NPG) Poids d'integration C XYZ(ND,NP) Coordonnees des noeuds de l'element C C Sortie C HR(NES,NP,NPG) Valeurs des gradients dans l'element courant C PGSQ(NPG) poids d'integration x element d'aire C ATTENTION en axi : poids d'integration x element d'aire x 2 pi R C RPG(NPG) Rayon des points d'integration C AIRE Aire de l'element C AJ(ND,ND,NPG) Valeurs du jacobien pour chaque point d'integration C SGN Pour les elements massif =1 si meme orientation element de Ref C SGN Pour les elements massif =-1 si orientation opposee C SGN 0 pour les autres types d'elements C C C CE SP TRAITE LES ELEMENTS VOLUMIQUES, SURFACIQUES et FILAIRES C DANS LES CAS 2D (PLAN et AXI) ET 3D C C Dans AJ on stoke l'inverse du Jacobien (AJ=1/J) si l'element est C volumique C Sinon pour les elements coques on stoke tangentes et normales C AJ=|tx ty| |tx ty tz| C |nx ny| ou |ux uy uz| C |nx ny nz| C C C CALCUL INTERMEDIAIRE DE L'ELEMENT D'AIRE SQ=DET(J) C C*********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C REAL*8 FN(NP,NPG),GR(NES,NP,NPG),HR(ND,NP,NPG) REAL*8 PG(NPG),XYZ(ND,NP),PGSQ(NPG),RPG(NPG),AJ(ND,ND,NPG) REAL*8 AG(3,3,25) -INC CCREEL C C*** SGN=0.D0 CALL INITD(RPG,NPG,1.D0) IF(NES.EQ.1.AND.ND.EQ.2)THEN IF(NPG.GT.25)CALL ARRET(0) AIRE=0.D0 DO 110 L=1,NPG AJX=0.D0 AJY=0.D0 DO 111 I=1,NP AJX=AJX+GR(1,I,L)*XYZ(1,I) AJY=AJY+GR(1,I,L)*XYZ(2,I) 111 CONTINUE AJN=(AJX*AJX+AJY*AJY)**0.5D0 C write(6,*)' AJX,AJY=',AJX,AJY,AJN DET=AJX*AJX+AJY*AJY AJ(1,1,L)=AJX/AJN AJ(2,1,L)=AJY/AJN AJ(1,2,L)=-AJY/AJN AJ(2,2,L)=AJX/AJN AG(1,1,L)=AJX/DET AG(2,1,L)=-AJY/DET AG(1,2,L)=AJY/DET AG(2,2,L)=AJX/DET PGSQ(L)=PG(L)*AJN AIRE=AIRE+PGSQ(L) 110 CONTINUE C write(6,*)' AJ ds CALJBR',AIRE C write(6,1002)AJ DO 131 L=1,NPG DO 131 I=1,NP DO 131 N=1,ND U=0.D0 DO 132 M=1,NES 132 U=U+AG(M,N,L)*GR(M,I,L) 131 HR(N,I,L)=U C write(6,*)'GR' C write(6,1002)gr C write(6,*)'HR' C write(6,1002)hr ELSEIF(NES.EQ.1.AND.ND.EQ.3)THEN IF(NPG.GT.25)CALL ARRET(0) AIRE=0.D0 DO 310 L=1,NPG AJX=0.D0 AJY=0.D0 AJZ=0.D0 BJX=0.D0 BJY=0.D0 BJZ=0.D0 DO 311 I=1,NP AJX=AJX+GR(1,I,L)*XYZ(1,I) AJY=AJY+GR(1,I,L)*XYZ(2,I) AJZ=AJZ+GR(1,I,L)*XYZ(3,I) C BJX=BJX+GR(2,I,L)*XYZ(1,I) C BJY=BJY+GR(2,I,L)*XYZ(2,I) C BJZ=BJZ+GR(2,I,L)*XYZ(3,I) 311 CONTINUE C write(6,*)ajx,ajy,ajz IF(ABS(AJX).NE.0.D0.AND.ABS(AJY).NE.0.D0)THEN BJX=AJY BJY=-AJX BJZ=AJZ ELSE BJX=1.D0 BJY=1.D0 BJZ=0.D0 ENDIF XB=AJY*BJZ-AJZ*BJY YB=AJZ*BJX-AJX*BJZ ZB=AJX*BJY-AJY*BJX AJN=(XB*XB+YB*YB*ZB*ZB)**0.5D0+1.D-5 C write(6,*)' AJN=',AJN PGSQ(L)=PG(L)*AJN AIRE=AIRE+PGSQ(L) AJ(1,1,L)=AJX/AJN AJ(2,1,L)=AJY/AJN AJ(3,1,L)=AJZ/AJN AJ(1,2,L)=BJX/AJN AJ(2,2,L)=BJY/AJN AJ(3,2,L)=BJZ/AJN AJ(1,3,L)=XB/AJN AJ(2,3,L)=YB/AJN AJ(3,3,L)=ZB/AJN AG(1,1,L)=AJX AG(2,1,L)=AJY AG(3,1,L)=AJZ AG(1,2,L)=BJX AG(2,2,L)=BJY AG(3,2,L)=BJZ AG(1,3,L)=XB AG(2,3,L)=YB AG(3,3,L)=ZB D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L) D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L) D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L) D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L) D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L) D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L) D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L) D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L) D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L) VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31 RVINT=1.D0/VINT AG(1,1,L)= RVINT*D11 AG(1,2,L)=-RVINT*D12 AG(1,3,L)= RVINT*D13 AG(2,1,L)=-RVINT*D21 AG(2,2,L)= RVINT*D22 AG(2,3,L)=-RVINT*D23 AG(3,1,L)= RVINT*D31 AG(3,2,L)=-RVINT*D32 AG(3,3,L)= RVINT*D33 310 CONTINUE C write(6,*)' AJ ds CALJBR' C write(6,1002)AJ DO 331 L=1,NPG DO 331 I=1,NP DO 331 N=1,ND U=0.D0 DO 332 M=1,NES 332 U=U+AG(M,N,L)*GR(M,I,L) 331 HR(N,I,L)=U C write(6,*)'GR' C write(6,1002)gr C write(6,*)'HR' C write(6,1002)hr ELSEIF(NES.EQ.2.AND.ND.EQ.3)THEN C write(6,*)' Je passe ICI' AIRE=0.D0 DO 210 L=1,NPG AJX=0.D0 AJY=0.D0 AJZ=0.D0 BJX=0.D0 BJY=0.D0 BJZ=0.D0 DO 211 I=1,NP AJX=AJX+GR(1,I,L)*XYZ(1,I) AJY=AJY+GR(1,I,L)*XYZ(2,I) AJZ=AJZ+GR(1,I,L)*XYZ(3,I) BJX=BJX+GR(2,I,L)*XYZ(1,I) BJY=BJY+GR(2,I,L)*XYZ(2,I) BJZ=BJZ+GR(2,I,L)*XYZ(3,I) 211 CONTINUE XB=AJY*BJZ-AJZ*BJY YB=AJZ*BJX-AJX*BJZ ZB=AJX*BJY-AJY*BJX AJN=(XB*XB+YB*YB+ZB*ZB)**0.5D0 PGSQ(L)=PG(L)*AJN AIRE=AIRE+PGSQ(L) AJ(1,1,L)=AJX/AJN AJ(2,1,L)=AJY/AJN AJ(3,1,L)=AJZ/AJN AJ(1,2,L)=BJX/AJN AJ(2,2,L)=BJY/AJN AJ(3,2,L)=BJZ/AJN AJ(1,3,L)=XB/AJN AJ(2,3,L)=YB/AJN AJ(3,3,L)=ZB/AJN AG(1,1,L)=AJX AG(2,1,L)=AJY AG(3,1,L)=AJZ AG(1,2,L)=BJX AG(2,2,L)=BJY AG(3,2,L)=BJZ AG(1,3,L)=XB AG(2,3,L)=YB AG(3,3,L)=ZB D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L) D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L) D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L) D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L) D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L) D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L) D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L) D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L) D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L) VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31 RVINT=1.D0/VINT AG(1,1,L)= RVINT*D11 AG(1,2,L)=-RVINT*D12 AG(1,3,L)= RVINT*D13 AG(2,1,L)=-RVINT*D21 AG(2,2,L)= RVINT*D22 AG(2,3,L)=-RVINT*D23 AG(3,1,L)= RVINT*D31 AG(3,2,L)=-RVINT*D32 AG(3,3,L)= RVINT*D33 210 CONTINUE C DO 231 L=1,NPG DO 231 I=1,NP DO 231 N=1,ND U=0.D0 DO 232 M=1,NES 232 U=U+AG(M,N,L)*GR(M,I,L) 231 HR(N,I,L)=U C write(6,*)'GR' C write(6,1002)gr C write(6,*)'HR' C write(6,1002)hr ELSE DO 10 L=1,NPG DO 10 M=1,ND DO 10 N=1,ND AJT=0.D0 DO 11 I=1,NP AJT=AJT+GR(M,I,L)*XYZ(N,I) 11 CONTINUE AJ(N,M,L)=AJT 10 CONTINUE C DO 20 L=1,NPG IF(ND.EQ.1)THEN VINT=AJ(1,1,L) C VINT=ABS(VINT) AJ(1,1,L)=1.D0/VINT ELSEIF(ND.EQ.2)THEN VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L) C VINT=ABS(VINT) RVINT=1.D0/VINT D11=AJ(2,2,L) D12=AJ(1,2,L) D21=AJ(2,1,L) D22=AJ(1,1,L) AJ(1,1,L)= RVINT*D11 AJ(1,2,L)=-RVINT*D12 AJ(2,1,L)=-RVINT*D21 AJ(2,2,L)= RVINT*D22 ELSEIF(ND.EQ.3)THEN D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L) D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L) D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L) D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L) D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L) D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L) D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L) D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L) D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L) VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31 C VINT=ABS(VINT) RVINT=1.D0/VINT AJ(1,1,L)= RVINT*D11 AJ(1,2,L)=-RVINT*D12 AJ(1,3,L)= RVINT*D13 AJ(2,1,L)=-RVINT*D21 AJ(2,2,L)= RVINT*D22 AJ(2,3,L)=-RVINT*D23 AJ(3,1,L)= RVINT*D31 AJ(3,2,L)=-RVINT*D32 AJ(3,3,L)= RVINT*D33 ELSE VINT=0.D0 ENDIF PGSQ(L)=VINT C*** CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER) C 20 CONTINUE C DO 31 L=1,NPG DO 31 I=1,NP DO 31 N=1,ND U=0.D0 DO 32 M=1,ND 32 U=U+AJ(M,N,L)*GR(M,I,L) 31 HR(N,I,L)=U U=0.D0 V=0.D0 DO 4 L=1,NPG W=ABS(PGSQ(L)) U=U+PG(L)*W V=V+PG(L)*PGSQ(L) 4 PGSQ(L)=PG(L)*W AIRE=U SGN=SIGN(1.D0,V) ENDIF IF(IAXI.EQ.0)RETURN C DO 45 L=1,NPG RPGT =0.D0 DO 46 I=1,NP RPGT=RPGT+XYZ(1,I)*FN(I,L) 46 CONTINUE RPG(L)=RPGT 45 CONTINUE C DEUPI=2.D0*XPI U=0.D0 DO 47 L=1,NPG C? PGSQ(L)=PGSQ(L)*DEUPI*RPG(L) U=U+PGSQ(L)*DEUPI*RPG(L) 47 CONTINUE AIRE=U C RETURN 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END