C JACONO SOURCE CB215821 20/11/04 21:18:27 10766 SUBROUTINE JACONO(IPMODL,INORM,IPCHE,IRET) C======================================================================= C ENTREES : C --------- C IPMODL= pointeur sur un MMODEL C INORM = 1 si les vecteurs doivent etre normes 0 sinon C C SORTIES : C -------- C C IPCHE = CHAMELEM contenant les JACOBIENS C ( = normale aux faces des elements dans le cas des coques) C ( = tangente a la fibre neutre dans le cas des poutres) C IRET = 1 si succes 0 sinon C C Passage au nouveau Chamelem PAR S.RAMAHANDRY le 11/09/90 C C 2013-01-02 (BP) : ajout zones cohesives (ZCO2,3 et 4 => coque mince) C + calcul de la tangente pour les poutres C C C===================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(UN=1.D0,XZER=0.D0) -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC SMCHAML -INC SMMODEL -INC SMELEME -INC SMCOORD -INC SMINTE C SEGMENT TRA REAL*8 XEL(3,NBNN) ,SHP(6,NBNN) ,XE(3,NBNN) ENDSEGMENT C SEGMENT TR1 REAL*8 TH(NBN1) ,TXR(3,3,NBN1) ,XJ(3,3) ENDSEGMENT C SEGMENT INFO INTEGER INFELL(JG) ENDSEGMENT C SEGMENT MPTVAL INTEGER IPOS(NS) , NSOF(NS) INTEGER IVAL(NCOSOU) CHARACTER*16 TYVAL(NCOSOU) ENDSEGMENT C DIMENSION BPSS(3,3) CHARACTER*8 CMATE C REAL*8 XU(3) REAL*8 XV(3) REAL*8 XW(3) C NHRM=NIFOUR IRET=1 C C ACTIVATION DU MODELE C MMODEL= IPMODL SEGACT MMODEL NSOUS=KMODEL(/1) C C CREATION DU MCHELM C N1=NSOUS L1=8 N3=6 SEGINI MCHELM IPCHE=MCHELM IF (INORM .EQ. 1) THEN TITCHE='NORMALES' ELSE TITCHE='VECTEURS SURFACE' ENDIF IFOCHE=IFOUR IRT8=1 C____________________________________________________________________ C C DEBUT DE LA BOUCLE SUR LES DIFFERENTES ZONES C____________________________________________________________________ C DO 500 ISOUS=1,NSOUS C C ON RECUPERE L INFORMATION GENERALE C IMODEL=KMODEL(ISOUS) SEGACT IMODEL IPMAIL=IMAMOD IMACHE(ISOUS)=IPMAIL CONCHE(ISOUS)=CONMOD C C TRAITEMENT DU MODELE C MELE=NEFMOD MELEME=IMAMOD NFOR=FORMOD(/2) NMAT=MATMOD(/2) C____________________________________________________________________ C C INFORMATION SUR L'ELEMENT FINI C____________________________________________________________________ C if(infmod(/1).lt.7) then CALL ELQUOI(MELE,0,5,IPINF,IMODEL) IF (IERR.NE.0) THEN SEGSUP MCHELM IRET=0 RETURN ENDIF INFO=IPINF MELE =INFELL(1) MFR =INFELL(13) MINTE=INFELL(11) MINTE1=INFELL(12) segsup info else MELE =INFELE(1) MFR =INFELE(13) MINTE=INFMOD(7) MINTE1=INFMOD(8) endif C INFCHE(ISOUS,1)=0 INFCHE(ISOUS,2)=0 INFCHE(ISOUS,3)=NHRM INFCHE(ISOUS,4)=MINTE INFCHE(ISOUS,5)=0 INFCHE(ISOUS,6)=5 C C INITIALISATION DE MINTE C SEGACT MINTE NBPGAU=POIGAU(/1) C C ACTIVATION DU MELEME C SEGACT MELEME NBNN =NUM(/1) NBELEM=NUM(/2) C C RECHERCHE DE LA TAILLE DES MELVAL A ALLOUER C IF (MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN N1PTEL=NBPGAU N1EL=NBELEM ELSEIF(MFR.EQ.7 .OR. MFR.EQ.13) THEN N1PTEL=NBPGAU N1EL=NBELEM ELSE N1PTEL = 0 N1EL = 0 ENDIF C C CREATION DU MCHAML DE LA SOUS ZONE C NJAC=IDIM N2 = IDIM SEGINI MCHAML ICHAML(ISOUS)=MCHAML NS=1 NCOSOU=NJAC SEGINI MPTVAL IVAJAC=MPTVAL C C 2 OU 3 COMPOSANTES C ICOMP=1 IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN NOMCHE(ICOMP)='VR ' ELSE NOMCHE(ICOMP)='VX ' ENDIF TYPCHE(ICOMP)='REAL*8' N2PTEL=0 N2EL=0 SEGINI MELVA1 IELVAL(ICOMP)=MELVA1 IVAL(ICOMP)=MELVA1 C ICOMP=2 IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN NOMCHE(ICOMP)='VZ ' ELSE NOMCHE(ICOMP)='VY ' ENDIF TYPCHE(ICOMP)='REAL*8' N2PTEL=0 N2EL=0 SEGINI MELVA2 IELVAL(ICOMP)=MELVA2 IVAL(ICOMP)=MELVA2 C IF (IDIM .EQ. 3) THEN ICOMP=3 NOMCHE(ICOMP)='VZ ' TYPCHE(ICOMP)='REAL*8' N2PTEL=0 N2EL=0 SEGINI MELVA3 IELVAL(ICOMP)=MELVA3 IVAL(ICOMP)=MELVA3 ENDIF C 44 CONTINUE C SEGINI TRA C C ================ FORMULATION MASSIVE ======================= C C IF(MFR.EQ.1.OR.MFR.EQ.33) THEN GOTO 520 C C ================ FORMULATION COQUE MINCE ===================== C C ELSE IF(MFR.EQ.3.OR.MFR.EQ.9 .OR. MFR.EQ.77) THEN IDI2=IDIM-1 DO 3000 IB=1,NBELEM *--------------Calcul de la normale a l'élément IREF1 = (IDIM + 1)*(MELEME.NUM(1, IB) - 1) IREF2 = (IDIM + 1)*(MELEME.NUM(2, IB) - 1) DO 28 IC = 1, IDIM 28 XU(IC) = XCOOR(IREF2+IC)-XCOOR(IREF1+IC) XNORU = 0. DO 31 IC = 1, IDIM 31 XNORU = XNORU + XU(IC)*XU(IC) XNORU = SQRT(XNORU) DO 32 IC = 1, IDIM 32 XU(IC) = XU(IC)/XNORU IF (IDIM .EQ. 2) THEN XW(1) = -XU(2) XW(2) = XU(1) ELSE IN = 3 33 IREF3 = (IDIM + 1)*(MELEME.NUM(IN, IB) - 1) DO 34 IC = 1, IDIM 34 XV(IC) = XCOOR(IREF3+IC)-XCOOR(IREF1+IC) XNORV = 0. DO 35 IC = 1, IDIM 35 XNORV = XNORV + XV(IC)*XV(IC) XNORV = SQRT(XNORV) DO 36 IC = 1, IDIM 36 XV(IC) = XV(IC)/XNORV XW(1) = XU(2)*XV(3) - XU(3)*XV(2) XW(2) = XU(3)*XV(1) - XU(1)*XV(3) XW(3) = XU(1)*XV(2) - XU(2)*XV(1) XNORW = 0. DO 37 IC = 1, IDIM 37 XNORW = XNORW + XW(IC)*XW(IC) IF (XNORW .LT. 1.E-4) THEN IN = IN + 1 if(IN.le.NBNN) GOTO 33 ENDIF XNORW = SQRT(XNORW) IF(XNORW .LT.1.E-4) call erreur(345) DO 38 IC = 1, IDIM 38 XW(IC) = XW(IC)/XNORW ENDIF *--------------Fin du calcul de la normale a l'élément CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C IF(IDIM.EQ.2)THEN CALL VPAST2(XE,BPSS) ELSE IF(IDIM.EQ.3) THEN CALL VPAST(XE,BPSS) ENDIF CALL VCORL1(XE,XEL,BPSS,NBNN) DO 3002 IC=1,NBPGAU IF (INORM .EQ. 1) THEN DJAC = 1. ELSE DO 3003 ID=1,6 DO 3003 IE=1,NBNN SHP(ID,IE)=SHPTOT(ID,IE,IC) 3003 CONTINUE CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC) ENDIF MPTVAL=IVAJAC MELVAL = IVAL(1) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1) MELVAL = IVAL(2) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2) IF (IDIM .EQ. 3) THEN MELVAL = IVAL(3) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3) ENDIF 3002 CONTINUE 3000 CONTINUE GOTO 520 C C ================ FORMULATION POUTRE ET TUYAU ==================== C C ELSE IF(MFR.EQ.7.OR.MFR.EQ.13) THEN IDI2=IDIM-1 DO 4000 IB=1,NBELEM *-----------Calcul de la tangente a l'élément IREF1 = (IDIM + 1)*(MELEME.NUM(1, IB) - 1) IREF2 = (IDIM + 1)*(MELEME.NUM(2, IB) - 1) DO 41 IC = 1, IDIM 41 XU(IC) = XCOOR(IREF2+IC)-XCOOR(IREF1+IC) XNORU = 0.D0 DO 42 IC = 1, IDIM 42 XNORU = XNORU + XU(IC)*XU(IC) XNORU = SQRT(XNORU) DO 43 IC = 1, IDIM 43 XU(IC) = XU(IC)/XNORU *-----------Fin du calcul de la tangente a l'élément *-----------Calcul de la tangente a l'élément en chaque point de Gauss c BP : On suppose le jacobien constant dans l'element (idem POUJAC.eso) C => on sort le calcul du jacobien de la boucle sur les points de Gauss. C Cela implique que la POUTre de Bernoulli n'est pas isoparamétrique... IF (INORM .EQ. 1) THEN DJAC = 1.D0 ELSE DJAC = 1.D0/DBLE(NBPGAU) ENDIF DO 4002 IC=1,NBPGAU MPTVAL=IVAJAC MELVAL = IVAL(1) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(1) MELVAL = IVAL(2) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(2) IF (IDIM .EQ. 3) THEN MELVAL = IVAL(3) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC,MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XU(3) ENDIF 4002 CONTINUE 4000 CONTINUE GOTO 520 C C ================ FORMULATION COQUE EPAISSE ==================== C C ELSE IF(MFR.EQ.5) THEN SEGACT MINTE1 NBPGA1=MINTE1.POIGAU(/1) NBN1 =MINTE1.SHPTOT(/2) SEGINI TR1 C C UNE PETITE HORREUR ON CONSIDERE LES EPAISSEURS CONSTANTES C DO 5010 IC=1,NBNN TH(IC)=UN 5010 CONTINUE DO 5000 IB=1,NBELEM *--------------Calcul de la normale a l'élément IREF1 = (IDIM + 1)*(MELEME.NUM(1, IB) - 1) * IREF2 = (IDIM + 1)*(MELEME.NUM(2, IB) - 1) * bp : les EF de coque epaisse etant quadratiques (coq6 et coq8), on * prend les noeuds "coins" pour eviter pb avec les noeuds 1,2,3 si courbures IREF2 = (IDIM + 1)*(MELEME.NUM(3, IB) - 1) DO 8 IC = 1, IDIM 8 XU(IC) = XCOOR(IREF2+IC)-XCOOR(IREF1+IC) XNORU = 0. DO 11 IC = 1, IDIM 11 XNORU = XNORU + XU(IC)*XU(IC) XNORU = SQRT(XNORU) DO 12 IC = 1, IDIM 12 XU(IC) = XU(IC)/XNORU IF (IDIM .EQ. 2) THEN XW(1) = -XU(2) XW(2) = XU(1) ELSE * IN = 3 IN = 5 13 IREF3 = (IDIM + 1)*(MELEME.NUM(IN, IB) - 1) DO 14 IC = 1, IDIM 14 XV(IC) = XCOOR(IREF3 + IC) - XCOOR(IREF1 + IC) XNORV = 0. DO 15 IC = 1, IDIM 15 XNORV = XNORV + XV(IC)*XV(IC) XNORV = SQRT(XNORV) DO 16 IC = 1, IDIM 16 XV(IC) = XV(IC)/XNORV XW(1) = XU(2)*XV(3) - XU(3)*XV(2) XW(2) = XU(3)*XV(1) - XU(1)*XV(3) XW(3) = XU(1)*XV(2) - XU(2)*XV(1) XNORW = 0. DO 17 IC = 1, IDIM 17 XNORW = XNORW + XW(IC)*XW(IC) XNORW = SQRT(XNORW) IF (XNORW .LT. 1.E-4) THEN if(IN.LT.NBNN) then IN = IN + 1 GOTO 13 else write(6,*) 'Difficultes pour etablir la normale de' write(6,*) 'l element',IB write(6,*) 'Verifiez votre maillage' GOTO 9990 endif ENDIF DO 18 IC = 1, IDIM 18 XW(IC) = XW(IC)/XNORW ENDIF *--------------Fin du calcul de la normale a l'élément IF (INORM .EQ. 0) THEN CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C CALL CQ8LOC(XE,NBN1,MINTE1.SHPTOT,TXR,IRR) ENDIF C DO 5002 IC=1,NBPGAU IF (INORM .EQ. 1) THEN DJAC = 1. ELSE E=DZEGAU(IC) CALL COQ8JC(IC,NBN1,E,XE,TH,TXR,SHPTOT,XJ,DJAC,IRR) ENDIF MPTVAL=IVAJAC MELVAL = IVAL(1) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC, MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(1) MELVAL = IVAL(2) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC, MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(2) IF (IDIM .EQ. 3) THEN MELVAL = IVAL(3) IBMN=MIN(IB, MELVAL.VELCHE(/2)) IGMN=MIN(IC, MELVAL.VELCHE(/1)) MELVAL.VELCHE(IGMN,IBMN)=DJAC*XW(3) ENDIF 5002 CONTINUE 5000 CONTINUE GOTO 520 ENDIF C--------------------------------------------------------------------- C DESACTIVATION DES SEGMENTS PROPRES A LA ZONE GEOMETRIQUE ISOUS C--------------------------------------------------------------------- C 520 CONTINUE MPTVAL=IVAJAC SEGSUP MPTVAL * SEGSUP INFO C SEGSUP TRA C 500 CONTINUE RETURN C C 9990 CONTINUE * C------------------------------------------------------------------- C ERREUR DANS UNE ZONE , DESACTIVATION ET RETOUR C------------------------------------------------------------------- * * * MPTVAL=IVAJAC SEGSUP MPTVAL * * * * SEGSUP INFO END